基于粗糙海面的最小二乘残差法变深度缆接收点鬼波压制技术*

李志鹏 何兵寿 杨佳佳 李凯瑞

(1. 中国海洋大学 山东青岛 266100； 2. 青岛海洋科学与技术国家实验室海洋矿产资源评价与探测技术功能实验室 山东青岛 266071；3. 海底科学与探测技术教育部重点实验室 山东青岛 266100； 4. 中国地质调查局青岛海洋地质研究所 山东青岛 266100)

海上拖缆采集的地震数据中，不仅包含来自海底反射界面的上行波场，也包含上行波场传播至海面再由海面反射向下传播的下行鬼波波场。鬼波的存在会导致明显的陷波效应，使地震数据的频带变窄，分辨率降低[1]；同时，连续发育的鬼波也会在地震剖面上产生虚假同相轴，给地震地质解释造成困扰[2]。

目前，业界进行鬼波压制的思路主要有2种：一是在野外采集中采用合适的震源组合、电缆形态或电缆组合减小鬼波的影响[3-9]；二是在资料处理过程中采用合适的算法从实测数据中消除鬼波[10-15]。第一种思路的常用方法包括：①通过立体枪阵的延时激发[3-4]压制震源鬼波；②采用双检采集[5]、上下缆采集[6]和变深度缆采集[7-9]等方式压制接收点鬼波。上述方法均在实际生产中取得了效果[16]，但实际采集中往往存在电缆垂向定位难度大、采集成本高[17]等缺陷，而且受采集设备的严格限制[18]，仅凭野外采集手段无法实现鬼波的完全压制。第二种思路的常用方法包括：①对于双检采集数据，首先进行水、陆检数据的一致性处理，然后对两种检波器接收到的信息进行特征项匹配，进而进行信息合并[10]，以达到消除鬼波的目的。②对上下缆数据，采用相位校正算法[11]或移位和减去算法[12]进行相位修正，然后加权合并去除鬼波。③对于变深度缆数据，常采用联合反褶积算法[13-14]或最小二乘反演算法[15]进行鬼波压制。

上述鬼波压制方法一般以海面水平为假设前提，此时海面的反射系数为-1；而事实上，受风速、风向、涌浪、风浪和人为活动(如航船等)[19]等因素的影响，海面往往是起伏的，此时海面的反射系数不一定是-1[20]。这种理论假设与实际情况的不相符会导致上述技术难以取得理想效果，因此必须研究基于起伏海面的鬼波压制技术。本文提出了一种基于粗糙海面的最小二乘残差法变深度缆接收点鬼波压制方法，并开展了不同模型的合成数据测试，从而实现了在粗糙海面条件下对接收点鬼波波场的有效压制。

1 粗糙海面建模与反射系数计算

真实海面是粗糙且不可预测的，对于海洋数据的鬼波压制，这种海面条件会导致后期数据处理中反射系数与鬼波延迟时差等参数的估算误差，从而会造成鬼波记录时差、能量及相位的计算误差，进而降低鬼波压制效果。理论上，基于起伏海面的鬼波压制技术可以消除上述误差，而实现起伏海面鬼波压制的前提是海面起伏量在激发、接收时刻的求取与海面成像，而事实上，受成本、工期和设备等因素的限制，野外采集时往往不对海面起伏量进行实时的密集测量，因此海面起伏的建模需要在室内进行。

一般来说，某个时刻的海面是在一定的高差范围内随机起伏的，且起伏量满足高斯分布[21]，因此可利用统计学原理对粗糙海面进行建模。

1.1 描述粗糙海面的统计参量

依据蒙特卡洛建模方法[22-24]，描述粗糙海面的统计参量主要包括均方根高度σ、相关度函数C、相关长度l、功率谱密度W、均方根斜率δ和曲率半径Rc，各参数的定义如下。

均方根高度σ是描述粗糙海面起伏程度的重要参数，其值越高，海面起伏程度越大。在拖缆地震勘探中可通过浮标数据或其它观测资料近似求取。在野外采集过程的时间跨度内，该值可视为一常数，若浮标等资料缺失，也可以通过最优值搜索[25]预估其值。其定义为

(1)

相关度函数C定义为

C(r)=〈h(x),h(x+r)〉

(2)

式(2)中：r是海面上任意两点之间的距离；〈〉是函数相关运算符。相关度函数是描述海面上任意两点高程相关程度的物理量，当r=0时，C(r)最大且有C(r)=σ2。

利用海面均方根高度将相关度函数归一化,可得到相关系数ρ

(3)

ρ(r)=1/e时所对应的r值即为相关长度，记为l。相关长度是判断粗糙海面上任意两点是否相互独立的标准之一。

粗糙海面起伏高度的功率谱密度W(k)是相关度的傅里叶变换结果[26],表示信号功率在不同频率的分布情况，是单位频率的平均功率量纲,即

(4)

功率谱密度与均方根高度存在以下关系[26]：

(5)

均方根斜率δ为粗糙海面上各点斜率的均方根值[20],其定义为

(6)

式(6)中：dx为x方向的采样大小;dh为相距dx两点对应的高度差值。

曲率半径R反映粗糙面上各点的弯曲程度[27],其定义为

(7)

1.2 粗糙海面成像与反射系数计算

粗糙海面可以看作是由大量不同振幅、频率、相位的谐波叠加而成的，各谐波的振幅可以视作独立的高斯随机变量[28]，由此可得海面任一点x处起伏量的表达式为

(8)

由于构成粗糙海面的各谐波的振幅是独立的高斯随机变量，且谐波函数的相位也满足正态分布[29]，因此由式(2)和式(4)可求得高斯密度的相关函数和功率谱密度函数。定义离散波数kn=2πn/L，可得到长度为L的一维粗糙表面的谐波谱函数为

(9)

式(9)中：N(0,1)代表均值为0、方差为1的标准正态分布的随机数；Δk为谱域相邻的谐波样本的空间波数差。再对F(kn)做傅里叶变换,可得到粗糙海面离散化后的函数关系为

(10)

式(10)中：xi=iΔx(i=-N/2+1,-N/2+2,…,N/2)，表示粗糙海面上第i个采样点。

图1是按上述理论的粗糙海面建模结果，其中图1a为二维建模成像结果，图1b是均方根值分别为0.1与0.2时抽取的y=20 m时的海面起伏情况。

图1 粗糙海面建模结果

在短时间尺度内，采用低阶小斜率近似方法[29-30]建立时变海面复反射系数的计算模型，并由此求出对应的粗糙海面反射系数[31]为

(11)

式(1)中：R0为理想光滑海面的反射系数，即-1；σ为粗糙海面的均方根振幅；s为一个和频率有关的常数，其值由入射到海面的地震波频率成分决定；α为地震波传播到海面时的入射角；v为海水声波速度。

2 基于粗糙海面的变深缆接收点鬼波压制方法

基于粗糙海面的变深缆接收点鬼波的压制主要分三步进行：①求取粗糙海面的反射系数；②在τ-p域利用求取出的粗糙海面反射系数进行上下行波波场分离，得到不包含鬼波能量的上行波波场；③将得到的海面上行波场延拓至预先定义的新的水平基准面，最终得到压制鬼波后的记录。其中，粗糙海面反射系数的计算方法已由式(11)给出，在此不再赘述。

2.1 基于最小二乘残差法的τ -p域波场分离

不考虑船舶反射等背景噪声，海洋拖缆地震数据所记录的波场p(xi,t)可以视为来自海底地层的上行反射波场u(xi,t)和该上行波传播到海面后再由海面反射向下传播的下行波场d(xi,t)(鬼波)之和[32]，即

p(xi,t)=u(xi,t)+d(xi,t)

(12)

如图2所示，D点在某时刻接收到的总波场p可视为深部反射波传播到D点的上行波场u1与由深部传播至A点再由A点反射至D点的下行波场d1之和。B点为D点在海面的投影，对应某炮点接收排列的第i道(i=1,2,3,…,n)。假设该时刻在B点同样存在一个上行波u2，其传播方向与u1相同，E点为B点在u1传播路径的投影，那么CE段即为u1比u2多传播的距离。

图2 接收点波场几何关系示意图

设u1到达C点与u2到达B点的时差为Δti,m，则有

(13)

式(13)中：αm为第m个出射角;pm为第m个射线参数;zi为第i道相对于粗糙海面的沉放深度。

再假设一维海面在A—C这一范围内无起伏，则u1从D点传播到C点的时间tui,m与d1从A点传播到D点的时间tdi,m相同，将此时间记为tsi,m，有

(14)

设B点的横坐标为xi，则A点与C点的横坐标分别为xi+Δxi,m与xi-Δxi,m。存在如下关系：t时刻传播到D点的上行波u1与t+tsi,m时刻传播到C点的上行波场一致，t时刻传播到D点的下行波场d1与t-tsi,m时刻从A点开始传播的下行波场一致。进一步可以得到如下关系：D点在t时刻接收到的总波场p可以视作C点在t+tsi,m时刻接收到的上行波场u1与A点在t-tsi,m时刻接收到的上行波场u3乘以海面反射系数之后求和,即

p(xi,t)=u1(xi+Δxi,m,t+tsi,m)+

Rmu3(xi-Δxi,m,t-tsi,m)

(15)

式(15)中：Rm为B点的海面反射系数。

假设海水弹性参数在横向上均匀不变，则A、C两点(即横坐标xi+Δxi,m与xi-Δxi,m处)所接收的上行波与B点接收到的上行波存在如下关系：

uB(xi,t-Δti,m)=uC(xi+Δxi,m,t)

(16)

uB(xi,t+Δti,m)=uC(xi-Δxi,m,t)

(17)

为了使公式简洁，下文中u(xi,t)均指B点处的波场uB。

联立式(15)～(17)可得

p(xi,t)=u(xi,t+tsi,m-

Δti,m)+Rmu(xi,t-tsi,m+Δti,m)

(18)

对式(18)进行τ-p变换可得p(pm,τ)，再将p(pm,τ)变换到频率域得p(pm,ω)，最后在频率域进行线性τ-p反变换得p(xi,ω)[33],即

-jω(pmxi+

exp[-jω(pmxi-Δti,m+tsi,m)]

(19)

式(19)中：U(pm,ω)为上行波场u(xi,t)进行τ-p变换后再变换到频率域的结果。

令

D1i,j=exp[-jω(pmxi+Δti,m-tsi,m)]

(20)

D2i,j=Rmexp[-jω(pmxi-Δti,m+tsi,m)]

(21)

Di,j=D1i,j+D2i,j

(22)

利用上式将式(19)写成矩阵形式,可得

P=DU

(23)

其中

通过求解线性方程组(23)即可得到B点上行波波场U。

本文采用最小二乘残差算法求解上行波波场U。将接收点波场P作为数据空间矢量，矩阵D作为模型空间，U作为解空间矢量。选择正常数α1、β1使得数据空间正交基矢量c1和解空间正交基矢量ω1的范数为1，由于矩阵的双对角化迭代过程须满足[34]：P=β1c1，DTc1=α1ω1，则可采用如下迭代算法[33]建立双对角矩阵元素与各空间正交基矢量的计算关系,即

βk+1ck+1=Dωk-αkck

(24)

αk+1ωk+1=DTck+1-βk+1ωk+1

(25)

k次迭代之后，有

(26)

其中

Wk=(ω1,ω2,…,ωk)

Ck=(c1,c2,…,ck)

式(26)中:ek+1代表单位矩阵中的第k+1列的向量。

令Lk+1=(Bk,αk+1ek+1)，结合式(26)可得

(27)

在每一步的迭代中计算残差rk=P-DUk，假定方程的近似解为

Uk=Wkyk

(28)

其中yk与残差存在如下关系，并代入双对角化关系P=β1c1，DTc1=α1ω1以及式(27)可得

(29)

2.2 海面起伏影响的消除

利用式(24)～(29)求得的海面上行波场本质是粗糙海面上各点的上行波场，需要将其校正至接收点所在的深度。本文采用波场延拓方法进行校正,采用式(30)所示的算子[36]进行波场延拓。

Pz+Δz=Pze±iθ

(30)

3 模型数据测试

3.1 粗糙海面层状模型

利用图3所示的层状模型验证本文方法的有效性，其中海面的起伏量如图4所示。利用有限差分技术对图3进行正演模拟，其中模拟子波为主频35 Hz的雷克子波，震源置于海面，记录道数为240道，道间距为12.5 m，记录长度为0.8 s。变深度缆形态如图5，沉放深度在10～20 m逐步加深，变深缆处采用亚网格技术[37]进行空间和时间的加密采样，粗糙海面的处理采用双变网格技术[38]进行。

图3 层状模型

图4 层状模型顶端一维粗糙海面

图5 变深度缆形态

图6a与图6b分别为粗糙海面和水平海面条件下正演记录。对比可看出，图6a中在接收到海底的一次反射波之后，由于海面起伏不平，一次反射波上行到海面并被粗糙海面反射后，记录中出现了大量的同相轴断裂、不连续等现象，信噪比也较低。利用本文方法进行鬼波压制和海面起伏量校正后的结果见图6c，图6d为图6b的鬼波压制结果。图6c与图6d中存在海面下行波导致的多次波干扰，本文利用鬼波压制过程中的下行波分量进行延拓模拟出多次波能量，并从记录中减去，从而只保留有效反射波的能量记录。对比表明，本文方法可以取得良好的鬼波压制效果，消除了粗糙海面反射的下行波带来的各种干扰。图6c中存在的少量窜扰噪声是由于海面的成像精度不足导致的，更为有效的海面成像方法正在研究之中。

图6 粗糙及水平海面鬼波压制前后地震记录

图7为图6a与图6c的振幅谱。从图7a中可以观察到明显的陷波效应的存在(如红色箭头所示)，而图7b所示的压制鬼波后的振幅谱，由于消除了鬼波引起的陷波效应，各道频谱的连续性增强，并且补偿了低频成分，表明本文方法可以有效压制粗糙海面的鬼波反射。

图7 去鬼波前后地震记录频谱

3.2 MarmousiII模型

采用二维MarmousiII模型数据进行试验，模型数据由有限差分计算得到，共50炮，每炮240道接收，接收缆形态与图5一致，正演所用的其他参数与方法和3.1节相同，模型上部海面的起伏量如图8所示。

图8 MarmousiII模型顶端粗糙海面

图9a与图9c分别是含鬼波炮集的逆时偏移剖面与共偏移距剖面(1 000 m偏移距)，图9b与图9d分别为采用本文方法压制鬼波后的逆时偏移剖面与共偏移距剖面(1 000 m偏移距)，图9e与图9f分别为图9a与图9b的局部放大，对比发现，进行鬼波压制前的深度偏移剖面，尤其是浅部与低速体附近分辨率较低，存在由鬼波形成的虚反射成像界面，信噪比远低于鬼波压制后的深度剖面。由于模型中部分反射层位相邻很近，这些反射层位的鬼波同相轴与相邻层位的反射波形成的同相轴位置非常接近甚至叠加在一起，因此在图9a中可见一些相邻连续出现的反射层以及图9b中原本较弱的反射界面在图9a中能量变强。由于同时存在多次波的干扰，为更好地确认鬼波是否被正确压制，从所绘制的去鬼波前后的共偏移距剖面(1 000 m偏移距)可以观察到，图9c中箭头①、④处为海底、低速体处的鬼波成像结果，但在图9d中得到了良好的压制；图9c中箭头②、⑤处分辨率很低，无法分辨出成像的层位，而鬼波压制后在图9d中相应位置能看到较为清晰的成像结果；图9d中箭头③处成像与图9c中相比很弱，因此判断此处为下行波产生的多次波能量；图9c中箭头⑥处约1 600 m和1 800 m两个成像界面之间可以观察到1 600 m处界面鬼波的成像，但在图9d中得到了压制，分辨率也得到了提高。由于海底反射层位分布较为复杂，偏移结果中多次波和鬼波的成像难以直接进行区分，但依据海底反射和其鬼波成像的相对位置关系，配合鬼波压制前后图件的对比，可对多次波和鬼波进行一个大致的区分。图9c、d中箭头⑦、⑧所指部分与图9a、b中红框处放大区域(图9e、f)中箭头标注部分是同一部分，在共偏移距剖面中可以更清晰地看到，箭头⑧所指部分明显远离上层能量较强的成像界面，因此判断为多次波；而对于箭头⑦处，相对于上层界面位置距离与海底和其鬼波的成像位置关系较为符合，同时其成像能量的横向强弱变化也与相邻上层的变化一致，而多次波则不应该存在这样的关系，因此判断箭头⑦处是鬼波成像。

图9 去鬼波前后地震剖面

4 结论

1) 本文方法使用蒙特卡洛原理对粗糙海面建模，并求取对应的海面反射系数，更符合海洋拖缆采集时的海面条件，具有实际应用价值。

2) 本文方法使用数学原理上更稳定、更精确、收敛速度更快的最小二乘残差算法，在τ-p域对变深度缆数据进行波场分离，从而有效压制鬼波，消除地震数据频谱中的陷波效应。

3) 不同模型的合成数据测试结果表明，本文方法能够在粗糙海面条件下对接收点鬼波波场实现准确、有效压制。