2026, Vol. 56
2. 海洋油气开发与安全保障教育部工程研究中心，山东 青岛 266100;
3. 青岛海洋科技中心海洋矿产资源评价与探测技术功能实验室，山东 青岛 266237
三分量检波器可以接收更为丰富的地震波场信息，但三分量检波器的方位对采集到的地震记录会产生很大的影响；实际采集过程中三分量检波器的放置姿态往往与期望不同，导致三个分量记录的地震信号相互泄露与耦合，严重影响后续数据处理；因此需要对三分量检波器进行重定向处理，将检波器校正至设计方向[1]。
对于海底节点(Ocean bottom nodes, OBN)勘探而言，投放时仪器容易受自由落体、海底洋流、海底底质及海底地形等因素的影响，三个分量的方向可能与设计方向有较大偏差。通常，OBN设备内置磁力计或倾角仪，可实时记录检波器的姿态参数，在数据切割时通过这些参数对其进行重定向处理，从而消除各分量之间由于姿态偏转导致的多分量波场泄露与耦合。然而即使存在记录定向参数的传感器，实际生产中还需在后续处理时采用数据驱动的方法验证其准确性以及精度，尤其是对于精度较差的检波器，必须采用数据驱动的方法对其进行准确重定向处理。
目前业内已经发展形成了众多利用被动源地震技术来确定台站中的水平分量的方位角，它们主要基于地震波的偏振特性，利用面波[2-6]、体波[7-9]、接收函数[10-11]以及环境噪声[12-14]。而对于主动源地震勘探而言，最常用的方法则是通过分析震源(气枪、电火花)激发的声波信号的粒子运动来估计重定向参数。Li等[15]分析了不同海底介质条件对初至波偏振特性的影响，提出了一种适用于有万向节和无万向节系统的检波器重定向方法，成功对北海海底地震仪进行了重定向，误差控制在5°以内。周建新等[16]提出基于三分量检波器X、Y分量的能量进行方位角误差校正，尤其在采用纵波震源时，X分量的SV波能量较大，Y分量最小。刘一峰等[17]提出一种基于振幅法的水平分量方位校正方法，解决了传统能量法中通过正切计算方位角时的困难。李鹏[18]对海底电缆进行水平面校正，并通过X、Y分量的线性关系，推算其方位角。张文波等[19]通过偏振分析获得P波的偏振方向，再结合炮检点坐标计算传播方位角，最终利用直达波偏振方向与传播方向的差异对检波器进行重定向。Jeong等[20]提出了一种使用无监督机器学习的三分量检波器重定向质控流程，以识别存在定向问题的检波器，随后通过数据驱动的方法对其进行重定向处理。闫凯鑫等[21]通过欧拉角描述坐标系旋转，并结合折射波偏振矢量的特性，建立了检波器姿态的误差函数，从而实现了基于折射波偏振特性的OBN检波器重定向。张文波等[22]利用OBN三分量记录上的初至波和反射波能量，对三分量检波器进行重定向。刘劲松等[23]提出了一种利用OBS初至波偏振方位进行定位定向的线性反演和非线性反演方法，并通过实际观测数据验证了非线性反演方法的稳定性与精度优势，为海底地震仪定位定向提供了有效的技术手段。
当前，基于直达波偏振特性的重定向方法应用广泛，但其仍存在一些问题。直达波的偏振方向与入射角及海底速度有关，且偏振方向与传播方向之间存在一定的偏差，该偏差随着偏移距的增大而增大，而基于直达波偏振特性的重定向方法认为两者是一致的，导致重定向结果产生一定的误差；另外在浅海环境中直达波初至较少，方位分布不均匀，此时利用直达波进行重定向会导致重定向结果出现偏差。
综上所述，本文在前人研究的基础上，通过将给定时窗内初至波(直达波、折射波)能量、偏振特性以及P、Z分量相关性有机结合构建了关于检波器定向参数的目标函数，适用于深水和浅水等不同水深环境。采用基于领导者机制的哈里斯鹰优化算法(Leader harris hawks optimization, LHHO)反演角度参数，显著加快了计算过程。通过相应的评价方法，进一步证明该方法的正确性与高效性。
1 方法原理 1.1 三分量地震记录偏振分析偏振分析作为研究波场特征的重要手段，对于理解波的传播特性具有意义。本研究采用协方差矩阵法对三分量地震数据进行偏振分析。对于三分量地震数据X、Y、Z，可以构造如下协方差矩阵：
| $ \boldsymbol{M}=\left(\begin{array}{ccc} \operatorname{Var}(\boldsymbol{X}) & \operatorname{Cov}(\boldsymbol{X}, \boldsymbol{X}) & \operatorname{Cov}(\boldsymbol{X}, \boldsymbol{Z}) \\ \operatorname{Cov}(\boldsymbol{Y}, \boldsymbol{X}) & \operatorname{Var}(\boldsymbol{Y}) & \operatorname{Cov}(\boldsymbol{Y}, \boldsymbol{Z}) \\ \operatorname{Cov}(\boldsymbol{Z}, \boldsymbol{X}) & \operatorname{Cov}(\boldsymbol{Z}, \boldsymbol{Y}) & \operatorname{Var}(\boldsymbol{Z}) \end{array}\right) 。$ | (1) |
在矩阵中，Var和Cov分别代表方差和协方差，其计算公式表示为：
| $ \operatorname{Cov}(\boldsymbol{X}, \boldsymbol{Y})=\frac{1}{N} \sum\limits_{i=-L}^L\left(X_i-\mu_x\right)\left(Y_i-\mu_y\right), $ | (2) |
| $ \operatorname{Var}(\boldsymbol{X})=\operatorname{Cov}(\boldsymbol{X}, \boldsymbol{X}), $ | (3) |
| $ \mu_x=\frac{1}{N} \sum\limits_{i=1}^N X_i 。$ | (4) |
式中：L=N-1/2表示时窗长度的一半，N为时窗内的采样点数；μ为数据分析窗口内的算术平均值。协方差矩阵M是一个实对称矩阵，具有非负特征值，同时矩阵的特征向量Vj和特征值λj满足如下关系：
| $ \boldsymbol{M} \boldsymbol{V}_j=\lambda_j \boldsymbol{V}_j \text { 。} $ | (5) |
对协方差矩阵M进行特征值分解，求得特征值$\lambda_1 \geqslant \lambda_2 \geqslant \lambda_3$及其相应的特征向量V1, V2, V3，其中特征向量对应特征值的方向，而特征值反映不同方向的能量大小。
基于协方差矩阵的偏振分析方法具有计算结果稳定的优点，但参与计算的样本数量(即时窗长度)对结果影响较大。若时窗长度过短，计算结果不稳定；若时窗长度过长，则可能导致多个记录段被平均，从而影响分辨率。因此本文参考Li提出的计算方法[24]，依据地震数据的主频计算相应的时窗长度。
| $ T W=\left\{\begin{array}{l} 200 \text { samples, if } f_c<0.5 \mathrm{~Hz} \\ \text { floor }\left(1 / f_c \times\left(f_s\right)\right) \times cycs \text { samples, otherwise } \end{array}\right. 。$ | (6) |
式中：fs表示数据采样率；fc表示地震记录的主频；cycs默认值为2。以图 1为例，选取虚线范围作为时窗长度进行分析。
|
( (a) X分量时窗选取；(b) Y分量时窗选取；(c) Z分量时窗选取。(a)Time window selection for X component; (b)Time window selection for Y component; (c)Time window selection for Z component. ) 图 1 三分量地震记录时窗的选取 Fig. 1 Selection of time window for three-component seismic record |
图 2为图 1所示三分量地震记录对应时窗内的质点振动轨迹和相应的偏振矢量。
|
图 2 时窗内质点振动轨迹(蓝)与偏振矢量(红) Fig. 2 Particle vibration trajectory (blue) and polarization vector (red) within the time window |
检波器的重定向任务首先是确定两个坐标系之间的对应关系，然后根据该关系将实际采集坐标系校正为预先设计的坐标系。在三维直角坐标系中，坐标系之间的变换可视为通过多次围绕不同坐标轴的旋转实现，旋转过程中的角度即为欧拉角。引入欧拉角可以更加直观地描述三分量检波器的重定向问题。未经重定向的三分量检波器通过三次旋转便可使X分量与测线方向平行，Y分量与测线方向垂直，Z分量指向垂直下方。
欧拉旋转每一次旋转均围绕相应的坐标轴进行，以依次绕Z、Y、X轴的旋转为例，其旋转过程如下：
| $ \left[\begin{array}{l} \boldsymbol{x} \\ \boldsymbol{y} \\ \boldsymbol{z} \end{array}\right]=\underbrace{\left[\begin{array}{ccc} 1 & 0 & 0 \\ 0 & \cos (\psi) & \sin (\psi) \\ 0 & -\sin (\psi) & \cos (\psi) \end{array}\right]}_{R_x(\psi)} \underbrace{\left[\begin{array}{ccc} \cos (\varphi) & 0 & -\sin (\varphi) \\ 0 & 1 & 0 \\ \sin (\varphi) & 0 & \cos (\varphi) \end{array}\right]}_{R_y(\varphi)} \underbrace{\left[\begin{array}{ccc} \cos (\theta) & \sin (\theta) & 0 \\ -\sin (\theta) & \cos (\theta) & 0 \\ 0 & 0 & 1 \end{array}\right]}_{R_z(\theta)}\left[\begin{array}{l} \boldsymbol{x}^{\prime} \\ \boldsymbol{y}^{\prime} \\ \boldsymbol{z}^{\prime} \end{array}\right] 。$ | (7) |
式中：$[\boldsymbol{x}, \boldsymbol{y}, \boldsymbol{z}]^{\mathrm{T}}$ 和$\left[\boldsymbol{x}^{\prime}, \boldsymbol{y}^{\prime}, \boldsymbol{z}^{\prime}\right]^{\mathrm{T}}$分别表示旋转后的数据和未旋转的数据；$\boldsymbol{R}_z, \boldsymbol{R}_y, \boldsymbol{R}_x$分别表示绕$Z、Y、X$轴的旋转；$\theta, \varphi, \psi$ 为各自绕旋转轴转动的角度。三个旋转矩阵按照从右到左的顺序进行计算。图 3展示了三个矩阵$\boldsymbol{R}_z, \boldsymbol{R}_y, \boldsymbol{R}_x$顺序旋转的过程。
|
( 彩色坐标轴是对黑色坐标轴应用不同旋转操作的结果。(a) θ=10°; (b) θ=10°, $\phi$=20°; (c) θ=10°, $\phi$=20°, ψ=90°。The colored coordinate frame is the result of applying different rotation operations to the black coordinate frame. ) 图 3 欧拉旋转过程 Fig. 3 Euler rotation process |
上述过程模拟检波器姿态发生变化时的旋转顺序，进行重定向时需要沿相反方向旋转得到校正后的地震记录。由于旋转矩阵$\left|\boldsymbol{R}_x\right|=\left|\boldsymbol{R}_y\right|=\left|\boldsymbol{R}_z\right|=1$，因此旋转矩阵$\boldsymbol{R}_z, \boldsymbol{R}_y, \boldsymbol{R}_x$均为非奇异矩阵，存在可逆矩阵，使得：
| $ \begin{gathered} \boldsymbol{R}_z^{-1}(\theta) \boldsymbol{R}_y^{-1}(\phi) \boldsymbol{R}_x^{-1}(\psi)\left[\begin{array}{l} \boldsymbol{x} \\ \boldsymbol{y} \\ \boldsymbol{z} \end{array}\right]=\boldsymbol{R}_z^{-1}(\theta) \boldsymbol{R}_y^{-1}(\phi), \\ \boldsymbol{R}_x^{-1}(\psi) \boldsymbol{R}_x(\psi) \boldsymbol{R}_y(\phi) \boldsymbol{R}_z(\theta)\left[\begin{array}{l} \boldsymbol{x}^{\prime} \\ \boldsymbol{y}^{\prime} \\ \boldsymbol{z}^{\prime} \end{array}\right]=\left[\begin{array}{l} \boldsymbol{x}^{\prime} \\ \boldsymbol{y}^{\prime} \\ \boldsymbol{z}^{\prime} \end{array}\right] 。\end{gathered} $ | (8) |
在运用欧拉旋转进行连续旋转时，旋转顺序对结果有显著影响，因此旋转顺序是不可交换的。对于同一组欧拉角，按照不同的旋转顺序进行旋转，得到的结果通常会有所不同。
在三维空间中将一个坐标系旋转至另一个坐标系存在两种旋转方式，即欧拉角$(\psi, \phi, \theta)$ 和$(\psi+\pi, \pi- \phi, \theta+\pi)$ 具有相同的旋转效果。
1.3 基于偏振特性约束的初至波重定向优化算法在多分量地震勘探中，由于地表低速层影响，地震波接近垂直地表方向到达检波器。各向同性介质假设条件下，P波和P-SV转换波主要分布在炮点与检波点连线所在的垂直面内，而SH波则分布在垂直于该平面的位置。
若布设的检波器方位与设计方位一致，可以根据炮点和检波点的位置关系利用公式(9)将其旋转至径向-切向坐标系，如图 4所示。
| $ \left(\begin{array}{l} \boldsymbol{R} \\ \boldsymbol{T} \\ \boldsymbol{Z} \end{array}\right)=\left(\begin{array}{ccc} \cos (\alpha) & -\sin (\alpha) & 0 \\ \sin (\alpha) & \cos (\alpha) & 0 \\ 0 & 0 & 1 \end{array}\right)\left(\begin{array}{l} \boldsymbol{X} \\ \boldsymbol{Y} \\ \boldsymbol{Z} \end{array}\right) 。$ | (9) |
旋转角度α由炮检位置关系决定。
| $ \alpha=\arctan \left(\frac{y_r-y_s}{x_r-x_s}\right) 。$ | (10) |
式中：$\left(x_r, y_r\right)$为检波点坐标；$\left(x_s, y_s\right)$为炮点位置坐标。
经过水平分量旋转后，P波和P-SV转换波的能量主要集中在Z分量和R分量上，而T分量上的能量接近于零。对于实际OBN采集的多分量地震数据，尽管T分量的深部存在明显的能量，但直达P波、浅部反射波和浅层折射波的能量通常较小。依据各分量在给定时窗内能量关系构造目标函数。
| $ F(\psi, \phi, \theta)=\frac{\sum\limits_{i=1}^m \sum\limits_{t=\tau_1}^{\tau_2} \boldsymbol{T}_i^2(t)}{\sum\limits_{i=1}^m \sum\limits_{t=\tau_1}^{\tau_2} \boldsymbol{R}_i^2(t)+\sum\limits_{i=1}^m \sum\limits_{t=\tau_1}^{\tau_2} \boldsymbol{Z}_i^2(t)} 。$ | (11) |
式中：t为地震数据采样时间；τ1, τ2分别为时窗的开始时间和结束时间；i=1, 2, 3…，为道序号；R和T分别为水平分量旋转后的径向和切向分量；Z为垂向分量；m为参与计算的总道数。
单纯以最小化切向分量上P波能量作为目标函数(见式(11))求解最优重定向参数，忽略了地震记录极性的影响，会导致多解性问题的出现。具体来说，可能出现校正后的坐标系中X轴和Y轴指向设计相反方向，Y轴和Z轴指向设计相反方向以及X轴和Z轴指向设计相反方向。
本文对上述情况展开讨论，S1, S2表示两侧对称炮点位置，O表示检波点位置，虚线表示传播方向，实线表示振动方向，利用直达波或折射波的偏振特性以及P、Z分量的相关性对多解性问题进行约束。
(1) 校正后坐标系X、Y轴反向，Y、Z轴反向
在各向同性介质假设条件下，P波粒子运动沿射线传播方向，其偏振特性为线性偏振。当检波器按照设计方位布设在海底时，直达波偏振矢量在XOZ平面内的投影如图 5所示。
|
图 5 直达波偏振矢量在XOZ平面内的投影 Fig. 5 Projection of the direct wave polarization vector onto the XOZ plane |
当X、Y轴反向(见图 6(a))，对校正后的数据进行水平分量旋转，相当于将数据旋转至径向-切向的相反方向；同理，当Y、Z轴反向(见图 6(b))，进行水平分量旋转，相当于将数据旋转至切向-垂向的相反方向。在这两种情况下，均可令公式(11)达到最小值。为避免这种情况发生，加入以下限定条件(以检波点左侧炮点S1为例)：
| $ \boldsymbol{V}_x \times \boldsymbol{V}_z>0 。$ | (12) |
|
( (a) X、Y轴反向；(b) Y、Z轴反向。(a) Reversal of X and Y axes; (b) Reversal of Y and Z axes. ) 图 6 X、Y轴反向与Y、Z轴反向时直达波偏振矢量在XOZ平面的投影 Fig. 6 Projection of the direct wave polarization vector onto the XOZ plane under X、Z- and Y、Z-axis reversals |
式中，Vx和Vz为校正坐标系中对应时窗内三分量地震数据直达波偏振矢量在XOZ平面的投影分量。当校正后的坐标系满足上述条件时，即可确保校正后的X、Y轴和Y、Z轴未反向。
折射波的射线传播方向无法通过炮点和检波点的坐标求得，但其出射角由上下两层介质的速度决定。当检波器按照设计方位放置时，检波点两侧炮点产生的折射波在XOZ平面内的投影如图 7所示。对于X、Y轴反向或Y、Z轴反向(见图 8)的情况，可以利用折射波的偏振特性加入以下限定条件(以检波点左侧炮点S1为例)：
| $ \boldsymbol{V}_{x^{\prime}} \times \boldsymbol{V}_{z^{\prime}}<0 \text { 。} $ | (13) |
|
图 7 折射波偏振矢量在XOZ平面内的投影 Fig. 7 Projection of the refracted wave polarization vector onto the XOZ plane |
|
( (a) X、Y反向；(b) Y、Z轴反向。(a) Reversal of X and Y axes; (b) Reversal of Y and Z axes. ) 图 8 X、Y轴反向与Y、Z轴反向时折射波偏振矢量在XOZ平面的投影 Fig. 8 Projection of the refracted wave polarization vector onto the XOZ plane under X、Y- and Y、Z-axis reversals |
式中，$\boldsymbol{V}_{x^{\prime}}$ 和$\boldsymbol{V}_{z^{\prime}}$为校正后坐标系对应时窗内三分量地震数据折射波偏振矢量在XOZ平面的投影分量。当校正后的坐标系满足上述条件时，即可确保校正后的X、Y轴和Y、Z轴未反向。
(2) X轴和Z轴反向
针对校正后的坐标系，如果X轴和Z轴分别指向设计方向的相反方向(见图 9)，此时对校正后的数据进行水平分量旋转，相当于将数据旋转至径向-垂向的相反方向。在这种情况下，公式(11)仍然能够取得最小值，并且满足上述约束条件。由于OBN检波器P分量为压力分量，不受检波器方位的影响，P分量和Z分量的波形相似，因此可以通过计算校正后Z分量与P分量的相关系数，选择相关系数为正的组合来排除该情况。
|
( (a) 直达波偏振矢量在XOZ平面的投影；(b) 折射波偏振矢量在XOZ平面的投影。(a) Projection of the direct wave polarization vector onto the XOZ plane; (b) Projection of the refracted wave polarization vector onto the XOZ plane. ) 图 9 X、Z轴反向时直达波和折射波偏振矢量在XOZ平面的投影 Fig. 9 Projections of direct and refracted wave polarization vectors onto the XOZ plane under X- and Z-axis reversal |
本文的目标是求解最佳检波器重定向角度$(\psi, \phi, \theta)$。当地震记录极性相反时，其对应的能量保持不变，因此反演结果存在多解性。为了克服这一问题，本文将能量函数、直达波或折射波的偏振特性以及P、Z分量相关性联合约束，构建最终的目标函数如下：
| $ \begin{gathered} F_{\text {new }}(\psi, \phi, \theta)=\lambda\left(K_1-1\right)^2+\lambda\left(K_2-1\right)^2+ \\ \frac{\sum\limits_{i=1}^m \sum\limits_{t=\tau_1}^{\tau_2} \boldsymbol{T}_i^2(\boldsymbol{t})}{\sum\limits_{i=1}^m \sum\limits_{t=\tau_1}^{\tau_2} \boldsymbol{R}_i^2(t)+\sum\limits_{i=1}^m \sum\limits_{t=\tau_1}^{\tau_2} \boldsymbol{Z}_i^2(t)} \end{gathered} 。$ | (14) |
式中，K1, K2分别为初至波偏振矢量约束条件的逻辑值和P、Z分量相关性约束条件的逻辑值。若校正后的X、Y、Z分量处于正确的坐标系，则K1, K2均为1；若校正后出现X、Y轴反向的情况：K1=0, K2=1；若校正后出现Y、Z轴反向的情况：K1=0, K2=0；若校正后出现X、Z轴反向的情况：K1=1, K2=0；λ为惩罚系数，本文中默认设为100。当目标函数Fnew最小时，对应的(ψ, $\phi$, θ)即为最佳重定向角度。
综合利用初至波(直达波、折射波)的信息进行重定向计算，适用于深水和浅水环境，可以提高重定向分析结果的稳定性和可靠性。需要指出，在利用初至波偏振特性进行约束时，对于深水OBN数据，可用于直达P波重定向分析的道数较多，且直达P波为初至波，信噪比较高，可利用直达波进行约束；而对于浅水OBN数据，直达波初至较少且较难拾取，可利用折射波初至进行约束计算。
1.4 LHHO优化算法基于领导者机制的哈里斯鹰优化算法(LHHO)是Naik[26]等在哈里斯鹰优化算法[27](Harris hawks optimization, HHO)的基础上进行改进的一种算法，其由二个阶段构成：勘探阶段、开发阶段。
1.4.1 勘探阶段这一阶段模拟哈里斯鹰在全域内对猎物进行搜索，其通过两种策略进行搜索：
| $ \begin{array}{c} {\boldsymbol{X}}(t + 1) = \\ \left\{ {\begin{array}{*{20}{l}} {{{\boldsymbol{X}}_{{\rm{rand }}}}(t) - {r_1}\left| {{{\boldsymbol{X}}_{{\rm{rand }}}}(t) - 2{r_2}{\boldsymbol{X}}(t)} \right|}&{,qp_{ap}^i}\\ {\left( {{{\boldsymbol{X}}_{{\rm{rabit }}}}(t) - {{\boldsymbol{X}}_{\rm{m}}}(t)} \right) - {r_3}\left( {{\rm{LB}} + {r_4}({\rm{UB}} - {\rm{LB}})} \right)}&{,q < p_{ap}^i} \end{array}} \right. \end{array} 。$ | (15) |
| $ p_{a p}^i=\frac{\left|f\left(\boldsymbol{X}_i\right)-f\left(\boldsymbol{X}_{\text {best }}\right)\right|}{\left|f\left(\boldsymbol{X}_{\text {worst }}\right)-f\left(\boldsymbol{X}_{\text {best }}\right)\right|}, i=1, 2, \cdots, M 。$ | (16) |
| $ \boldsymbol{X}_{\mathrm{m}}(t)=\frac{1}{M} \sum\limits_{i=1}^M \boldsymbol{X}_i(t) \text { 。} $ | (17) |
式中：$\boldsymbol{X}(t+1), \boldsymbol{X}(t)$ 分别表示相应迭代进程中鹰的位置；$\boldsymbol{X}_{\mathrm{rand}}(t), \boldsymbol{X}_{\mathrm{rabbit}}(\mathrm{t}), \boldsymbol{X}_{\mathrm{m}}(t)$ 表示当前种群的随机位置、猎物的位置以及当前种群的平均位置；$\boldsymbol{X}_{\text {best }}(t)$，$\boldsymbol{X}_{\text {worst }}(t)$ 表示当前种群中鹰的最优位置和最差位置；$p_{a p}^i$ 表示适应性栖息概率；$M$ 为种群个数；$\mathrm{UB}, \mathrm{LB}$ 为取值的上下界；$r_1, r_2, r_3, r_4$ 为介于0和1之间的随机数；$f(x)$ 为适应度函数。
1.4.2 勘探与开发的转换算法依据猎物的逃逸能量实现从勘探到开发的平滑过渡：
| $ E=2 E_0\left(1-\frac{t}{T}\right) \text { 。} $ | (18) |
式中：E为猎物的逃逸能量；E0为猎物的初始能量，其取值为(-1, 1)之间的随机数；t为当前迭代次数；T为最大迭代次数，当$|E| \geqslant 1$时算法处于勘探阶段，当$|E| < 1$时算法处于开发阶段。
1.4.3 开发阶段哈里斯鹰在开发阶段依据猎物的逃逸能量和逃脱概率(r, r < 0.5表示猎物成功逃脱, r≥0.5表示猎物未成功逃脱)采取四种策略进行位置更新：
(1) 软包围($0.5 \leqslant|E| < 1, r \geqslant 0.5$)
| $ \boldsymbol{X}(t+1)=\boldsymbol{X}_{\text {rabbit }}(t)-\boldsymbol{X}(t)-E\left|J \cdot \boldsymbol{X}_{\text {rabbit }}(t)-\boldsymbol{X}(t)\right|, $ | (19) |
| $ J=2\left(1-r_5\right) \text { 。} $ | (20) |
式中：$\Delta \boldsymbol{X}(t)$ 表示第t次迭代中个体当前位置与猎物位置的差异；J表示猎物随机跳跃的强度；r5是(0, 1)内的随机数。
(2) 硬包围($|E| < 0.5, r \geqslant 0.5$)
| $ \boldsymbol{X}(t+1)=\boldsymbol{X}_{\text {rabbit }}(t)-E\left|\boldsymbol{X}_{\text {rabbit }}(t)-\boldsymbol{X}(t)\right| 。$ | (21) |
(3) 快速俯冲的软包围$(|E| \geqslant 0.5, r < 0.5)$
| $ \begin{gathered} \boldsymbol{X}(t+1)= \\ \left\{\begin{array}{l} \boldsymbol{Y}=\boldsymbol{X}_{\text {rabbit }}(t)-E\left|J \cdot \boldsymbol{X}_{\text {rabit }}(t)-\boldsymbol{X}(t)\right|, f(\boldsymbol{Y})<f(\boldsymbol{X}(t)) \\ \boldsymbol{Z}=\boldsymbol{Y}+\boldsymbol{S} \times \operatorname{LF}(D), f(\boldsymbol{Z})<f(\boldsymbol{X}(t)) \end{array}\right. \end{gathered} 。$ | (22) |
式中：S为随机的一维向量；D为问题的维度；LF为Levy飞行策略。
(4) 快速俯冲的硬包围($|E| < 0.5, r < 0.5$)
| $ \begin{gathered} \boldsymbol{X}(t+1)= \\ \left\{\begin{array}{l} \boldsymbol{Y}=\boldsymbol{X}_{\text {rabbit }}(t)-E\left|J \cdot \boldsymbol{X}_{\text {rabbit }}(\mathrm{t})-\boldsymbol{X}_m(t)\right|, f(\boldsymbol{Y})<f(\boldsymbol{X}(t)) \\ \boldsymbol{Z}=\boldsymbol{Y}+\boldsymbol{S} \times \operatorname{LF}(D), f(\boldsymbol{Z})<f(\boldsymbol{X}(t)) \end{array}\right. \end{gathered} 。$ | (23) |
为了增强算法的全局寻优能力，利用种群中表现最优的三个个体对当前个体进行变异：
| $ \begin{gathered} \boldsymbol{X}(\mathrm{mut})= \\ 2\left(1-\frac{t}{T}\right)(2 \mathrm{rand}-1)\left(2 \boldsymbol{X}_{\text {best }}^t-\left(\boldsymbol{X}_{\text {best-1 }}^t+\boldsymbol{X}_{\text {best-2 }}^t\right)\right)+ \\ \boldsymbol{X}(t)+(2 \mathrm{rand}-1)\left(\boldsymbol{X}_{\text {best }}^t-\boldsymbol{X}(t)\right)。\end{gathered} $ | (24) |
式中：$\boldsymbol{X}$（mut）为变异后的个体位置；$\boldsymbol{X}_{\text {best }}^t, \boldsymbol{X}_{\text {best－1 }}^t, \boldsymbol{X}_{\text {best－2 }}^t$表示当前种群中表现最优的三个个体的位置。
在下一次迭代中，个体位置的更新遵循以下规则：
| $ {\boldsymbol{X}}(t + 1) = \left\{ {\begin{array}{*{20}{c}} {{\boldsymbol{X}}({\rm{mut}})}&{f({\boldsymbol{X}}({\rm{mut}})) < f({\boldsymbol{X}}(t))}\\ {{\boldsymbol{X}}(t)}&{f({\boldsymbol{X}}({\rm{mut}})) \geqslant f({\boldsymbol{X}}(t))} \end{array}} \right. 。$ | (25) |
猎物位置的更新方式如下：
| $ {{\boldsymbol{X}}_{{\rm{rabbit }}}} = \left\{ {\begin{array}{*{20}{c}} {{\boldsymbol{X}}({\rm{ mut }})}&{f({\boldsymbol{X}}({\rm{ mut }})) < f\left( {{{\boldsymbol{X}}_{{\rm{rabbit }}}}} \right)}\\ {{\boldsymbol{X}}(t)}&{f({\boldsymbol{X}}(t)) < f \left( {{{\boldsymbol{X}}_{{\rm{rabbit }}}}} \right)} \end{array}} \right. 。$ | (26) |
应用LHHO算法进行重定向的流程如图 10所示。
|
图 10 LHHO算法重定向流程图 Fig. 10 Reorientation flowchart of the LHHO algorithm |
为验证检波器重定向算法的正确性，本文利用三维弹性波有限差分模拟OBN采集过程中各分量接收到的地震记录。模型为三层水平层状结构，具体模型参数见表 1。
|
表 1 模型参数 Table 1 Model Parameters |
震源参数为：雷克子波主频20 Hz，子波延迟时间200 ms，时间采样间隔为1 ms。观测系统设计如图 11所示，共有三条检波线(红色、绿色、蓝色)布设于150 m深的海底，检波线间距为50 m，每条检波线包含64个检波点，检波点间距为20 m，炮线(黑色)位于绿色检波线正上方，炮点沉放深度为10 m，炮点间距10 m，共128个炮点。
|
图 11 三维观测系统 Fig. 11 Three-dimensional observation system |
第一条检波线(红色)第32个检波器接收的三分量地震记录如图 12(a)所示。纵波和横波在三个分量上均有分布，但由于偏振方向的不同，能量强度存在差异。第二层反射界面的反射P波在X分量上的能量强度明显小于Z分量，而X分量上的反射S波能量则强于Z分量。此外，在X分量上，检波点两侧的地震记录出现了极性反转现象。
|
( (a) XYZ分量原始共检波点道集; (b) 随机旋转的共检波点道集；(c) 重定向后共检波点道集；(d) 重定向结果(c)与原始记录(a)的残差。(a) Original CRG (XYZ components); (b) Randomly rotated CRG; (c) Reoriented CRG; (d) Errors between reoriented result (c) and original record (a). ) 图 12 三分量共检波点道集 Fig. 12 Three-component common-receiver gather (CRG) |
对共检波点道集数据进行随机旋转，以模拟实际海底三分量检波器在采集过程中可能发生的姿态变化。旋转结果如图 12(b)所示，原本Y分量深部的能量较弱，经过旋转后，Y分量记录中的直达波能量较强，深部也出现了较为明显的同相轴，且在检波点两侧部分同相轴呈现极性相反现象；原来X分量检波点两侧存在明显极性反转现象，经过随机旋转后部分同相轴不再呈现极性反转现象；随机旋转后Z分量记录检波点两侧部分同相轴也呈现极性反转现象。三分量共检波点道集旋转前后的地震记录的变化，反映了姿态变化导致的波场能量泄露与耦合。
传统角度扫描法计算精度取决于步长大小，步长过大可能错过最优解，过小则增加计算负担。为了提高计算效率，本文利用LHHO优化算法对OBN节点进行重定向。给定旋转角度(37.6°，-30.8°，-53.7°)，初始种群数量为30，共进行200次迭代，由于三维空间中存在两组欧拉角经过相同的旋转顺序可以达到相同的旋转效果，故搜索空间上下界设置为：LB=[-180°, -90°, -180°]，UB=[180°, 90°, 180°]；利用本文提出的目标函数进行反演计算。各参数迭代过程如图 13所示，随着迭代次数的增加，各参量逐渐稳定，目标函数Fnew收敛至一个较小的值，最终重定向参数稳定在(36.52°，-30.89°，-53.35°)。
|
( (a) ψ迭代过程；(b) $\phi$迭代过程；(c) θ迭代过程；(d) 目标函数迭代过程。(a) ψ iterative process; (b) $\phi$ iterative process; (c) θ iterative process; (d) Objective function iterative process. ) 图 13 迭代过程 Fig. 13 Iterative process |
利用本文方法求取的重定向参数对随机旋转后的数据进行重定向，结果如图 12(c)所示，可见之前相互泄露在各分量上的能量已基本归位。重定向结果与原始记录的残差如图 12(d)所示，仍然有少量能量残留在各分量上。这是由于重定向角度与随机旋转角度之间存在一定的误差。
图 14为随机旋转后以及重定向后各道记录偏振矢量在XOZ以及XOY平面的投影，可见随机旋转后偏振矢量在XOZ平面以及XOY的投影并不对称，而重定向后偏振矢量在XOZ平面的投影呈现关于检波点基本对称的状态，在XOY平面的投影表现出随偏移距的减小与X轴的夹角逐渐增大的特征，更加契合观测系统。
|
( (a) 重定向前各道记录偏振矢量在XOZ平面的投影；(b) 重定向前各道记录偏振矢量在XOY平面的投影；(c) 重定向后各道记录偏振矢量在XOZ平面的投影；(d) 重定向后各道记录偏振矢量在XOY平面的投影。(a) Projection of polarization vectors of each trace onto the XOZ plane before reorientation; (b) Projection of polarization vectors of each trace onto the XOY plane before reorientation; (c) Projection of polarization vectors of each trace onto the XOZ plane after reorientation; (d) Projection of polarization vectors of each trace onto the XOY plane after reorientation. ) 图 14 重定向前后各道记录偏振矢量在XOZ平面以及XOY平面的投影 Fig. 14 Projections of polarization vectors of each trace before and after reorientation onto the XOZ and XOY planes |
实际资料采集中不可避免会记录各种噪声，为了检验算法的抗噪性，对模拟数据加入不同程度的随机噪声进行测试，分析方法的抗噪能力。试验中含噪数据的信噪比分别为：10、5、2、1、0.5。控制相同参数设置(旋转角度：(37.6°，-20.8°，-53.7°)，加入随机噪声后的试算结果如表 2所示，可以看出绕X轴以及Y轴的旋转角度受信噪比的影响较小，而绕Z轴的旋转角度，随着信噪比的降低，误差逐渐增大，总体而言算法具有较好的抗噪性。
|
|
表 2 不同信噪比误差测试 Table 2 Testing under different signal-to-noise ratio conditions |
本文应用的实测资料是胜利油田某浅海三维OBN资料，水深范围在3~24 m，平均水深15 m左右，该资料已经基于OBN自带姿态传感器测定参数完成现场重定向处理，故使用本文方法对其进行剩余重定向计算，并对该重定向参数的准确性进行验证。
图 15为原始X、Y、Z三分量共检波点道集，已经过现场重定向处理，利用本文方法进行剩余重定向分析，通过炮检点坐标计算方位角，根据地震记录主频选择相应的时窗长度，利用LHHO优化算法求解剩余重定向参数，计算结果如表 3所示。利用估算的三分量检波器剩余重定向参数对OBN数据进行剩余重定向处理(见图 15)，对比剩余重定向前(见图 16)后(见图 15)地震记录，结合表 3参数，说明本文方法对剩余重定向具有一定的效果，同时也能验证原始的仪器重定向效果。
|
( (a) Z分量；(b) Y分量；(c) X分量。(a) Z component; (b) Y component; (c) X component. ) 图 15 剩余重定向后的共检波点道集 Fig. 15 Residual reorientation CRG |
|
|
表 3 剩余重定向参数 Table 3 Residual reorientation parameters |
|
( (a) Z分量；(b) Y分量；(c) X分量。(a) Z component; (b) Y component; (c) X component. ) 图 16 原始共检波点道集 Fig. 16 Original CRG |
为了更好地展示本文方法效果，假设部分检波器姿态传感器参数失灵(对其进行随机旋转)然后进行重定向分析。图 17为对剩余重定向后的数据(见图 17)进行随机旋转后的某共检波点道集，可以看到Z分量上包含了大量的横波信息(图中红框所示)，纵横波相互混叠；X分量检波点两侧部分同相轴极性并不相反。
|
( (a) Z分量；(b) Y分量；(c) X分量。(a) Z component; (b) Y component; (c) X component. ) 图 17 随机旋转后共检波点道集 Fig. 17 CRG after random rotation |
利用本文方法对其进行重定向分析，计算结果如表 4所示，根据计算得到的重定向参数进行重定向处理，结果如图 18所示，Z分量上的横波能量明显衰减(图中红框所示)，纵波能量得到一定的增强；Y分量上小炮检距处能量较强，大炮检距处能量较弱。
|
|
表 4 重定向参数 Table 4 Reorientation parameters |
|
( (a) Z分量；(b) Y分量；(c) X分量。(a) Z component; (b) Y component; (c) X component. ) 图 18 重定向后共检波点道集 Fig. 18 CRG after reorientation |
从随机旋转后各道记录偏振矢量在XOZ平面(见图 19(a))和XOY平面(见图 19(b))的投影可以看出，在XOZ平面上检波点两侧的偏振矢量投影并不关于检波点呈对称状态，在XOY平面的投影与X轴夹角并不随偏移距增大而减小。重定向后各道记录偏振矢量在XOZ平面(见图 19(c))在检波点两侧呈对称状态，在XOY平面(见图 19(d))的投影表现出随偏移距的增大与X轴的夹角逐渐减小。重定向后的偏振矢量与观测系统更为契合。
|
( 黑色箭头代表检波点的位置。(a) 重定向前各道记录偏振矢量在XOZ平面的投影；(b) 重定向前各道记录偏振矢量在XOY平面的投影；(c) 重定向后各道记录偏振矢量在XOZ平面的投影；(d) 重定向后各道记录偏振矢量在XOY平面的投影。The black arrows indicate the positions of the receiver. (a) Projection of polarization vectors of each trace onto the XOZ plane before reorientation; (b) Projection of polarization vectors of each trace onto the XOY plane before reorientation; (c) Projection of polarization vectors of each trace onto the XOZ plane after reorientation; (d) Projection of polarization vectors of each trace onto the XOY plane after reorientation. ) 图 19 重定向前后各道记录偏振矢量在XOZ平面以及XOY平面的投影。 Fig. 19 Projections of polarization vectors of each trace before and after reorientation onto the XOZ and XOY planes |
利用炮点与检波点坐标求取传播方位角，对重定向前后的数据进行水平分量旋转得到径向-切向分量，如图 20、图 21所示。重定向后的共检波点道集T分量上初至P波和浅层反射波的能量明显减弱(图中红框所示)，其主要集中在Z分量和R分量上；反射P波的能量主要集中在Z分量上，P-SV转换波的能量主要集中在R分量上。重定向处理后由于检波器姿态偏转导致的波场泄露问题得到有效地解决。
|
( (a) Z分量；(b) R分量；(c) T分量。(a) Z component; (b) R component; (c) T component. ) 图 20 随机旋转后共检波点道集 Fig. 20 CRG after random rotation |
|
( (a) Z分量；(b) R分量；(c) T分量。(a) Z component; (b) R component; (c) T component. ) 图 21 重定向后共检波点道集 Fig. 21 CRG after reorientation |
在处理过程中，为了检验重定向的效果，抽取重定向前后的Z分量共偏移距(760~780 m)道集如图 22、图 23所示，可以看到重定向前Z分量共偏移距道集(见图 22)波场一致性较差，这是由于检波器姿态发生偏离导致的；重定向后Z分量共偏移距道集(见图 23)波场一致性得到明显改善，同相轴连续性更好(如图中红框所示)。
|
图 22 随机旋转后共偏移距道集(Z分量) Fig. 22 Common-offset gather after randomly rotation (Z component) |
|
图 23 重定向后共偏移距道集(Z分量) Fig. 23 Common-offset gather after reorientation (Z component) |
三分量检波器重定向工作是OBN地震资料处理的基础。在浅海和高速海底环境中，OBN资料直达波初至数据很少，因此仅依赖直达波进行重定向工作可能会导致较大的误差，故可以利用较为发育的折射波初至进行重定向工作。此外，单纯通过最小化切向分量上的P波能量进行检波器定向时，忽略了地震记录极性的影响，容易导致多解性问题。因此，本文综合分析了直达波或折射波偏振矢量在XOZ平面上的投影，以及P、Z分量的相关性，重新构建了目标函数，并将LHHO优化算法应用于重定向角度的估计。模拟数据和实际资料测试结果表明该方法具有较高的精度和实用性。
致谢: 感谢中国海洋大学宋鹏、谭军等老师对本文研究的大力支持与帮助！感谢中石化胜利油田物探研究院对本文研究的大力支持与帮助！感谢油气勘探计算机软件国家工程研究中心提供的GeoEast地震数据处理解释一体化软件系统的支持！
| [1] |
吴志强, 张训华, 赵维娜, 等. 海底节点(OBN)地震勘探: 进展与成果[J]. 地球物理学进展, 2021, 36(1): 412-424. Wu Z Q, Zhang X H, Zhao W N, et al. Ocean bottom station nodes (OBN): Progress and achievement[J]. Progress in Geophysics, 2021, 36(1): 412-424. (  0) |
| [2] |
Laske G. Global observation of off-great-circle propagation of long-period surface waves[J]. Geophysical Journal International, 1995, 123(1): 245-259. DOI:10.1111/j.1365-246X.1995.tb06673.x ( 0) |
| [3] |
Larson E W, Ekström G. Determining surface wave arrival angle anomalies[J]. Journal of Geophysical Research: Solid Earth, 2002, 107(B6): ESE 7-1-ESE 7-15. ( 0) |
| [4] |
Stachnik J, Sheehan A F, Zietlow D W, et al. Determination of new zealand ocean bottom seismometer orientation via rayleigh-wave polarization[J]. Seismological Research Letters, 2012, 83(4): 704-713. DOI:10.1785/0220110128 ( 0) |
| [5] |
Doran A K, Laske G. Ocean-bottom seismometer instrument orientations via automated rayleigh-wave arrival-angle measurements[J]. Bulletin of the Seismological Society of America, 2017, 107(2): 691-708. DOI:10.1785/0120160165 ( 0) |
| [6] |
Scholz J R, Barruol G, Fontaine F R, et al. Orienting ocean-bottom seismometers from P-wave and Rayleigh wave polarizations[J]. Geophysical Journal International, 2016, 208(3): 1277-1289. ( 0) |
| [7] |
Schulte-Pelkum V, Masters G, Shearer P. Upper mantle anisotropy from long-period P polarization[J]. Journal of Geophysical Research: Solid Earth, 2001, 106(B10): 21917-21934. DOI:10.1029/2001JB000346 ( 0) |
| [8] |
Fontaine F R, Barruol G, Kennett B L, et al. Upper mantle anisotropy beneath Australia and Tahiti from P wave polarization: Implications for real-time earthquake location[J]. Journal of Geophysical Research: Solid Earth, 2009, 114(B3): B03306. ( 0) |
| [9] |
Niu F, Li J. Component azimuths of the CEArray stations estimated from P-wave particle motion[J]. Earthquake Science, 2011, 24(1): 3-13. DOI:10.1007/s11589-011-0764-8 ( 0) |
| [10] |
Janiszewski H A, Abers G A. Imaging the plate interface in the cascadia seismogenic zone: New constraints from offshore receiver functions[J]. Seismological Research Letters, 2015, 86(5): 1261-1269. DOI:10.1785/0220150104 ( 0) |
| [11] |
Lim H, Kim Y H, Song T R A, et al. Measurement of seismometer orientation using the tangential P-wave receiver function based on harmonic decomposition[J]. Geophysical Journal International, 2018, 212(3): 1747-1765. DOI:10.1093/gji/ggx515 ( 0) |
| [12] |
Zha Y, Webb S C, Menke W. Determining the orientations of ocean bottom seismometers using ambient noise correlation[J]. Geophysical Research Letters, 2013, 40(14): 3585-3590. DOI:10.1002/grl.50698 ( 0) |
| [13] |
Xu H, Luo Y, Tang C C, et al. Systemic comparison of seismometer horizontal orientations based on teleseismic earthquakes and ambient-noise data[J]. Bulletin of the Seismological Society of America, 2018, 108(6): 3576-3589. ( 0) |
| [14] |
Ensing J X, van Wijk K. Estimating the orientation of borehole seismometers from ambient seismic noise[J]. Bulletin of the Seismological Society of America, 2019, 109(1): 424-432. DOI:10.1785/0120180118 ( 0) |
| [15] |
Li X Y, Yuan J. Geophone orientation and coupling in three-component sea-floor data: A case study[J]. Geophysical Prospecting, 1999, 47(6): 995-1013. DOI:10.1046/j.1365-2478.1999.00160.x ( 0) |
| [16] |
周建新, 姚姚. 海上多波地震勘探检波器方位误差校正[J]. 中国海上油气(地质), 1999, 13(5): 355-358. Zhou J X, Yao Y. Orientation correction for geophone in offshore multi-component seismic exploration[J]. China offshore Oil and Gas (Geology), 1999, 13(5): 355-358. ( 0) |
| [17] |
刘一峰, 傅旦丹. 基于振幅法的海上多分量地震勘探方位校正研究[J]. CT理论与应用研究, 2012, 21(2): 231-238. Liu Y F, Fu D D. Orientation correction study for horizontal components in marine multi-component seismic exploration based on amplitude method[J]. CT Theory and Applications, 2012, 21(2): 231-238. ( 0) |
| [18] |
李鹏. 海上全方位观测系统的采集设计方法技术研究[D]. 武汉: 中国地质大学(武汉), 2014. Li P. Study on the Geometries Design of Marine Full-Azimuth Acquisition[D]. Wuhan: China University of Geosciences, 2014. ( 0) |
| [19] |
张文波, 李建峰, 孙鹏远, 等. 基于直达波偏振分析的三分量检波器定向方法[J]. 石油地球物理勘探, 2017, 52(S2): 19-25. Zhang W B, Li J F, Sun P Y, et al. A three-component geophone orientation method based on direct P-wave polarization[J]. Oil Geophysical Prospecting, 2017, 52(S2): 19-25. ( 0) |
| [20] |
Jeong W, Almubarak M S, Tsingas C. Quality control for the geophone reorientation of ocean bottom seismic data using k-means clustering[J]. Geophysical Prospecting, 2021, 69(7): 1487-1502. DOI:10.1111/1365-2478.13127 ( 0) |
| [21] |
闫凯鑫, 王芷琪, 李庆春. 利用折射波的海底节点地震三分量检波器重定向方法[J]. 石油地球物理勘探, 2023, 58(2): 305-314. Yan K X, Wang Z Q, Li Q C. Three-component geophone reorientation method based on refracted waves for ocean bottom nodes seismic survey[J]. Oil Geophysical Prospecting, 2023, 58(2): 305-314. ( 0) |
| [22] |
张文波, 张少华, 孙鹏远, 等. OBN数据三分量重定向分析与质控[J]. 石油地球物理勘探, 2023, 58(3): 610-616. Zhang W B, Zhang S H, Sun P Y, et al. Three-component reorientation analysis and quality control for OBN data[J]. Oil Geophysical Prospecting, 2023, 58(3): 610-616. ( 0) |
| [23] |
刘劲松, 刘健, 李小平, 等. 利用水波初至偏振方位反演海底地震仪位置及方位[J]. 工程地球物理学报, 2024, 21(2): 305-313. Liu J S, Liu J, Li X P, et al. OBS location and orientation correction using polarization azimuth of near source water breaks[J]. Chinese Journal of Engineering Geophysics, 2024, 21(2): 305-313. ( 0) |
| [24] |
Li H, Qu K, Rong W, et al. PolarGUI: A MATLAB-based tool for polarization analysis of the three-component seismic data using different algorithms[J]. Seismological Research Letters, 2021, 92(6): 3821-3831. DOI:10.1785/0220200439 ( 0) |
| [25] |
Gaiser J E. Applications for vector coordinate systems of 3-D converted-wave data[J]. The Leading Edge, 1999, 18(11): 1290-1300. DOI:10.1190/1.1438202 ( 0) |
| [26] |
Naik M K, Panda R, Wunnava A, et al. A leader Harris hawks optimization for 2-D Masi entropy-based multilevel image thresholding[J]. Multimedia Tools and Applications, 2021, 80: 35543-35583. DOI:10.1007/s11042-020-10467-7 ( 0) |
| [27] |
Heidari A A, Mirjalili S, Faris H, et al. Harris hawks optimization: Algorithm and applications[J]. Future Generation Computer Systems, 2019, 97: 849-872. DOI:10.1016/j.future.2019.02.028 ( 0) |
2. Engineering Research Center of Marine Petroleum Development and Security Safeguard, Ministry of Education, Qingdao 266100, China;
3. Marine Mineral Resources Evaluation and Exploration Technology Functional Laboratory, Qingdao Marine Science and Technology Center, Qingdao 266237, China