CN103675910B - 一种水陆检波器地震数据标定因子反演方法 - Google Patents

Abstract

本发明提供了一种水陆检波器地震数据标定因子反演方法，该方法包括：采集水中检波器地震数据及陆地检波器地震数据进行预处理；计算水中检波器地震数据及陆地检波器地震数据的均方根振幅比、自相关函数及互相关函数的加权和；根据加权和计算相关函数的乘积及相关函数的乘积的乘积；根据相关函数的乘积及相关函数的乘积的乘积计算最大方差模方程系数，并以此计算标定因子特征方程系数；构造标定因子特征方程并求解标定因子特征方程；根据最大方差模方程系数及标定因子计算最大方差模；根据最大方差模确定最佳标定因子；根据最佳标定因子计算反射系数；绘制最佳标定因子及海底反射系数的剖面并将最佳标定因子和海底反射系数存储于地震数据道头中。

Description

一种水陆检波器地震数据标定因子反演方法

技术领域

本发明是关于石油勘探、开发及开采技术，具体的讲，是关于一种水陆检波器地震数据标定因子反演方法。

背景技术

随着地震勘探技术的发展，海上油气勘探的难度和深度越来越大，对地震资料的信噪比和分辨率要求也越来越高。利用OBC(Ocean Bottom Cable，海底电缆)进行数据采集的技术是一种联合海上和陆地地震数据采集技术。OBC使用一个固定的排列，把检波器固定在海底接收地震波，而用航船拖拽震源进行地震波激发。在OBC数据采集过程中，至少使用三艘航船：一艘震源船，用来拖拽着气枪震源排列，进行地震波激发；一艘接收船，固定不动，连接着海底电缆，接收地震波；一艘航船或者几艘航船，铺设海底电缆和回收海底电缆。OBC目前仅仅局限于水深不超过150米的海域，随着电缆和电缆回收系统的改进，OBC技术正在向更深的海域推进。三维地震可以精确地描述储层，为了得到地下介质高精度的三维图像，必须精确地知道所有炮点和接收点的位置，同时也要精确定位电缆船，把每一个检波器放置到预先设计的位置。

在OBC数据采集时，由于海底和海面都是较强的反射界面。随着一个反射地震子波从地下达到海底，海底电缆中的检波器感应并记录下这个反射地震子波。这个反射地震子波继续向上前进达到海面，受到海面的反射，然后改变方向向下传播，达到海底。海底电缆中的检波器，再一次感应并记录下这个地震子波。同时这个地震子波受到海底的反射，然后改变方向向上传播，达到海面，受到海面的反射，然后改变方向向下传播，达到海底，反射地震子波如上的循环重复进行，形成海水鸣震多次波(交混回响)，海水鸣震多次波是海上地震勘探数据中最大的噪声干扰，消除海水鸣震多次波噪声干扰，是海上地震数据处理中最为重要的步骤。

在浅水区，使用反褶积方法，可以有效消除海水鸣震多次波干扰，以恢复一次反射地震子波。但是对于海水深度超过10米的地区，反射地震子波与后续海水鸣震多次波干扰之间的时差很大，使得反褶积算法不能有效去除海水鸣震多次波干扰产生的同相轴。如果不能有效地从地震数据中去除海水鸣震多次波干扰，那么几个子波代表着一个反射界面，这样就会模糊地质断层界面。

OBC采集的数据提供了同一位置水中检波器地震数据和陆地检波器地震数据两种数据。这两种数据分别使用水中检波器和陆地检波器记录。水中检波器是一种压力检波器，记录的是地震波产生的压力变化；陆地检波器是一种质点速度检波器，记录的是质点速度变化。由于这两种检波器的记录机理不同，对于同一位置处海水鸣震多次波干扰，表现出不同特征。与水中检波器记录的海水鸣震多次波干扰相比，陆地检波器记录的海水鸣震多次波干扰表现出极性和振幅特征差异。两种检波器记录的海水鸣震多次波干扰，其相位是不同的，振幅也是不同的，且振幅相差一个与海底反射系数成比例的常数，这个常数值就是标定因子。因此，利用这种振幅和相位特征差异，可以有效消除海水鸣震多次波干扰。消除海水鸣震多次波干扰的步骤包括：(1)在每一个接收点位置处，记录水中检波器地震数据和陆地检波器地震数据两种数据；(2)利用两种检波器的传感器灵敏度(传导常数)，调整陆地检波器地震数据的振幅，以匹配水中检波器地震数据的振幅；(3)利用标定因子，标定调整振幅后的陆地检波器地震数据；(4)把标定后的陆地检波器地震数据，与对应的水中检波器地震数据相加，得到消除海水鸣震多次波干扰后的海上地震数据。这样，计算标定因子，构成了室内海上地震数据处理消除海水鸣震多次波干扰方法的根本方法和关键步骤。

在包含海水鸣震多次波的OBC数据中，由于数据的能量是由地下反射波的能量和海水鸣震多次波的能量构成。数据的自相关包含着地下反射波的自相关和海水鸣震多次波的自相关。由于地下反射波的自相关函数能量集中在自相关函数零点附近，而海水鸣震多次波的自相关能量位于远离自相关函数零点处。消除海水鸣震多次波后，由于地震数据中已经消除了多次波能量，剩下地下反射波的能量，即零点附近自相关函数值增加，而远离零点处自相关函数值减小。自相关函数方差模可以反映这种关系，因此，为了有效估算标定因子，通常使用自相关函数方差模最小作为确定标定因子的准则。

常规处理方法，多采用扫描方法计算确定标定因子。采用预先设定一个标定因子范围值和扫描步长，采用扫描方法给出一系列的标定因子值，然后计算水中检波器地震数据和陆地检波器地震数据的数据和，再对数据和计算自相关函数，由自相关函数计算最大方差模，最后由最大的最大方差模值，确定出标定因子数值。该方法需要大量的自相关计算和最大方差模计算，因此计算非常费时。

在通常情况下，水陆两种检波器的传感器灵敏度(传导常数)数值，仅仅野外采集人员知道，而室内处理人员并不知道水陆两种检波器的传感器灵敏度数值，因此不可能使用水陆检波器的传感器灵敏度数值来调整陆地检波器地震数据的振幅。另一方面，由于各个检波器之间和连接电缆之间也存在制造和与海底接触耦合的差异，使得各个检波器之间表现出的传感器灵敏度数值也不相同，采集时仅仅知道制造商提供的制造上的传感器灵敏度数值，而真正检波器接收时的传感器灵敏度数值不可能知道。

发明内容

本发明提供一种水陆检波器地震数据标定因子反演方法，以快速、准确地进行水中检波器数据和陆地检波器数据的合并处理，达到消除地震数据中海水鸣震多次波干扰的目的。

为了实现上述目的，本发明提供一种水陆检波器地震数据标定因子反演方法，所述的方法包括：

采集利用地震人工震源激发的水中检波器地震数据及陆地检波器地震数据，并对所述的水中检波器地震数据及陆地检波器地震数据进行预处理；

计算所述水中检波器地震数据及陆地检波器地震数据的均方根振幅比；

计算所述水中检波器地震数据及陆地检波器地震数据的自相关函数及互相关函数；

计算所述自相关函数及互相关函数的加权和；

根据所述的加权和计算相关函数的乘积及所述相关函数的乘积的乘积；

根据所述相关函数的乘积及所述相关函数的乘积的乘积计算最大方差模方程系数；

根据所述的最大方差模方程系数计算标定因子特征方程系数；

根据所述的标定因子特征方程系数构造标定因子特征方程并求解所述的标定因子特征方程；

根据所述的最大方差模方程系数及所述标定因子特征方程的标定因子计算最大方差模；

根据所述的最大方差模确定最佳标定因子；

根据所述的最佳标定因子计算反射系数；

绘制最佳标定因子及海底反射系数的剖面并将所述最佳标定因子和海底反射系数存储于地震数据道头中。

在一实施例中，对所述的水中检波器地震数据及陆地检波器地震数据进行预处理，包括：定义观测系统、对所述的水中检波器地震数据及陆地检波器地震数据置标签、对所述的水中检波器地震数据及陆地检波器地震数据进行分离、去噪、滤波、速度分析及叠加处理。

在一实施例中，计算所述水中检波器地震数据及陆地检波器地震数据的均方根振幅比，包括：

按照下述公式计算水中检波器与陆地检波器均方根振幅比

$W_{HG}^i=\frac{\sqrt{\underset{j=1}{\overset{N}{\Σ}}{H}_{i,j}^2}}{\sqrt{\underset{j=1}{\overset{N}{\Σ}}{G}_{i,j}^2}}$

式中，H_i,j为水中检波器地震数据；G_i,j为陆地检波器地震数据；i为时窗参数道的顺序号，i＝1,2,…,L，L为时窗参数总道数；j为时窗参数时间样点的顺序号，j＝1,2,…,N，N为时窗参数时间样点数。

在一实施例中，计算所述水中检波器地震数据及陆地检波器地震数据的自相关函数及互相关函数，包括：

按照下述公式计算水中检波器地震数据自相关函数

$B_{i,k1}^0=\underset{j=1}{\overset{N}{\Σ}}{H}_{ij}{H}_{i,j+k1}$

式中，i为时窗参数道的顺序号，i＝1,2,…,L；k1为自相关函数样点顺序号，k1＝0,1,2,…,M1，(M1+1)是自相关函数时间样点数；

按照下述公式计算水中检波器地震数据和陆地检波器地震数据互相关函数

$B_{i,k2}^1=\underset{j=1}{\overset{N}{\Σ}}(H_{i,j}{G}_{i,j+k2}+G_{i,j}{H}_{i,j+k2})$

式中，i为时窗参数道的顺序号，i＝1,2,…,L；k2为互相关函数样点顺序号，k2＝0,1,2,…,M2，(M2+1)是互相关函数时间样点数；

按照以下公式计算陆地检波器地震数据自相关函数

$B_{i,k1}^2=\underset{j=1}{\overset{N}{\Σ}}{G}_{i,j}{G}_{i,j+k1}$

式中，i为时窗参数道的顺序号，i＝1,2,…,L。

在一实施例中，计算所述自相关函数及互相关函数的加权和，包括：

按照下述公式计算水中检波器数据自相关函数加权和

$C_{k1}^0=\underset{i=1}{\overset{L}{\Σ}}{B}_{i,k1}^0,k1=0,1,2,\mathrm{...},M1;$

按照下述公式计算计算水中检波器数据和陆地检波器互相关函数加权和

$C_{k2}^1=\underset{i=1}{\overset{L}{\Σ}}{W}_{HG}^i{B}_{i,k2}^1,k2=0,1,2,\mathrm{...},M2;$

按照以下公式计算计算陆地检波器自相关函数加权和

$C_{k1}^2=\underset{i=1}{\overset{L}{\Σ}}{\left(W_{HG}^i\right)}^2{B}_{i,k1}^2.$

在一实施例中，根据所述的加权和计算相关函数的乘积及所述相关函数的乘积的乘积，包括：

按照下述公式计算所述相关函数的乘积

$D_k^n=\underset{m=0}{\overset{n}{\Σ}}{C}_k^m{C}_k^{n-m},n=0,1,2,3,4,k=0,1,2,\mathrm{...},M,$

式中，k为相关函数样点顺序号，k＝0,1,2,…,M，(M+1)是相关函数时间样点数，且 $C_k^n=0,n=3,4,k=0,1,2,\mathrm{...},M;$

按照下述公式计算所述相关函数乘积的乘积

$F_k^n=\underset{m=0}{\overset{n}{\Σ}}{D}_k^m{D}_k^{n-m},n=0,1,2,\mathrm{...},8,k=0,1,2,\mathrm{...},M,$

其中， $D_k^n=0,n=5,6,7,8,k=0,1,2,\mathrm{...},M.$

在一实施例中，根据所述相关函数的乘积及所述相关函数的乘积的乘积计算最大方差模方程系数，包括：

按照下述公式计算最大方差模方程分子的系数Pⁿ：

$P^n=\underset{k=0}{\overset{M}{\Σ}}{F}_k^n,n=0,1,2,\mathrm{...},8;$

按照下述公式计算最大方差模方程分母的系数Qⁿ：

$Q^n=\underset{m=0}{\overset{n}{\Σ}}{E}^m{E}^{n-m},n=0,1,2,\mathrm{...},8;$

其中， $E^n=\underset{k=0}{\overset{M}{\Σ}}{D}_k^n,n=0,1,2,3,4,E^n=0,n=5,6,7,8.$

在一实施例中，根据所述的最大方差模方程系数计算标定因子特征方程系数，包括：

按照下述公式计算标定因子特征方程系数T_n：

$T_n=\underset{m=1}{\overset{n+1}{\Σ}}m(P^m{Q}^{n-m+1}-Q^m{P}^{n-m+1}),n=0,1,2,\mathrm{...},15,$

其中，Pⁿ＝0，Qⁿ＝0，n＝9,10,11,…,16。

在一实施例中，根据所述的标定因子特征方程系数构造标定因子特征方程并求解所述的标定因子特征方程，包括：

按照下述公式构造所述标定因子特征方程：

T₀+T₁α+T₂α²+…+T₁₅α¹⁵＝0；

求解所述标定因子特征方程，得到标定因子α特征方程对应的15个根α_n，n＝1,2,3,…,15。

在一实施例中，根据所述的最大方差模方程系数及所述标定因子特征方程的标定因子计算最大方差模，包括：

按照下述公式计算最大方差模Varm(α)：

$Varm\left(\mathrm{\α}\right)=\frac{\underset{n=0}{\overset{8}{\Σ}}{P}^n{\mathrm{\α}}^n}{\underset{n=0}{\overset{8}{\Σ}}{Q}^n{\mathrm{\α}}^n};$

将所述的15个根α_n，n＝1,2,3,…,15，分别代入上式中，计算出15个最大方差模值Varm(α_n)Pⁿ为最大方差模方程分子的系数，Qⁿ为最大方差模方程分母的系数。

在一实施例中，根据所述的最大方差模确定最佳标定因子，包括：

确定15个所述最大方差模值Varm(α_n)中的最大值，所述最大值对应的特征方程的根为最佳标定因子α_best＝Max{Varm(α₁),Varm(α₂),…,Varm(α₁₅)}。

在一实施例中，根据所述的最佳标定因子计算反射系数，包括：

按照下述公式计算反射系数R：

$R=\frac{{\mathrm{\α}}_{best}-1}{{\mathrm{\α}}_{best}+1}.$

本发明实施例的有益效果在于，本发明可以快速、准确地进行水中检波器数据和陆地检波器数据的合并处理，消除地震数据中海水鸣震多次波干扰，并且数据计算、快捷。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例水陆检波器地震数据标定因子反演方法流程图；

图2(a)是水陆检合并共检波点道集水检数据示意图，图2(b)是水陆检合并共检波点道集陆检数据示意图，图2(c)是现有水陆检合并数据示意图，图2(d)是本发明水陆检合并共检波点道集水陆检合并数据示意图；

图3(a)是水陆检合并共检波点道集水检数据频谱示意图，图3(b)是水陆检合并共检波点道集陆检数据频谱示意图，图3(c)是现有水陆检合并数据频谱示意图，图3(d)是本发明水陆检合并共检波点道集水陆检合并数据频谱示意图；

图4(a)是水陆检合并共炮点道集水检数据示意图，图4(b)是水陆检合并共炮点道集陆检数据示意图，图4(c)是现有水陆检合并数据，图4(d)是本发明水陆检合并共炮点道集水陆检合并数据示意图；

图5(a)是水陆检合并共炮点道集水检数据频谱示意图，图5(b)是水陆检合并共炮点道集陆检数据频谱示意图，图5(c)是现有水陆检合并数据频谱示意图，图5(d)是本发明水陆检合并共炮点道集水陆检合并数据频谱示意图；

图6(a)是水陆检合并共检波点道集水检数据自相关函数示意图，图6(b)是水陆检合并共检波点道集陆检数据自相关函数示意图，图6(c)是水陆检合并共检波点道集水陆检数据互相关函数示意图，图6(d)是本发明水陆检合并共检波点道集水陆检合并数据自相关函数示意图；

图7(a)是水陆检合并叠加水检数据示意图，图7(b)是水陆检合并叠加陆检数据示意图，图7(c)是现有水陆检合并数据示意图，图7(d)是本发明水陆检合并叠加水陆检合并数据示意图；

图8(a)是水陆检合并叠加水检数据示意图，图8(b)是水陆检合并叠加陆检数据示意图，图8(c)是现有水陆检合并数据示意图，图8(d)是本发明水陆检合并叠加水陆检合并数据示意图；

图9(a)是反演标定因子示意图，图9(b)是海底反射系数示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，本发明实施例提供一种水陆检波器地震数据标定因子反演方法，包括以下步骤：

步骤101：采集利用地震人工震源激发的水中检波器地震数据及陆地检波器地震数据，并对所述的水中检波器地震数据及陆地检波器地震数据进行预处理；

步骤102：计算所述水中检波器地震数据及陆地检波器地震数据的均方根振幅比；

步骤103：计算所述水中检波器地震数据及陆地检波器地震数据的自相关函数及互相关函数；

步骤104：计算所述自相关函数及互相关函数的加权和；

步骤105：根据所述的加权和计算相关函数的乘积及所述相关函数的乘积的乘积；

步骤106：根据所述相关函数的乘积及所述相关函数的乘积的乘积计算最大方差模方程系数；

步骤107：根据所述的最大方差模方程系数计算标定因子特征方程系数；

步骤108：根据所述的标定因子特征方程系数构造标定因子特征方程并求解所述的标定因子特征方程；

步骤109：根据所述的最大方差模方程系数及所述标定因子特征方程的标定因子计算最大方差模；

步骤110：根据所述的最大方差模确定最佳标定因子；

步骤111：根据所述的最佳标定因子计算反射系数；

步骤112：绘制最佳标定因子及海底反射系数的剖面并将所述最佳标定因子和海底反射系数存储于地震数据道头中。

由图1的流程可知，本发明的水陆检波器地震数据标定因子反演方法，通过计算水中检波器数据和陆地检波器数据的自相关函数以及它们之间的互相关函数，然后建立标定因子特征方程，通过求解标定因子特征方程，确定最佳标定因子，再由最佳标定因子计算确定海底反射系数，最后对水中检波器数据和陆地检波器数据进行合并处理，可以快速、准确地进行水中检波器数据和陆地检波器数据合并处理，消除地震数据中海水鸣震多次波干扰，有效提高地震数据信噪比和分辨率，并且计算精确、快捷。

步骤101具体实施时，用地震人工震源进行激发，产生水中检波器与陆地检波器地震数据，然后采集水中检波器与陆地检波器地震数据。数据采集完成之后，还需要对采集的水中检波器地震数据及陆地检波器地震数据进行预处理，包括：定义观测系统、对所述的水中检波器地震数据及陆地检波器地震数据置标签、对所述的水中检波器地震数据及陆地检波器地震数据进行分离、去噪、滤波、速度分析及叠加处理。

步骤102具体实施时，按照下述公式计算水中检波器与陆地检波器均方根振幅比

$W_{HG}^i=\frac{\sqrt{\underset{j=1}{\overset{N}{\Σ}}{H}_{i,j}^2}}{\sqrt{\underset{j=1}{\overset{N}{\Σ}}{G}_{i,j}^2}}---\left(1\right)$

式中，H_i,j为水中检波器地震数据；G_i,j为陆地检波器地震数据；i为时窗参数道的顺序号，i＝1,2,…,L，L为时窗参数总道数；j为时窗参数时间样点的顺序号，j＝1,2,…,N，N为时窗参数时间样点数。

步骤103具体实施时，具体包括：

1)按照下述公式计算水中检波器地震数据自相关函数

$B_{i,k1}^0=\underset{j=1}{\overset{N}{\Σ}}{H}_{i,j}{H}_{i,j+k1}---\left(2\right)$

式中，i为时窗参数道的顺序号，i＝1,2,…,L；k1为自相关函数样点顺序号，k1＝0,1,2,…,M1，(M1+1)是自相关函数时间样点数。

2)按照下述公式计算水中检波器地震数据和陆地检波器地震数据互相关函数

$B_{i,k2}^1=\underset{j=1}{\overset{N}{\Σ}}(H_{i,j}{G}_{i,j+k2}+G_{i,j}{H}_{i,j+k2})---\left(3\right)$

式中，i为时窗参数道的顺序号，i＝1,2,…,L；k2为互相关函数样点顺序号，k2＝0,1,2,…,M2，(M2+1)是互相关函数时间样点数；

3)按照以下公式计算陆地检波器地震数据自相关函数

$B_{i,k1}^2=\underset{j=1}{\overset{N}{\Σ}}{G}_{i,j}{G}_{i,j+k1}---\left(4\right)$

式中，i为时窗参数道的顺序号，i＝1,2,…,L。

步骤104具体实施时，包括：

1)按照下述公式计算水中检波器数据自相关函数加权和

$C_{k1}^0=\underset{i=1}{\overset{L}{\Σ}}{B}_{i,k1}^0,k1=0,1,2,\mathrm{...},M1;---\left(5\right)$

2)按照下述公式计算计算水中检波器数据和陆地检波器互相关函数加权和

$C_{k2}^1=\underset{i=1}{\overset{L}{\Σ}}{W}_{HG}^i{B}_{i,k2}^1,k2=0,1,2,\mathrm{...},M2;---\left(6\right)$

3)按照以下公式计算计算陆地检波器自相关函数加权和

$C_{k1}^2=\underset{i=1}{\overset{L}{\Σ}}{\left(W_{HG}^i\right)}^2{B}_{i,k1}^2---\left(7\right).$

步骤105具体实施时，包括：

1)按照下述公式计算所述相关函数的乘积

$D_k^n=\underset{m=0}{\overset{n}{\Σ}}{C}_k^m{C}_k^{n-m}---\left(8\right)$

式中，n＝0,1,2,3,4，k为相关函数样点顺序号，k＝0,1,2,…,M，(M+1)是相关函数时间样点数，且

$C_k^n=0---\left(9\right)$

其中，n＝3,4，k＝0,1,2,…,M；

2)按照下述公式计算所述相关函数乘积的乘积

$F_k^n=\underset{m=0}{\overset{n}{\Σ}}{D}_k^m{D}_k^{n-m}---\left(10\right)$

$D_k^n=0---\left(11\right)$

其中，n＝5,6,7,8，k＝0,1,2,…,M。

步骤106具体实施时，包括：

1)按照下述公式计算最大方差模方程分子的系数Pⁿ：

$P^n=\underset{k=0}{\overset{M}{\Σ}}{F}_k^n---\left(12\right)$

其中，n＝0,1,2,…,8；

2)按照下述公式计算最大方差模方程分母的系数Qⁿ：

$Q^n=\underset{m=0}{\overset{n}{\Σ}}{E}^m{E}^{n-m},n=0,1,2,\mathrm{...},8;---\left(13\right)$

其中， $E^n=\underset{k=0}{\overset{M}{\Σ}}{D}_k^{n}n=0,1,2,3,4---\left(14\right)$

Eⁿ＝0，n＝5,6,7,8。 (15)

步骤107具体实施时，按照下述公式计算标定因子特征方程系数T_n：

$T_n=\underset{m=1}{\overset{n+1}{\Σ}}m(P^m{Q}^{n-m+1}-Q^m{P}^{n-m+1}),n=0,1,2,\mathrm{...},15,---\left(16\right)$

其中，Pⁿ＝0，Qⁿ＝0，n＝9,10,11,…,16。 (17)

步骤108具体实施时，包括：

1)按照下述公式构造所述标定因子特征方程：

T₀+T₁α+T₂α²+…+T₁₅α¹⁵＝0； (18)

2)求解上述标定因子特征方程，得到标定因子α特征方程对应的15个根α_n，n＝1,2,3,…,15。

步骤109具体实施时，包括：

1)按照下述公式计算最大方差模Varm(α)：

$Varm\left(\mathrm{\α}\right)=\frac{\underset{n=0}{\overset{8}{\Σ}}{P}^n{\mathrm{\α}}^n}{\underset{n=0}{\overset{8}{\Σ}}{Q}^n{\mathrm{\α}}^n};---\left(19\right)$

2)将所述的15个根α_n，n＝1,2,3,…,15，分别代入上式中，计算出15个最大方差模值Varm(α_n)；

上式中，Pⁿ为最大方差模方程分子的系数，Qⁿ为最大方差模方程分母的系数。

步骤110具体实施时，在15个最大方差模值Varm(α_n)中，最大值对应的特征方程的根就是最佳标定因子α_best。确定最佳标定因子即：

α_best＝Max{Varm(α₁),Varm(α₂),…,Varm(α₁₅)}。 (20)

步骤110具体实施时，按照下述公式计算反射系数R：

$R=\frac{{\mathrm{\α}}_{best}-1}{{\mathrm{\α}}_{best}+1}.---\left(21\right)$

在海底电缆(Ocean Bottom Cable，OBC)数据采集中，可以在同一位置同时记录水中检波器数据与陆地检波器数据两种数据。由于两种数据位于相同的位置，它们同时记录到海水鸣震产生的多次波。通过两种数据的合理组合，可以有效地消除地震数据中海水鸣震产生的多次波。水陆检波器组合的计算公式是

$S_{i,j}=H_{i,j}+{\mathrm{\αW}}_{HG}^i{G}_{i,j}---\left(22\right)$

这里，S_i,j是水陆检波器组合数据，H_i,j是水中检波器(压力检波器)数据，G_i,j是陆地检波器(质点速度检波器)数据，α是标定因子，i是时窗参数道的顺序号，j是时窗参数时间样点的顺序号，是水中检波器与陆地检波器均方根振幅比。

为了计算标定因子(scale factor)，使用最大方差模(Varimax)准则。最大方差模计算公式为

$Varm=\frac{\underset{k1=0}{\overset{M1}{\Σ}}{A}_{k1}^4}{{\left(\underset{k1=0}{\overset{M1}{\Σ}}{A}_{k1}^2\right)}^2}---\left(23\right)$

这里，A_k1是水陆检波器组合的自相关函数，M是自相关函数长度，k1是自相关函数样点顺序号。

每一道水陆检波器组合自相关函数的计算公式是

$A_{i,k1}=\underset{j=1}{\overset{N}{\Σ}}{S}_{i,j}{S}_{i,j+k1},k1=0,1,2,\mathrm{...},M---\left(24\right)$

将方程(1)代入方程(24)，然后展开，有

$\begin{array}{c}{A}_{i,k1}=\underset{j=1}{\overset{N}{\Σ}}{H}_{i,j}{H}_{i,j+k1}+{\left({\mathrm{\αW}}_{HG}^i\right)}^2\underset{j=1}{\overset{N}{\Σ}}{G}_{i,j}{G}_{i,j+k1}\\ +{\mathrm{\αW}}_{HG}^i\underset{j=1}{\overset{N}{\Σ}}(H_{i,j}{G}_{i,j+k1}+G_{i,j}{H}_{i,j+k1})\end{array}---\left(25\right)$

令

$B_{i,k1}^0=\underset{j=1}{\overset{N}{\Σ}}{H}_{i,j}{H}_{i,j+k1}---\left(2\right)$

$B_{i,k2}^1=\underset{j=1}{\overset{N}{\Σ}}(H_{i,j}{G}_{i,j+k2}+G_{i,j}{H}_{i,j+k2})---\left(3\right)$

$B_{i,k1}^2=\underset{j=1}{\overset{N}{\Σ}}{G}_{i,j}{G}_{i,j+k1},k1=0,1,2,\mathrm{...},M1---\left(4\right)$

$C_{k1}^0=\underset{i=1}{\overset{L}{\Σ}}{B}_{i,k1}^0---\left(5\right)$

$C_{k2}^1=\underset{i=1}{\overset{L}{\Σ}}{W}_{HG}^i{B}_{i,k2}^1---\left(6\right)$

$C_{k1}^2=\underset{i=1}{\overset{L}{\Σ}}{\left(W_{HG}^i\right)}^2{B}_{i,k1}^2,k=0,1,2,\mathrm{...},M---\left(7\right)$

则有

$A_k=\underset{n=0}{\overset{2}{\Σ}}{C}_k^n{\mathrm{\α}}^n,k=0,1,2,\mathrm{...},M---\left(26\right)$

这样，有

$A_k^2=\underset{n=0}{\overset{4}{\Σ}}{D}_k^n{\mathrm{\α}}^n,k=0,1,2,\mathrm{...},M---\left(27\right)$

其中，

$D_k^n=\underset{m=0}{\overset{n}{\Σ}}{C}_k^m{C}_k^{n-m},n=0,1,2,3,4,k=0,1,2,\mathrm{...},M---\left(8\right)$

且

$C_k^n=0,n=3,4,k=0,1,2,\mathrm{...},M---\left(9\right)$

这样，有

$\underset{k=0}{\overset{M}{\Σ}}{A}_k^2=\underset{n=0}{\overset{4}{\Σ}}{E}^n{\mathrm{\α}}^n,k=0,1,2,\mathrm{...},M---\left(28\right)$

$A_k^4=\underset{n=0}{\overset{8}{\Σ}}{F}_k^n{\mathrm{\α}}^n,k=0,1,2,\mathrm{...},M---\left(29\right)$

其中，

$E^n=\underset{k=0}{\overset{M}{\Σ}}{D}_k^n,n=0,1,2,3,4---\left(14\right)$

$F_k^n=\underset{m=0}{\overset{n}{\Σ}}{D}_k^m{D}_k^{n-m},n=0,1,2,\mathrm{...},8---\left(10\right)$

且

$D_k^n=0,n=5,6,7,8---\left(11\right)$

因此，有

$\underset{k=0}{\overset{M}{\Σ}}{A}_k^4=\underset{n=0}{\overset{8}{\Σ}}{P}^n{\mathrm{\α}}^n---\left(30\right)$

${\left(\underset{k=0}{\overset{M}{\Σ}}{A}_k^2\right)}^2=\underset{n=0}{\overset{8}{\Σ}}{Q}^n{\mathrm{\α}}^n---\left(31\right)$

其中

$P^n=\underset{k=0}{\overset{M}{\Σ}}{F}_k^n,n=0,1,2,\mathrm{...},8---\left(12\right)$

$Q^n=\underset{m=0}{\overset{n}{\Σ}}{E}^m{E}^{n-m},n=0,1,2,\mathrm{...},8---\left(13\right)$

且

Eⁿ＝0，n＝5,6,7,8 (15)

这样，方程(3)的最大方差模为

$Varm\left(\mathrm{\α}\right)=\frac{\underset{n=0}{\overset{8}{\Σ}}{P}^n{\mathrm{\α}}^n}{\underset{n=0}{\overset{8}{\Σ}}{Q}^n{\mathrm{\α}}^n}---\left(32\right)$

为了计算标定因子α，方程(22)对α求导数，并令其导数为0，有

T₀+T₁α+T₂α²+…+T₁₅α¹⁵＝0 (18)

这是一个关于标定因子α的15阶特征方程。其中系数T_n为

$T_n=\underset{m=1}{\overset{n+1}{\Σ}}m(P^m{Q}^{n-m+1}-Q^m{P}^{n-m+1}),n=0,1,2,\mathrm{...},15---\left(16\right)$

且

Pⁿ＝0，Qⁿ＝0，n＝9,10,11,…,16 (33)

求解方程(18)，可以得到标定因子α特征方程对应的15个根α_n，n＝1,2,3,…,15。将这15个特征方程对应的根分别代入方程(32)，可以计算出15个最大方差模值Varm(α_n)，其中最大值对应的特征方程的根就是最佳标定因子α_best。即

α_best＝Max{Varm(α₁),Varm(α₂),…,Varm(α₁₅)} (20)

确定出最佳标定因子α_best，可以计算反射系数R，即

$R=\frac{{\mathrm{\α}}_{best}-1}{{\mathrm{\α}}_{best}+1}---\left(21\right)$

反射系数R和最佳标定因子α_best参数，用于水中检波器数据与陆地检波器数据合并。为了清楚的说明本发明的有益效果，下面结合附图进行说明：

结合图2(a)、图2(b)、图2(c)及图2(d)进行水陆检合并共检波点道集数据对比，图2(a)是水陆检合并共检波点道集水检数据示意图，图2(b)是水陆检合并共检波点道集陆检数据示意图，图2(c)是现有水陆检合并数据示意图，图2(d)是本发明水陆检合并共检波点道集水陆检合并数据示意图。

结合图3(a)、图3(b)、图3(c)及图3(d)进行水陆检合并共检波点道集数据频谱对比，图3(a)是水陆检合并共检波点道集水检数据频谱示意图，图3(b)是水陆检合并共检波点道集陆检数据频谱示意图，图3(c)是现有水陆检合并数据频谱示意图，图3(d)是本发明水陆检合并共检波点道集水陆检合并数据频谱示意图。

结合图4(a)、图4(b)、图4(c)及图4(d)进行水陆检合并共炮点道集数据对比，图4(a)是水陆检合并共炮点道集水检数据示意图，图4(b)是水陆检合并共炮点道集陆检数据示意图，图4(c)是现有水陆检合并数据，图4(d)是本发明水陆检合并共炮点道集水陆检合并数据示意图。

结合图5(a)、图5(b)、图5(c)及图5(d)进行水陆检合并共炮点道集数据频谱对比，图5(a)是水陆检合并共炮点道集水检数据频谱示意图，图5(b)是水陆检合并共炮点道集陆检数据频谱示意图，图5(c)是现有水陆检合并数据频谱示意图，图5(d)是本发明水陆检合并共炮点道集水陆检合并数据频谱示意图。

结合图6(a)、图6(b)、图6(c)及图6(d)进行水陆检合并共检波点道集数据平均相关函数对比，图6(a)是水陆检合并共检波点道集水检数据自相关函数示意图，图6(b)是水陆检合并共检波点道集陆检数据自相关函数示意图，图6(c)是水陆检合并共检波点道集水陆检数据互相关函数示意图，图6(d)是本发明水陆检合并共检波点道集水陆检合并数据自相关函数示意图。

结合图7(a)、图7(b)、图7(c)及图7(d)进行水陆检合并叠加数据对比，图7(a)是水陆检合并叠加水检数据示意图，图7(b)是水陆检合并叠加陆检数据示意图，图7(c)是现有水陆检合并数据示意图，图7(d)是本发明水陆检合并叠加水陆检合并数据示意图。

结合图8(a)、图8(b)、图8(c)及图8(d)进行水陆检合并叠加数据频谱对比，图8(a)是水陆检合并叠加水检数据示意图，图8(b)是水陆检合并叠加陆检数据示意图，图8(c)是现有水陆检合并数据示意图，图8(d)是本发明水陆检合并叠加水陆检合并数据示意图。

结合图9(a)及图9(b)进行反演标定因子和海底反射系数结果，图9(a)是反演标定因子示意图，图9(b)是海底反射系数示意图。

由图2(a)至图2(d)、图3(a)至图3(d)、图4(a)至图4(d)、图5(a)至图5(d)、图6(a)至图6(d)、图7(a)至图7(d)和图8(a)至图8(d)可以看出，利用本发明反演的标定因子和海底反射系数，有效消除了水层产生的多次波干扰，提高了OBC数据信噪比。

由图3(a)至图3(d)、图5(a)至图5(d)、和图8(a)至图8(d)可以看出，利用本发明反演的标定因子和海底反射系数，有效拓宽了海底电缆数据有效频带，提高了OBC数据分辨率。因此，本发明反演的标定因子和海底反射系数，满足了实际数据处理的要求。

本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

Claims (12)

1.一种水陆检波器地震数据标定因子反演方法，其特征在于，所述的方法包括：

采集利用地震人工震源激发的水中检波器地震数据及陆地检波器地震数据，并对所述的水中检波器地震数据及陆地检波器地震数据进行预处理；

计算所述水中检波器地震数据及陆地检波器地震数据的均方根振幅比；

计算所述水中检波器地震数据及陆地检波器地震数据的自相关函数及互相关函数；

计算所述自相关函数及互相关函数的加权和；

根据所述的加权和计算相关函数的乘积及所述相关函数的乘积的乘积；

根据所述相关函数的乘积及所述相关函数的乘积的乘积计算最大方差模方程系数；

根据所述的最大方差模方程系数计算标定因子特征方程系数；

根据所述的标定因子特征方程系数构造标定因子特征方程并求解所述的标定因子特征方程；

根据所述的最大方差模方程系数及所述标定因子特征方程的标定因子计算最大方差模；

根据所述的最大方差模确定最佳标定因子；

根据所述的最佳标定因子计算反射系数；

绘制最佳标定因子及海底反射系数的剖面并将所述最佳标定因子和海底反射系数存储于地震数据道头中。

2.根据权利要求1所述的方法，其特征在于，对所述的水中检波器地震数据及陆地检波器地震数据进行预处理，包括：定义观测系统、所述的水中检波器地震数据及陆地检波器地震数据置标签、对所述的水中检波器地震数据及陆地检波器地震数据进行分离、去噪、滤波、速度分析及叠加处理。

3.根据权利要求2所述的方法，其特征在于，计算所述水中检波器地震数据及陆地检波器地震数据的均方根振幅比，包括：

按照下述公式计算水中检波器与陆地检波器均方根振幅比

$W_{HG}^i=\frac{\sqrt{\underset{j=1}{\overset{N}{\Σ}}{H}_{i,j}^2}}{\sqrt{\underset{j=1}{\overset{N}{\Σ}}{G}_{i,j}^2}}$

式中，H_i,j为水中检波器地震数据；G_i,j为陆地检波器地震数据；i为时窗参数道的顺序号，i＝1,2,…,L，L为时窗参数总道数；j为时窗参数时间样点的顺序号，j＝1,2,…,N，N为时窗参数时间样点数。

4.根据权利要求3所述的方法，其特征在于，计算所述水中检波器地震数据及陆地检波器地震数据的自相关函数及互相关函数，包括：

按照下述公式计算水中检波器地震数据自相关函数

$B_{i,k1}^0=\underset{j=1}{\overset{N}{\Σ}}{H}_{i,j}{H}_{i,j+k1}$

式中，i为时窗参数道的顺序号，i＝1,2,…,L；k1为自相关函数样点顺序号，k1＝0,1,2,…,M1，(M1+1)是自相关函数时间样点数；

按照下述公式计算水中检波器地震数据和陆地检波器地震数据互相关函数

$B_{i,k2}^1=\underset{j=1}{\overset{N}{\Σ}}(H_{i,j}{G}_{i,j+k2}+G_{i,j}{H}_{i,j+k2})$

式中，i为时窗参数道的顺序号，i＝1,2,…,L；k2为互相关函数样点顺序号，k2＝0,1,2,…,M2，(M2+1)是互相关函数时间样点数；

按照以下公式计算陆地检波器地震数据自相关函数

$B_{i,k1}^2=\underset{j=1}{\overset{N}{\Σ}}{G}_{i,j}{G}_{i,j+k1}$

式中，i为时窗参数道的顺序号，i＝1,2,…,L。

5.根据权利要求4所述的方法，其特征在于，计算所述自相关函数及互相关函数的加权和，包括：

按照下述公式计算水中检波器数据自相关函数加权和

$C_{k1}^0=\underset{i=1}{\overset{L}{\Σ}}{B}_{i,k1}^0;$

按照下述公式计算计算水中检波器数据和陆地检波器互相关函数加权和

$C_{k2}^1=\underset{i=1}{\overset{L}{\Σ}}{W}_{HG}^i{B}_{i,k2}^1;$

按照以下公式计算计算陆地检波器自相关函数加权和

$C_{k1}^2=\underset{i=1}{\overset{L}{\Σ}}{\left(W_{HG}^i\right)}^2{B}_{i,k1}^2.$

6.根据权利要求5所述的方法，其特征在于，根据所述的加权和计算相关函数的乘积及所述相关函数的乘积的乘积，包括：

按照下述公式计算所述相关函数的乘积

$D_k^n=\underset{m=0}{\overset{n}{\Σ}}{C}_k^m{C}_k^{n-m},n=0,1,2,3,4,k=0,1,2,...,M,$

式中，k为相关函数样点顺序号，k＝0,1,2,…,M，(M+1)是相关函数时间样点数，且 $C_k^n=0,n=3,4,k=0,1,2,...,M;$

按照下述公式计算所述相关函数乘积的乘积

$F_k^n=\underset{m=0}{\overset{n}{\Σ}}{D}_k^m{D}_k^{n-m},n=0,1,2,...,8,k=0,1,2,...,M,$

其中， $D_k^n=0,n=5,6,7,8,k=0,1,2,...,M.$

7.根据权利要求6所述的方法，其特征在于，根据所述相关函数的乘积及所述相关函数的乘积的乘积计算最大方差模方程系数，包括：

按照下述公式计算最大方差模方程分子的系数Pⁿ：

$P^n=\underset{k=0}{\overset{M}{\Σ}}{F}_k^n,n=0,1,2,...,8;$

按照下述公式计算最大方差模方程分母的系数Qⁿ：

$Q^n=\underset{m=0}{\overset{n}{\Σ}}{E}^m{E}^{n-m},n=0,1,2,...,8;$

其中， $E^n=\underset{k=0}{\overset{M}{\Σ}}{D}_k^n,n=0,1,2,3,4,E^n=0,n=5,6,7,8.$

8.根据权利要求7所述的方法，其特征在于，根据所述的最大方差模方程系数计算标定因子特征方程系数，包括：

按照下述公式计算标定因子特征方程系数T_n：

$T_n=\underset{m=1}{\overset{n+1}{\Σ}}m(P^m{Q}^{n-m+1}-Q^m{P}^{n-m+1}),n=0,1,2,...,15,$

其中，Pⁿ＝0，Qⁿ＝0，n＝9,10,11,…,16。

9.根据权利要求8所述的方法，其特征在于，根据所述的标定因子特征方程系数构造标定因子特征方程并求解所述的标定因子特征方程，包括：

按照下述公式构造所述标定因子特征方程：

T₀+T₁α+T₂α²+…+T₁₅α¹⁵＝0；

求解所述标定因子特征方程，得到标定因子α特征方程对应的15个根α_n，n＝1,2,3,…,15。

10.根据权利要求9所述的方法，其特征在于，根据所述的最大方差模方程系数及所述标定因子特征方程的标定因子计算最大方差模，包括：

按照下述公式计算最大方差模Varm(α)：

$Varm\left(\mathrm{\α}\right)=\frac{\underset{n=0}{\overset{8}{\Σ}}{P}^n{\mathrm{\α}}^n}{\underset{n=0}{\overset{8}{\Σ}}{Q}^n{\mathrm{\α}}^n};$

将所述的15个根α_n，n＝1,2,3,…,15，分别代入上式中，计算出15个最大方差模值Varm(α_n)Pⁿ为最大方差模方程分子的系数，Qⁿ为最大方差模方程分母的系数。

11.根据权利要求10所述的方法，其特征在于，根据所述的最大方差模确定最佳标定因子，包括：

确定15个所述最大方差模值Varm(α_n)中的最大值，所述最大值对应的特征方程的根为最佳标定因子α_best＝Max{Varm(α₁),Varm(α₂),…,Varm(α₁₅)}。

12.根据权利要求11所述的方法，其特征在于，根据所述的最佳标定因子计算反射系数，包括：

按照下述公式计算反射系数R：

$R=\frac{{\mathrm{\α}}_{best}-1}{{\mathrm{\α}}_{best}+1}.$