基于水波初至偏振方位的海底地震仪位置及方位反演方法

摘要

本发明公开了一种基于水波初至偏振方位的海底地震仪位置及方位反演方法，采用非线性反演方法，对于一个待修正的海底地震仪，利用水波初至偏振方位残差之差构造目标函数，并求取令目标函数最小的地震仪位置，在求得精确的地震仪位置之后，再以炮点方位与偏振方位之差求得地震仪的实际方位指向。尤其是以多对水波初至偏振方位残差之差的平方和作为目标函数，应用DE最优化方法求取使目标函数取最小值的地震仪位置的方案，在目标函数中消去了地震仪方位的影响，具有较快的收敛速度和较好的稳定性，相比线性方法具有较好的精度。

说明书

技术领域

本发明涉及海底地震仪的位置及方法反演方法，具体涉及一种基于水波初至偏振方位的海底地震仪位置及方位反演方法。

背景技术

OBS海底地震仪在布放到海底之后，其方位指向通常是随机的，由于海流等因素的影响，其实际位置通常也会偏离预定位置，因此需要对其位置及方位进行反演处理。传统方法采用线性反演的方式，在一些特殊情况下，存在不收敛的情况。

发明内容

针对现有技术的不足，本发明旨在提供一种基于水波初至偏振方位的海底地震仪位置及方位反演方法，采用非线性反演方法，对于一个待修正的海底地震仪，以多对水波初至偏振方位残差之差的平方和作为目标函数，应用DE最优化方法求取使目标函数取最小值的地震仪位置，在目标函数中消去了地震仪方位的影响，具有较快的收敛速度和较好的稳定性。

为了实现上述技术目的，本发明采用如下技术方案：

基于水波初至偏振方位的海底地震仪位置及方位反演方法，具体为：

S1、根据待修正的海底地震仪的设计布放位置，利用该地震仪接收N个不同位置的炮点的三分量地震波形信号；

S2、对于一个待修正的海底地震仪，利用信号的水波初至偏振方位残差构造目标函数，并求取使目标函数取最小值的待修正的海底地震仪位置；

S3、利用步骤S1求得使目标函数取最小值的待修正的海底地震仪位置计算待修正的海底地震仪位置指向炮点方向的方位角，再以该待修正的海底地震仪位置指向炮点的方位角与炮点信号的水波初至偏振方位角之差求得待修正的海底地震仪的实际方位指向，并据此得到海底地震仪方位的修正量。

进一步地，步骤S2中，以多对信号的水波初至偏振方位残差之差的平方和作为目标函数，应用DE最优化方法求取使目标函数取最小值的待修正的海底地震仪位置。

进一步地，步骤S2的具体过程为：

设待修正的海底地震仪在海底的设计布放位置为(x,y)，实际位置为(x+δx,y+δy)，(δx,δy)为偏移量；

用多个炮点对的方位残差之差的平方和构造目标函数S₂：

和别为待修正的海底地震仪位置指向第i号和第j号炮点的方位角的观测值，i≠j，表达式如下：

和分别为待修正的海底地震仪位置指向第i号和第j号炮点方向的方位角的理论值，i≠j，表达式如下：

φ₀表示待修正的海底地震仪的方位的修正量；和分别为根据第i号和第j号炮点的信号的三分量水波初至波形计算得到的偏振方位角，和分别表示待修正的海底地震仪在初始位置偏离(δx,δy)时指向第i号和第j号炮点的理论方位角；

由此可得：

采用DE非线性最优化方法求取使S₂取最小值的偏移量(δx,δy)。

作为另一种实施方式，步骤S2中，设待修正的海底地震仪在海底的设计布放位置为(x,y)，实际位置为(x+δx,y+δy)，(δx,δy)为偏移量；

构造目标函数S₁，求取使S₁取最小值的(φ₀，δx，δy)：

其中，为待修正的海底地震仪位置指向第i号炮点的方位角的观测值，表达式如下：

为待修正的海底地震仪位置指向第i号炮点方向的方位角的理论值，表达式如下：

S₁是所有炮点的方位残差的平方和，表示在初始位置偏离(δx,δy)时指向i号炮点的理论方位角，N为炮点的数量；φ₀表示待修正的海底地震仪的方位的修正量；为根据第i号炮点的信号的三分量水波初至波形计算得到的偏振方位角。

更进一步地，步骤S3具体为：

针对第i号炮点计算的方位修正量φ₀为：

对于一个待修正的海底地震仪，针对每个炮点i，1≤i≤N都可计算得到一个方位修正量φ₀，对所有结果求平均值，即得到最终的海底地震仪方位修正量φ'₀：

本发明的有益效果在于：采用非线性反演方法，对于一个待修正的海底地震仪，利用水波初至偏振方位残差之差构造目标函数，并求取令目标函数最小的地震仪位置，在求得精确的地震仪位置之后，再以炮点方位与偏振方位之差求得地震仪的实际方位指向，准确性较好。尤其是以多对水波初至偏振方位残差之差的平方和作为目标函数，应用DE最优化方法求取使目标函数取最小值的地震仪位置的方案，在目标函数中消去了地震仪方位的影响，具有较快的收敛速度和较好的稳定性。

附图说明

图1为本发明实施例3中实际数据采用第二种方法的迭代收敛情况示意图；

图2为本发明实施例3中实际数据采用第一种方法的迭代收敛情况示意图。

具体实施方式

以下将结合附图对本发明作进一步的描述，需要说明的是，本实施例以本技术方案为前提，给出了详细的实施方式和具体的操作过程，但本发明的保护范围并不限于本实施例。

实施例1

本实施例提供一种基于水波初至偏振方位的海底地震仪位置及方位反演方法，具体为：

S1、根据待修正的海底地震仪的设计布放位置，利用该地震仪接收N个不同位置的炮点的三分量地震波形信号

S2、设待修正的海底地震仪在海底的设计布放位置为(x,y)，实际位置为(x+δx,y+δy)，(δx,δy)为偏移量。

根据纵波的性质，纵波的偏振方向与传播方向相同，因此炮点放出的信号的水波初至偏振方位与待修正的海底地震仪位置指向炮点的方向一致，待修正的海底地震仪位置指向第i号炮点方向的观测值可以表示为：

其中为根据第i号炮点的三分量水波初至波形计算得到的偏振方位角，φ₀表示待修正的海底地震仪的方位的修正量；待修正的海底地震仪位置指向第i号炮点方向的理论值可由待修正的海底地震仪位置与炮点位置计算得到：

本实施例提供两种方式计算待修正的海底地震仪的方位的修正量φ₀。

1、作为一种实施方式，直接的方法是，构造目标函数S₁，求取使S₁取最小值的(φ₀，δx，δy)：

S₁是所有炮点方位残差的平方和，表示在初始位置偏离(δx,δy)时指向i号炮点的理论方位角，N为炮点的数量。

2、作为另一种实施方式，可以用多个炮点对的方位残差之差的平方和构造目标函数S₂：

将前述各式代入上式，得到：

φ₀表示待修正的海底地震仪的方位的修正量；和分别为根据第i号和第j号炮点的信号的三分量水波初至波形计算得到的偏振方位角，和分别表示待修正的海底地震仪在初始位置偏离(δx,δy)时指向第i号和第j号炮点的理论方位角；

式(5)中，对于同一待修正的海底地震仪，方位修正量相同，因而被消掉，因此只含有两个未知数。

采用DE(deferential evolution)非线性最优化方法求取使S₂取最小值的偏移量(δx,δy)，针对第i号炮点计算的方位修正量φ₀为：

对于一个待修正的海底地震仪，针对每个炮点i，1≤i≤N都可计算得到一个方位修正量φ₀，对所有结果求平均值，即得到最终的海底地震仪方位修正量φ仩₀：

实施例2

本实施例旨在提供一种如实施例1所述方法所得反演结果的误差计算方法。

本实施例采用bootstrap方法给出反演结果的误差估计，具体做法是根据给定的实际初至偏振数据和方位角误差信息，按实际的炮点/台站位置生成一组理论观测数据，并将原有数据的残差随机加到理论数据中，然后对这些理论观测数据进行反演计算，以结果的统计方差统作为误差。

对于一个海底地震仪，假定有N个炮点数据参与计算，所涉及数据包括台站水平位置及方位(lat，lon，φ₀)，炮点坐标、方位角及残差(lat_i，lon_i，θ_i，r_i)，根据这组数据生成一组“新的理论数据”的具体步骤如下：

1.生成随机序列R_i，i＝1，N。先对R_i序列赋初值R_i＝i；

2.对R_i中的每个元素进行循环，每次循环生成[1，N]区间内的一个随机数j，然后R_i和R_j的值对调。循环完成，则随机序列R_i确定。实际程序中采用C语言中的系统调用函数srand()/rand()生成随机数，并以调用时的系统时间作为seed。

3.则新的一组观测数据为

重复上述步骤生成M组“观测数据”，对每组数据进行反演得到M个结果，这些结果的统计方差即是反演结果的误差。

实施例3

本实施例旨在对实施例1所述的方法进行进一步的检验。

为检验上述方法，本实施例分别进行了理论模型试验和实际数据的处理，所用数据来自2003年在我国南海某区域进行的OBS观测，这次观测有5台OBS台站。理论模型采用其中的obs4的实际观测系统，在理论计算得到的方位角中加上均值为0并且符合高斯分布的随机扰动，随机扰动的方差为6°，(φ₀，δx，δy)各值的真实值均为0。表1是进行了1000次这种计算的统计结果，从表1中的结果可以看出，非线性反演方法的结果好于线性反演的结果，各反演量的平均值更接近真实值，方差也较小。

表2是5个台站实际数据的反演结果。从结果可以看出，两种方法得到的地震仪方位φ₀非常接近，而地震仪的位置有些台站两种方法偏差较大，达到30米左右，但结合表1理论模型的结果，这个偏差也在正常范围内。

表1

表2实际数据的反演结果

理论模型实验和实际数据的结果表明，采用公式(3)(5)的两种方法得到的结果几乎完全相同，但迭代收敛情况公式(5)明显好于公式(3)的方法。如图1－2所示，图1－2是一个实际的obs台站数据分别采用公式(5)(3)的迭代收敛情况，公式(5)在迭代到600多次时即达到最优稳定状态，而采用公式(3)的方法，在迭代1000多次之后才达到稳定最优状态。

相关参考文献：

[1]Christeson，G.L.，K.D.McIntosh，T.H.Shipley，E.R.Flueh，and H.Goedde(1999)，Structure of the costa rica convergent margin，offshore nicoyapeninsula，JOURNAL OF GEOPHYSICAL RESEARCH－SOLID EARTH，104(B11)，25，443–25，468.

[2]JURKEVICS，A.(1988)，Polarization analysis of 3－component arraydata，B.SEISMOL.SOC.AM.，78(5)，1725–1743.

[3]Li，X.Y.，andJ.X.Yuan(1999)，Geophone orientation and coupling inthree－component sea－floor data：a case study，GEOPHYS.PROSPECT.，47(6)，995–1013.

[4]NAKAMURA，Y.，P.L.DONOHO，P.H.ROPER，and P.M.MCPHERSON(1987)，Large－offset seismic surveying using ocean－bottom seismographs and air guns－instrumentation and field technique，GEOPHYSICS，52(12)，1601–1611.

[5]PARK，J.，F.L.VERNON，andC.R.LINDBERG(1987)，Frequency－dependentpolarization analysis of high－frequency seismograms，J GEOPHYS RES－SOLID，92(B12)，12，664–12，674.

[6]Ruzek，B.，and M.Kvasnicka(2001)，Differential evolution algorithm inthe earthquake hypocenter location，PURE APPL GEOPHYS，158(4)，667–693.

[7]Ruzek，B.，P.Hrubcova，M.Novotny，A.Spicak，and O.Karousova(2007)，Inversion of travel times obtained during active seismic refractionexperiments celebration 2000，alp 2002 and sudetes 2003，STUD GEOPHYS GEOD，51(1)，141–164.

对于本领域的技术人员来说，可以根据以上的技术方案和构思，给出各种相应的改变和变形，而所有的这些改变和变形，都应该包括在本发明权利要求的保护范围之内。