同频同向的连续海底地震波场波动成分分离方法

1.本发明涉及水中目标航行诱发的海底地震波场波动成分分离技术领域，具体地涉及同频同向的连续海底地震波场波动成分分离方法。


背景技术：

2.水中目标产生的信号包括脉冲信号和连续信号，二者均可在海底诱发海底地震波。
3.对于脉冲信号而言，其诱发的海底地震波场，根据波动成分出现的先后，其波动成分分离相对较易；但对于连续信号而言，由其诱发的连续海底地震波场多种波动成分，包括但不限于水声信号、scholte波、海底纵波、海底横波，混合在一起，难以分离；
4.尤其针对是同频同向的连续海底地震波信号，含多阶简正波，scholte波等波动成分，使得在时域、频域等方面难以实现多种波动成分的分离。
5.现有技术的研究，大多为针对脉冲声源，如水中气枪，诱发的海底地震波进行波动成分离；现有技术的缺陷如下：
6.1.由于现有技术尚不完善，从而造成对连续点源激发的地震波信号分离精度不高，进一步导致无法应用；
7.2.由于现有技术尚不完善，从而造成对连续点源激发的地震波信号的波动成分能量选取缺乏有效原则，影响应用工作展开。
8.但是连续声源诱发的地震波信号的波动成分分离，是后续利用水中目标航行诱发的连续海底地震波场中的scholte波、水声信号等波动成分，开展水中目标探测以及基于水声信号和地声信号联合探测等研究的基础性工作；
9.为了解决上述第二个问题，即能量选取原则等问题，目前还有极少数现有技术采用频率-波数变换，即二维傅里叶变换的形式，将传感器采集得到的时间
‑ꢀ
空间域的数据变换到频率-波数域。这种方法的缺陷在于：
10.由于该方法未考虑海底地震波场中波动成分能量的相对大小，从而导致其分辨率依然不高，以至于不可进行实际应用。
11.综上所述，当前尚无，但亟需一种有效的且具备可应用性能的连续声源激发的海底地震波信号波动成分分离方法；


技术实现要素：

12.本发明针对上述问题，提供同频同向的连续海底地震波场波动成分分离方法，其目的在于自动精准分离水中目标诱发的多种同频同向的海底地震波波动成分；极大的提高波动成分的分辨率，尤其适合于多种波动成分的波速相近的情形。
13.为解决上述问题，本发明提供的技术方案为：
14.一种同频同向的连续海底地震波场波动成分分离方法，包含以下步骤：
15.s100.获取海底地震波场的时间-空间域数据；所述时间-空间域数据为关于时间
和空间的二维函数；
16.s200.将所述时间-空间域数据转化为频率-波数域图像，具体包含以下步骤：
17.s210.对所述时间-空间域数据进行频率-波数变换操作，获得海底地震波场的频率-波数域数据；所述频率-波数域数据包含有用且需提取的有用信号和无用且需去除的干扰噪声；
18.s220.然后记录下海底地震波场的所述频率-波数域数据的实部和虚部；
19.s230.然后将所述海底地震波场的所述频率-波数域数据进行信号离散化处理，得到所述频率-波数域图像；
20.s300.将所述频率-波数域图像进行像素化处理，得到像素化后的频率-波数域图像；所述频率-波数域图像的行为频率，列为波数；所述频率-波数域图像的每个像素点的像素值的对应范围为[0-255]；
[0021]
s400.按所述频率-波数域图像的行进行阈值划分操作，计算得到所述频率
‑ꢀ
波数域图像的每一行的频率-波数域图像能量均值和图像能量标准差；
[0022]
s500.根据所述频率-波数域图像能量均值和所述图像能量标准差，将所述频率-波数域图像的每一行的图像能量阈值设为该行的所述图像能量均值与人工预设的标准差倍数乘以该行的所述图像能量标准差之和；
[0023]
s600.对所述频率-波数域图像进行阈值分割操作，具体包含以下步骤：
[0024]
s610.逐一将每一行中的每个像素点，与该像素点所在的行的图像能量阈值进行大小对比；
[0025]
s620.根据s610的对比结果，做出如下操作：
[0026]
如果该像素点的像素值比该像素点所在的行的图像能量阈值小，则将该像素点的像素值置0；
[0027]
否则，保持该像素点的像素值的值不变；
[0028]
s630.再次执行s610，直至比对完所述频率-波数域图像中所有的像素点；
[0029]
s700.对经过所述阈值分割操作后的所述频率-波数域图像进行连通性分析操作，搜索得到所述有用信号的所在的连通域和所述干扰噪声的所在的连通域；
[0030]
s800.对s700中搜索得到所述特定波动成分的所在的连通域进行频率-波数域滤波操作；然后将s220中的海底地震波场的所述频率-波数域数据的虚部补充进入所述频率-波数域滤波操作的结果中，得到滤波后的特定波动成分的频率
‑ꢀ
波数域数据；
[0031]
s900.对s800中得到的所述滤波后的频率-波数域数据进行频率-波数域反变换操作，得到只包含特定波动成分的时间-空间域海底地震波场数据；所述只包含特定波动成分的时间-空间域海底地震波场数据即为本分离方法的最终结果。
[0032]
优选地，所述时间-空间域数据由人工预置的地震检波器阵列上的阵元接收得到，按下式表达：
[0033]
f(t,x)＝u(t,x) n(t,x)
[0034]
其中：u(t,x)为所述有用信号；n(t,x)为所述干扰噪声；t表示时间；x表示空间。
[0035]
优选地，s210中采用二维傅里叶变换对所述时间-空间域数据进行频率-波数变换操作，将信号从时间-空间域信号转换到频率-波数域中；所述频率-波数变换操作按下式表达：
[0036][0037]
其中：f为所述时间-空间域数据的信号频率；e为自然对数；j为计数器，且j＝0,1,2...n；n为空域傅里叶变换的采样点数，由人工预设；k为所述时间-空间域数据的波数，按下式表达：
[0038]
k＝1/λ
[0039]
其中：λ为所述时间-空间域数据的波长；
[0040]
所述频率-波数域数据按下式表达：
[0041]
f(f,k)＝u(f,k) n(f,k)
[0042]
其中：f(f,k)为时间-空间域数据f(t,x)的频率-波数谱；u(f,k)为有用信号u(t,x)的频率-波数谱；n(f,k)为干扰噪声n(t,x)的频率-波数谱。
[0043]
优选地，s230中将所述海底地震波场的所述频率-波数域数据进行信号离散化处理，得到频率-波数域图像，按下式表达：
[0044]
f(iδf,jδk)＝u(iδf,jδk) n(iδf,jδk)
[0045]
其中：f(iδf,jδk)为离散化的时间-空间域数据的频率-波数谱；u(iδf,jδk) 为离散化的有用信号的频率-波数谱；n(iδf,jδk)为离散化的干扰噪声的频率
‑ꢀ
波数谱；，δf为频域分辨率，按下式表达：
[0046]
δf＝1/(mδt)
[0047]
其中：m为时域傅里叶变换变换的采样点数，由人工预设；δt为时间采样间隔，由人工预设，且满足按下式表达的采样定律：
[0048][0049]
其中：f
max
为时间-空间域数据的频率的最大值；
[0050]
δk为波数域分辨率，按下式表达：
[0051]
δk＝1/(nδx)
[0052]
其中：n为空域傅里叶变换的采样点数，由人工预设；δx为空间采样间隔，由人工预设，且满足按下式表达的采样定律：
[0053][0054]
其中：k
max
为时间-空间域数据中波数的最大值；
[0055]
i为计数器，且i＝0,1,2...m。
[0056]
优选地，s400中所述频率-波数域图像能量均值按下式表达：
[0057][0058]
其中：e
′
为所述频率-波数域图像能量均值；s
′
为所述频率-波数域图像的像素矩
阵，按下式表达：
[0059][0060]
其中：s为所述频率-波数域图像所对应的应力矩阵；
[0061]s′
(i,j)为s
′
在(i,j)处的像素值；
[0062]
所述图像能量标准差按下式表达：
[0063][0064]
其中：δ
′
为所述图像能量标准差。
[0065]
优选地，所述图像能量阈值按下式表达：
[0066]vth
′
＝e
′
α
′
δ
′
[0067]
其中：v
th
′
为所述图像能量阈值；α
′
为所述标准差倍数，由人工预设。
[0068]
优选地，s700中采用八连通区域法或四连通区域法进行所述连通性分析操作。
[0069]
优选地，s800中的所述频率-波数域滤波操作中采用二维滤波器，按下式表达：
[0070][0071]
其中：为所述滤波后的频率-波数域数据；y(f,k)为海底地震波数据进行二维傅里叶变换得到的矩阵；h(f,k)为滤波函数，按下式表达：
[0072][0073]
其中：b为波动成分对应的连通区域；
[0074]
优选地，s900中所述频率-波数域反变换操作采用二维逆傅里叶变换，按下式表达：
[0075][0076]
本发明与现有技术对比，具有以下优点：
[0077]
1.由于本发明利用了连续海底地震波场中波动成分的能量特性，通过设置阈值和选择连通域，将连续海底地震波场的主要波动成分自动分离开来，从而可以自动精准分离水中目标诱发的多种同频同向的连续海底地震波波动成分；
[0078]
2.由于本发明采用了阈值分割和连通域选择的方式，从而极大的提高了波动成分的分辨率，尤其适合于多种波动成分的波速相近的情形。
附图说明
[0079]
图1为本发明具体实施例的流程示意图；
[0080]
图2a为本发明具体实施例的八连通区域示意图；
[0081]
图2b为本发明具体实施例的四连通区域示意图；
[0082]
图3为本发明具体实施例的频率-波数滤波原理示意图；
[0083]
图4为本发明具体实施例的点声源诱发的海底地震波场模型；
[0084]
图5为本发明具体实施例的海底为砾岩时连续点声源诱发的海底地震波场数据的频率-波数域图像；
[0085]
图6为本发明具体实施例的海底为砾岩时海底地震波场中波动成分组成示意图；
[0086]
图7为本发明具体实施例的海底为砾岩时海底地震波场中主要波动成分的波速示意图；
[0087]
图8为本发明具体实施例的海底为砾岩时连续点声源诱发的海底地震波场数据的频率-波数域图像；
[0088]
图9为本发明具体实施例的海底为砾岩时海水-海底分界面处30hz时法向应力的切片示意图；
[0089]
图10a为本发明具体实施例的当海底为砾岩时，α
′
值为2时第一次阈值分割及连通域提取后的图像；
[0090]
图10b为本发明具体实施例的当海底为砾岩时，α
′
值为3时第一次阈值分割及连通域提取后的图像；
[0091]
图10c为本发明具体实施例的当海底为砾岩时，α
′
值为4时第一次阈值分割及连通域提取后的图像；
[0092]
图10d为本发明具体实施例的当海底为砾岩时，α
′
值为5时第一次阈值分割及连通域提取后的图像；
[0093]
图11a为本发明具体实施例的在图7的30hz处，当α
′
值为2时的切片示意图；
[0094]
图11b为本发明具体实施例的在图7的30hz处，当α
′
值为3时的切片示意图；
[0095]
图11c为本发明具体实施例的在图7的30hz处，当α
′
值为4时的切片示意图；
[0096]
图11d为本发明具体实施例的在图7的30hz处，当α
′
值为5时的切片示意图；
[0097]
图12a为本发明具体实施例的将图11a的结果与图9的结果，统一为像素值后的对比示意图；
[0098]
图12b为本发明具体实施例的将图11b的结果与图9的结果，统一为像素值后的对比示意图；
[0099]
图12c为本发明具体实施例的将图11c的结果与图9的结果，统一为像素值后的对比示意图；
[0100]
图12d为本发明具体实施例的将图11d的结果与图9的结果，统一为像素值后的对比示意图；
[0101]
图13a为本发明具体实施例的海底为砾岩时滤波的海底地震波场时间-空间域图像的原有地震波信号部分截图；
[0102]
图13b为本发明具体实施例的海底为砾岩时滤波的海底地震波场时间-空间域图像的分离出的scholte波信号部分截图。
[0103]
其中：1.无反射型的声压阻抗边界；2.弹性单元的无限元边界。
具体实施方式
[0104]
下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本技术所附权利要求所限定的范围。
[0105]
一种同频同向的连续海底地震波场波动成分分离方法，包含以下步骤：
[0106]
s100.获取海底地震波场的时间-空间域数据；时间-空间域数据为关于时间和空间的二维函数。
[0107]
本具体实施例中，时间-空间域数据由人工预置的地震检波器阵列上的阵元接收得到，按式(1)表达：
[0108]
f(t,x)＝u(t,x) n(t,x)
ꢀꢀꢀꢀꢀꢀ
(1)
[0109]
其中：u(t,x)为有用信号；n(t,x)为干扰噪声；t表示时间；x表示空间。
[0110]
需要说明的是，实际应用中，除本具体实施例中采用的实际采集的方式获得时间-空间域数据，还可以通过有限元软件comsol multiphysics仿真海底地震波场，也可以通过开展海上试验测量得到。
[0111]
s200.将时间-空间域数据转化为频率-波数域图像，具体包含以下步骤：
[0112]
s210.对时间-空间域数据进行频率-波数变换操作，获得海底地震波场的频率-波数域数据；频率-波数域数据包含有用且需提取的有用信号和无用且需去除的干扰噪声。
[0113]
本具体实施例中，在s210中采用二维傅里叶变换对时间-空间域数据进行频率-波数变换操作，将信号从时间-空间域信号转换到频率-波数域中；频率-波数变换操作按式(2)表达：
[0114][0115]
其中：f为时间-空间域数据的信号频率；e为自然对数；j为计数器，且 j＝0,1,2...n；n为空域傅里叶变换的采样点数，由人工预设；k为时间-空间域数据的波数，按式(3)表达：
[0116]
k＝1/λ
ꢀꢀꢀꢀ
(3)
[0117]
其中：λ为时间-空间域数据的波长。
[0118]
频率-波数域数据按式(4)表达：
[0119]
f(f,k)＝u(f,k) n(f,k)
ꢀꢀꢀꢀ
(4)
[0120]
其中：f(f,k)为时间-空间域数据f(t,x)的频率-波数谱；u(f,k)为有用信号u(t,x)的频率-波数谱；n(f,k)为干扰噪声n(t,x)的频率-波数谱。
[0121]
s220.然后记录下海底地震波场的频率-波数域数据的实部和虚部。
[0122]
需要说明的是，在s220中，还需要根据频率-波数域数据的波动成分的速度在频率-波数域数据中指示出对应波动成分在频率-波数域数据的大致分布范围，为后续分析缩小搜索范围。这里的波动成分，就是指的有用信号和干扰噪声。
[0123]
需要进一步说明的是，如果时间-空间域数据来自于仿真数据对应的波动成分速度则通过理论公式进行求解；而如果时间-空间域数据通过海上试验测量得到，则试验中的
波动成分速度则是通过在该海域的气枪声源试验获得，即水中气枪为声源，以布放于海底的地震波阵列接收由气枪声源激发的海底地震波信号，采用tao-p法获得scholte波和水声信号等波动成分的传播速度。
[0124]
s230.然后将海底地震波场的频率-波数域数据进行信号离散化处理，得到频率-波数域图像。
[0125]
本具体实施例中，在s230中将海底地震波场的频率-波数域数据进行信号离散化处理，得到频率-波数域图像，按式(5)表达：
[0126]
f(iδf,jδk)＝u(iδf,jδk) n(iδf,jδk)
ꢀꢀꢀꢀ
(5)
[0127]
其中：f(iδf,jδk)为离散化的时间-空间域数据的频率-波数谱；u(iδf,jδk) 为离散化的有用信号的频率-波数谱；n(iδf,jδk)为离散化的干扰噪声的频率
‑ꢀ
波数谱；，δf为频域分辨率，按式(6)表达：
[0128]
δf＝1/(mδt)
ꢀꢀꢀꢀꢀꢀꢀꢀ
(6)
[0129]
其中：m为时域傅里叶变换变换的采样点数，由人工预设；δt为时间采样间隔，由人工预设，且满足按式(7)表达的采样定律：
[0130][0131]
其中：f
max
为时间-空间域数据的频率的最大值。
[0132]
δk为波数域分辨率，按式(8)表达：
[0133]
δk＝1/(nδx)
ꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0134]
其中：n为空域傅里叶变换的采样点数，由人工预设；δx为空间采样间隔，由人工预设，且满足按式(9)表达的采样定律：
[0135][0136]
其中：k
max
为时间-空间域数据中波数的最大值。
[0137]
i为计数器，且i＝0,1,2...m。
[0138]
s300.将频率-波数域图像进行像素化处理，得到像素化后的频率-波数域图像；频率-波数域图像的行为频率，列为波数；频率-波数域图像的每个像素点的像素值的对应范围为[0-255]。
[0139]
需要说明的是，如果将海底地震波场中待提取的波动成分看作是有用信号，其余波动成分看作是干扰噪声，对经过频率-波数域正变换中的海底地震波场的频率-波数域图像进行阈值分割处理，从而就可以实现波动成分的提取；其基本思路是：将频率-波数域图像像素化处理，将原图像中的点转换为像元，将对应的像元赋予对应的像素值，对波动成分包括的频率范围内逐一扫频做分割处理，从而将波动成分的所有能量提取出来。
[0140]
s400.按频率-波数域图像的行进行阈值划分操作，计算得到频率-波数域图像的每一行的频率-波数域图像能量均值和图像能量标准差。
[0141]
本具体实施例中，在s400中频率-波数域图像能量均值按式(10)表达：
[0142][0143]
其中：e
′
为频率-波数域图像能量均值；s
′
为频率-波数域图像的像素矩阵，按式(11)表达：
[0144][0145]
其中：s为频率-波数域图像所对应的法向应力矩阵。
[0146]s′
(i,j)为s
′
在(i,j)处的像素值。
[0147]
图像能量标准差按式(12)表达：
[0148][0149]
其中：δ
′
为图像能量标准差。
[0150]
需要说明的是，通过式(11)可将频率-波数域图像中的能量对应转换成0-255 范围内的像素值。
[0151]
需要进一步说明的是，将海底地震波场中scholte波视为有用信号，其余信号，如第一阶、第二阶简正波、纵波等，则为干扰噪声；或者也可以降简正波视为有用信号，scholte波等波动成分为干扰噪声；从本质上来说没有区别，因为二者是相对的；甚至还可以是其他波动成分的组合当成有用信号，其他的成分作为干扰噪声。为根据频率-波数图中各波动成分整体能量起伏地计算分割阈值，以分离频率-波数域图像中的主要波动成分。
[0152]
需要进一步说明的是，s400的阈值划分方法是通过以下过程得到的：
[0153]
首先，为了实现对频率-波数域图像的阈值划分，其阈值划分公式按式(13)、式(14)和式(15)表达：
[0154]vth
＝e αδ
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(13)
[0155][0156][0157]
其中：v
th
为选取的阈值；e为频率-波数域图像整体能量的平均值，由对整幅图像所有像元值s(i,j)求加权平均值得到；δ为图像能量的方差；α为一个加权系数，其选取原则遵从正态分布中的3δ原则，即假设频率-波数域图像中各点的能量服从高斯分布。
[0158]
具体来说：当波场中大于某阈值的能量占整个波场能量的0.1％时，将α设置为3，若使波场中大于该阈值的能量占整个波场能量的0.001％时，将α设置为5，即该值的选取根据要提取的波动成分的能量与波场整体能量的占比来确定。
[0159]
于是，根据上述说明，就可以按照下述方法进行阈值划分：
[0160]
在获得海底地震波场的频率-波数域图像后，将频率-波数域图像中的频率看做行，波数看做列，按照图像的行进行阈值分割，计算每一行图像能量的均值和标准差δ
′
，将这一行的阈值v
th
′
设置为该行均值与α
′
倍的该行标准差之和，若该行中像素点处的能量值
小于该阈值v
th
′
，则将该像素点对应的能量值则设为0。相应的计算表达式即为式(13)、式(14)和式(15)。
[0161]
以上就是本发明的阈值划分的推理过程。
[0162]
s500.根据频率-波数域图像能量均值和图像能量标准差，将频率-波数域图像的每一行的图像能量阈值设为该行的图像能量均值与人工预设的标准差倍数乘以该行的图像能量标准差之和。
[0163]
本具体实施例中，图像能量阈值按式(16)表达：
[0164]vth
′
＝e
′
α
′
δ
′ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(16)
[0165]
其中：v
th
′
为图像能量阈值；α
′
为标准差倍数，由人工预设。
[0166]
s600.对频率-波数域图像进行阈值分割操作，具体包含以下步骤：
[0167]
s610.逐一将每一行中的每个像素点，与该像素点所在的行的图像能量阈值进行大小对比。
[0168]
s620.根据s610的对比结果，做出如下操作：
[0169]
如果该像素点的像素值比该像素点所在的行的图像能量阈值小，则将该像素点的像素值置0。
[0170]
否则，保持该像素点的像素值的值不变。
[0171]
s630.再次执行s610，直至比对完频率-波数域图像中所有的像素点。
[0172]
需要说明的是，s600到s630就是所谓的阈值分割操作。
[0173]
s700.对经过阈值分割操作后的频率-波数域图像进行连通性分析操作，搜索得到有用信号的所在的连通域和干扰噪声的所在的连通域；
[0174]
本具体实施例中，在s700中采用八连通区域法或四连通区域法进行连通性分析操作。
[0175]
需要说明的是，通常，海底地震波中波动成分的能量以一个中心区域向外扩散，在频率-波数域图像上形成一定的展宽，为了将这个区域中的能量从频率
‑ꢀ
波数域图像上提取出来，可以采用以阈值分割为基础的连通区域来实现，即认为远离中心区域、能量过低的地方，不是需要提取的波动成分的重要能量分布区域。
[0176]
在划分连通域的过程中，通常采用八连通区域法与四连通区域法。
[0177]
如图2a所示，八连通区域法提取按照遍历的方式逐个能量点进行阈值判断，对高于阈值的能量点，则同时判断该点上、下、左、右、左上、右上、左下、右下等八个方向是否存在超过阈值的能量点。如果存在，则认为该点属于要提取的波动成分的能量分布区域。
[0178]
同理，如图2b所示，四连通区域是在该点的上、下、左、右四个方向进行判断。
[0179]
需要进一步说明的是，八连通区域适合于能量向多方向分散型的波动成分提取，由于需要向八个方向搜索，故计算量较大；四连通区域则适合于在水平向和垂直向分散的波动成分提取，其计算量相对前者明显减小。
[0180]
需要进一步说明的是，实际研究表明，海底地震波信号分布在一定的甚低频段内，表现为一簇簇的线谱，可将这些频带切割为多个单频点，提取每个单频点对应的波动成分能量，这是本方法解决同频同向信号分离的重要体现，最终将分离方法处理后的每个频点能量合成为频带内的波动成分能量，因此对于每个单频点的处理，选用四连通区域即可，故计算量也随之减少。
[0181]
由上可知，在实际应用过程中，需要根据使用环境以及用户需求来灵活调整，而不是一成不变的固定模式。
[0182]
s800.对s700中搜索得到特定波动成分，例如水声信号或scholte波，所在的连通域进行频率-波数域滤波操作；然后将s220中的海底地震波场的频率-波数域数据的虚部补充进入频率-波数域滤波操作的结果中，得到滤波后的特定波动成分的频率-波数域数据。
[0183]
本具体实施例中，在s800中的频率-波数域滤波操作中采用二维滤波器，按式(17)表达：
[0184][0185]
其中：为滤波后的频率-波数域数据，亦即滤波后的结果；y(f,k)为海底地震波数据进行二维傅里叶变换得到的矩阵；h(f,k)为滤波函数，通常采用区域滤波的形式，按式(18)表达：
[0186][0187]
其中：b为波动成分对应的连通区域。
[0188]
需要说明的是，按波数的定义，并由式(3)可推得将时间-空间域的信号变换到频率-波数域的公式，按式(19)表达：
[0189][0190]
其中：v为波动成分的波速。
[0191]
由式(19)可知波动成分波数、频率和传播速度三者之间的关系：在如图3 所示的频率-波数平面上，某种波动成分对应的斜率越大，其传播速度越高，如果波动成分的频率分布在f1和f2之间、速度分布在v1和v2之间时，假设区域a 为高速波动成分干扰区，区域c为低速波动成分干扰区，区域b为有效信号区。
[0192]
s900.对s800中得到的滤波后的频率-波数域数据进行频率-波数域反变换操作，得到只包含特定波动成分的时间-空间域海底地震波场数据；只包含特定波动成分的时间-空间域海底地震波场数据即为本分离方法的最终结果。
[0193]
本具体实施例中，在s900中频率-波数域反变换操作采用二维逆傅里叶变换，按式(20)表达：
[0194][0195]
需要说明的是，式(20)是通过二维逆傅里叶变换，将海底地震波场的频率
ꢀ‑
波数数据变换到时间-空间域，得到滤波后的波动成分的时间域-空间域数据。
[0196]
需要进一步说明的是，在对式(18)滤波的结果进行傅里叶反变换时，需要进行补充虚部的操作，否则会导致频域-波数域数据的背景丢失。
[0197]
为了进一步展示本发明方法的有效性和技术优势，本具体实施了还包含了以下的
算例，进行进一步的说明：
[0198]
需要事先说明的是，以浅海中连续点声源诱发的海底地震波场为算例，采用有限元软件comsol multiphysics建立海水、点声源、海底三者的耦合动力学模型。其中，其中，模型示意图如图4所示，海水、海底对应的介质参数见表1：
[0199]
表1介质类型及参数
[0200]
介质类型纵波速度(m/s)横波速度(m/s)密度(g/cm3)海水150001.0砾岩20008002.1
[0201]
如图4所示，无反射型的声压阻抗边界设置在海水层两侧以实现声压的吸收，从而模拟水平方向上的无限海域空间。考虑海水-空气界面为压力释放型边界，另外点声源设置在网格节点处，激励类型设置为声学压强。
[0202]
将图4中点声源的深度设置为90m，加载30hz的正弦载波，根据式(5)，将时间采样率设置为1000hz，采样时长为4s，空间采样率设置为10m，采样距离为2000m，等间隔采集连续点声源诱发的海底地震波场数据(海水-海底分界面处的法向应力)，将其汇总为第1节中时间-空间的形式。对采集得到的海底地震波场的时间-空间域的应力做频率-波数变换，得到的结果如图5所示。
[0203]
如图6所示，为当海底为砾岩时海底地震波场中的波动成分。
[0204]
其中，海底地震波场中的横波、侧面波以及泄露瑞利波由于衰减较快，其能量较小，考虑海底地震波场的主要波动成分时通常将其忽略。
[0205]
海底为砾岩时，海底地震波场中schole波与简正波波速的计算公式按式(21) 和式(22)表达：
[0206][0207][0208]
其中：c为海底地震波场中波动成分的波速；cf为海水介质中声波的波速；ρf为海水介质的密度；c
p2
为海底介质中纵波的波速；c
s2
为海底介质中横波的波速。
[0209]
如图7所示，通过求解式(21)可得到砾岩海底，30hz时各波动成分的波速。
[0210]
根据理论对图6中主要波动成分的速度进行计算，其对应的频率与波数值见表2：
[0211]
表2仿真条件下海底地震波场中主要波动成分对应的波速和波数
[0212]
对比项传播速度(m/s)波数(m-1
)scholte波6990.0429第一阶简正波15600.0192第二阶简正波17130.0175纵波20000.0150
[0213]
如图8所示，为根据表2中的数据将对应的波动成分在图5中标记得到的结果。
[0214]
从图8中可清晰看出，海底为砾岩时点声源激发的海底地震波场中包含纵波、第一、二阶简正波以及scholte波等波动成分，且各波动成分在频率和波数上存在一定展宽(即其能量不能简单的视为一个“点”)。
[0215]
如图9所示，为了分析其相对能量的大小，选取频率-波数域图像中30hz的切片。
[0216]
由图9可知，各波动成分对应的波数与表2中的理论值一致，侧面验证了仿真的正确性。从波动成分的能量角度分析，海底为砾岩时海底地震波场中第一阶简正波能量最大，scholte波次之，第二阶简正波与纵波能量最小，且在各波动成分中间掺杂着海洋环境噪声，导致各波动成分的分辨率较低。为此，采用本文提出的基于阈值滤波的频率-波数分离方法进行改进。
[0217]
如图10a～图10d所示，为当海底为砾岩时，不同α
′
值分别为2、3、4、5时第一次阈值分割及连通域提取后的图像。
[0218]
由图10a～图10d可知，α
′
分别取值为2、3、4、5时，硬海底条件下，四连通区域的数量对应为3、2、1、1；海底为砾岩条件下，α
′
对应为5、3、2、2；两种海底条件下，随着α
′
值的增大，连通区域数量均逐渐减小，最终趋于稳定。
[0219]
本具体实施例中，提供了如图11a～图11d所示的在图7的30hz处，当α
′
分别取值为2、3、4、5时的切片图。
[0220]
从图11a～图11d可知，海底为砾岩时，当α
′
＝2时，海底地震波场中的背景噪声由于能量较小，自动被滤除，只剩下纵波、第二阶简正波，第一阶简正波与scholte波；随着α
′
值增大，第二阶简正波由于能量过小被自动滤除，在α
′
＝4、α
′
＝5时波场中只存在第一阶简正波与scholte波，即二者被自动分离出来。
[0221]
为了更好地对比两种方法的结果，基于式(7)和式(9)，将图11中的结果与图9中的结果均统一为像素值，得到的结果对比图如图12a～图12d所示。
[0222]
由图12可知，基于阈值分割的频率-波数分离方法的分离精度与分辨率较高，还可将波场中的背景噪声及杂波滤除。
[0223]
海底地震波场中，水声信号(通常指第一阶简正波)与scholte波包含海底底质与目标位置的相关信息，通常将二者视为式(4)中的有用信号u(t,x)，其余信号则是被视为干扰n(t,x)，根据s700和s800得到滤波前以及滤波后的归一化海底地震波场的时间-空间域数据。以scholte波的分离为例，作出滤波前与滤波后海底地震波场的时间-空间域数据如图13a和图13b所示，其余波动成分的分离与之类似。
[0224]
从图13a和图13b中可以看出，海底为砾岩时，本方法将海底地震波场中的 scholte波与其他波动成分分离，同样用该方法还可对水声信号进行有效分离。
[0225]
至此，本具体实施例用一个具体的算例，证明了本发明的有效性以及相对现有技术的优势。
[0226]
在上述的详细描述中，各种特征一起组合在单个的实施方案中，以简化本公开。不应该将这种公开方法解释为反映了这样的意图，即，所要求保护的主题的实施方案需要比清楚地在每个权利要求中所陈述的特征更多的特征。相反，如所附的权利要求书所反映的那样，本发明处于比所公开的单个实施方案的全部特征少的状态。因此，所附的权利要求书特此清楚地被并入详细描述中，其中每项权利要求独自作为本发明单独的优选实施方案。
[0227]
为使本领域内的任何技术人员能够实现或者使用本发明，上面对所公开实施例进行了描述。对于本领域技术人员来说；这些实施例的各种修改方式都是显而易见的，并且本文定义的一般原理也可以在不脱离本公开的精神和保护范围的基础上适用于其它实施例。因此，本公开并不限于本文给出的实施例，而是与本技术公开的原理和新颖性特征的最广范围相一致。
[0228]
上文的描述包括一个或多个实施例的举例。当然，为了描述上述实施例而描述部件或方法的所有可能的结合是不可能的，但是本领域普通技术人员应该认识到，各个实施例可以做进一步的组合和排列。因此，本文中描述的实施例旨在涵盖落入所附权利要求书的保护范围内的所有这样的改变、修改和变型。此外，就说明书或权利要求书中使用的术语“包含”，该词的涵盖方式类似于术语“包括”，就如同“包括，”在权利要求中用作衔接词所解释的那样。此外，使用在权利要求书的说明书中的任何一个术语“或”是要表示“非排它性的或者”。
[0229]
以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
[0230]
最后，应当指出，以上实施例仅是本发明较有代表性的例子。显然，本发明不限于上述实施例，还可以有许多变形。凡是依据本发明的技术实质对以上实施例所做的任何简单修改、等同变化及修饰，均应认为属于本发明的保护范围。
[0231]
在此，需要说明的是，上述技术方案的描述是示例性的，本说明书可以以不同形式来体现，并且不应被解释为限于本文阐述的技术方案。相反，提供这些说明将使得本发明公开将是彻底和完整的，并且将向本领域技术人员充分传达本说明书所公开的范围。此外，本发明的技术方案仅由权利要求的范围限定。
[0232]
用于描述本说明书和权利要求的各方面公开的形状、尺寸、比率、角度和数字仅仅是示例，因此，本说明书和权利要求的不限于所示出的细节。在以下描述中，当相关的已知功能或配置的详细描述被确定为不必要地模糊本说明书和权利要求的重点时，将省略详细描述。
[0233]
在使用本说明书中描述的“包括”、“具有”和“包含”的情况下，除非使用否则还可以具有另一部分或其他部分，所用的术语通常可以是单数但也可以表示复数形式。
[0234]
应该指出，尽管在本说明书可能出现并使用术语“第一”、“第二”、“顶部”、“底部”、“一侧”、“另一侧”、“一端”、“另一端”等来描述各种不同的组件，但是这些成分和部分不应受这些术语的限制。这些术语仅用于区分一个成分和部分和另一个成分和部分。例如，在不脱离本说明书的范围的情况下，第一部件可以被称为第二部件，并且类似地，第二部件可以被称为第一部件，顶部和底部的部件在一定情况下，也可以彼此对调或转换；一端和另一端的部件可以彼此性能相同或者不同。
[0235]
此外，在构成部件时，尽管没有其明确的描述，但可以理解必然包括一定的误差区域。
[0236]
在描述位置关系时，例如，当位置顺序被描述为“在...上”、“在...上方”、“在... 下方”和“下一个”时，除非使用“恰好”或“直接”这样的词汇或术语，此外则可以包括它们之
间不接触或者接触的情形。如果提到第一元件位于第二元件“上”，则并不意味着在图中第一元件必须位于第二元件的上方。所述部件的上部和下部会根据观察的角度和定向的改变而改变。因此，在附图中或在实际构造中，如果涉及了第一元件位于第二元件“上”的情况可以包括第一元件位于第二元件“下方”的情况以及第一元件位于第二元件“上方”的情况。在描述时间关系时，除非使用“恰好”或“直接”，否则在描述“之后”、“后续”、“随后”和“之前”时，可以包括步骤之间并不连续的情况。本发明的各种实施方案的特征可以部分地或全部地彼此组合或者拼接，并且可以如本领域技术人员可以充分理解的以各种不同地构造来执行。本发明的实施方案可以彼此独立地执行，或者可以以相互依赖的关系一起执行。