一种海底硫化物矿体结构宽频带地震约束反演成像方法与流程

1.本发明涉及海洋地震探测领域，具体是一种海底硫化物矿体结构宽频带地震约束反演成像方法。


背景技术：

2.海底热液循环系统伴生富含铜、锌和金的热液硫化物这一重要海底资源，是当前国际海洋科学领域的前沿。海底硫化物矿体结构的精准成像是亟待解决的重要问题，是进行资源评价的最为重要的基础。
3.当前国际上直接研究矿体结构很大程度上只能依赖于大深度钻探，在钻探数量和分布上具有很高的要求，成本和作业方式受到很大限制。由于海底热液硫化物矿体具有分布面积小、水深大和结构复杂的特点，其地球物理响应剖面仅有几百米甚至几十米，基于深海拖曳平台或各类潜水器(auv、rov等)的近底地球物理方法在硫化物矿体勘查中作用显著。然而，在硫化物矿体勘查中受方法本身分辨率的限制，近底磁法只能大致刻画蚀变带范围，近底电法只能大致刻画矿体横向范围(zhu et al.,jgr,2020，zl201910589392.2，zl202010801649.9)，但无法准确获得矿体深度信息，矿体结构也无法准确成像。
4.近底地震是当近年来用于硫化物矿体结构研究的新颖方法(zl201821239127.9)，通过将震源和/或接收器阵列拖曳于近海底海水中或垂直固定于海底进行探测，配合中高频震源还可以有效提高探测分辨率。已有硫化物地震探测结果表明，垂直地震可能是获取硫化物矿体内部结构的一种有效方式，但基于等效偏移距(eom)的垂直地震速度建模方法(asakawa and mcintyre,seg annual meeting,2015)是一种基于kirchhoff时间偏移的方法，速度建模精度不足造成矿体成像结果中同相轴没有正确归位(产生“横条”状假象)，制约了矿体结构解译。
5.从理论上讲，地震全波形反演(full wavefield inversion，fwi)是解决速度建模问题的理想方案，对于海底硫化物近底地震成像具有应用前景。但受海底硫化物分布区域较小、矿体结构复杂、物性特征不清晰等方面的限制，直接使用fwi进行硫化物矿体结构的反演成像不太现实，同时含有高频信息和低频信息的宽频带地震才能适用于硫化物矿体结构成像。因此，发展一种海底硫化物矿体结构的宽频带地震反演成像方法具有重要的实用价值，将有力推动我国海底硫化物资源评价重大任务的执行，同时也对海底探测能力的提高具有重要的意义。


技术实现要素：

6.针对上述问题中存在的不足之处，本发明提供一种海底硫化物矿体结构宽频带地震约束反演成像方法。
7.本发明的技术方案如下：
8.本发明提供了一种海底硫化物矿体结构宽频带地震约束反演成像方法，其包括如下步骤：
9.步骤1：对硫化物矿体样品进行物性测试，获取硫化物矿体有限范围的声速、密度、孔隙度；
10.步骤2：对垂直电缆地震数据进行倾斜状态估计、定位信息融合、直达波压制、数据去噪、振幅恢复、数据反褶积处理；
11.步骤3：对垂直电缆地震数据进行压缩感知高精度重建处理；
12.步骤4：多源数据约束初始速度建模；
13.依据步骤1物性测试获取的围岩、底质和岩心的声速特征，构建一个平滑速度模型va，平滑速度模型包括海底以上的水体部分和海底以下的固体部分；
14.步骤5：多步频率阻尼多尺度全波形反演；
15.全波形反演是通过优化方法寻找目标函数的极小值，满足目标函数最小值的速度即为最终反演结果；利用近底地震多方向和宽频带的特点，分步使用频率阻尼波场全波形反演方法从低频到高频进行频率域全波形反演，实现由轮廓到细节的逐步重建，建立准确速度模型vb。
16.作为本发明的优选方案，所述的步骤1为：
17.以已有资料大致确定的硫化物矿体的展布范围为参考，选择在硫化物矿体上不同位置和不同深度的底质取样和钻探岩心样品，开展常压和高压环境的物性测试，获取硫化物矿体有限范围的包括声速、密度、孔隙度在内的物性特征。
18.作为本发明的优选方案，所述的步骤2为：
19.利用超短基线数据确定测线作业的震源和垂直电缆的水下位置，并基于姿态仪数据进行垂直电缆倾斜状态估计，获取各接收点的位置；对近底地震数据进行定位信息融合、直达波压制、数据去噪、振幅恢复、数据反褶积处理。
20.作为本发明的优选方案，所述的步骤3中：
21.使用基于迭代阈值的压缩感知稀疏约束重建方法对垂直电缆地震数据进行压缩重建处理。由于震源和接收两个方向上均可能有数据缺失，需要在两个方面进行数据的压缩感知重建，对缺失炮和道的数据进行恢复：
22.a)地震采集的数据是成百上千炮数据，对于某些缺失或无效的炮数据，首先进行共接收点道集数据的重建，获取规则共接收点道集数据；
23.b)将上一步重建的数据转换为共炮道集数据，并对每个共炮道集数据进行重建，获取规则的共炮道集数据。
24.作为本发明的优选方案，所述的步骤3中，基于迭代阈值的压缩感知重建方法具体为：
25.根据压缩感知理论，数据缺失的过程表示为：
26.b＝rf
27.其中，f为规则网格完整的地震数据，设共有m道数据，b为含缺失道的观测地震数据，设共有n道数据，缺失了m-n道数据，r为缺失采样矩阵，设矩阵维度为n
×
m；对于地震数据而言，其波形使用curvelet变换或者seislet变换进行稀疏表示，设变换矩阵为s，则：
28.b＝ax，with a＝rs
*
29.其中，x为数据f在s域的稀疏表示，维度为m，a＝rs
*
为测量矩阵，上标
*
表示复数的共轭；这样数据的重建问题表示为：
[0030][0031]
其中，变量上面的波浪线～表示对变量的估计，λ为拉格朗日乘子；
[0032]
由于数据缺失具有随机性，r为随机采样矩阵，在满足稀疏度要求(kdlog2(m/n)≤n，其中，k为稀疏度，d为一个常数)的情况下，使用稀疏约束反演方法求解上式，获取最终数据重建结果
[0033]
使用迭代阈值方法进行稀疏约束反演，迭代式如下：
[0034][0035]
其中，i为迭代次数，表示阈值运算，γi为迭代步长，θ为使用的阈值。
[0036]
作为本发明的优选方案，所述的步骤4中，模型的固体部分和水体部分建模方法如下：
[0037]
1)平滑速度模型的固体部分，其背景速度为围岩速度，矿体范围以已有其他资料(如电法探测)勾勒的硫化物矿体轮廓为基础，块状硫化物使用岩心速度，网脉状硫化物的速度根据已有慢性扩张洋中脊硫化物区两层之间的关系和目标区实测块状硫化物的速度确定；
[0038]
2)平滑速度模型中水体部分，通过在探测区域开展若干站位的ctd探测作业获取的ctd数据建立水体声速剖面。使用的方法为：以若干站位的ctd获取的一定深度的声速数据作为控制点，通过数据拟合和插值建立整个探测剖面的声速数据。
[0039]
作为本发明的优选方案，所述的步骤5中，多步频率阻尼多尺度全波形反演的实现方法如下：
[0040]
a)基于对地震探测使用震源的一般性认识，通过对地震数据进行频谱分析，确定地震数据频率成分所集中的范围，根据地震数据的频率分布带数量确定需要执行反演的步数；
[0041]
b)根据地震数据的频率特征，确定每个频率分布带的起始频率fb和终止频率fe，要保证fb～fe之间涵盖该频率分布带的所有有效频率成分；
[0042]
c)对地震数据进行频率域波场阻尼处理，获取每个频率分布带的阻尼波场u
α
(f)，计算方法如下：
[0043][0044]
其中，u(f)为频率域波场，γd(f)为频率阻尼算子，f为频率，α为和反演尺度相关的阻尼因子，用于控制反演尺度，0＜α＜1表示波场按照指数曲线压制高频，α＝0表示波场
未阻尼；
[0045]
d)分步进行全波形反演，每步对应1个地震数据的频率分布带，且前者反演的结果作为下一步反演的初始速度模型。
[0046]
作为本发明的优选方案，所述的步骤5中，多步频率阻尼多尺度全波形反演中，在进行波形反演的过程中增加波场照明，以实现对深部信息的高分辨成像；每次迭代计算一次总的波场照明矩阵，计算公式如下：
[0047][0048]
其中，imi为第i次迭代的波场照明矩阵，vi为第i次迭代生成的速度模型，为模型vi下的第j个震源在位置r处的正演波场，为模型vi下第k个检波点在位置r处的逆传播波场，ns为震源总数，ng为接收点总数。全波形反演中使用以下带有照明信息的模型修改方法：
[0049][0050]
其中，δvi为通过第i次迭代反演出来的模型修改量，模型v
i+1
为第i+1次迭代生成的速度模型。当迭代次数达到预定值或满足预设迭代终止条件时停止迭代，获取最终速度模型vb。预设终止条件为：
[0051][0052]
其中，e(i)为经过i次迭代后目标函数的最小二乘误差，η为一个不小于零的正数。
[0053]
与现有技术相比，本发明所具有的有益效果至少包括：
[0054]
(1)本发明在反演中充分利用了包含低频和高频信息的混合震源数据，这种宽频带数据十分有利于提高反演的可靠性和精准度；
[0055]
(2)本发明通过岩石样品物性测量获取海底物性参数信息，并充分利用已有电法探测、岩心资料、若干站位的ctd数据等多源数据构建用于全波形反演的更为可靠的初始速度模型，大大提高了速度反演的有效性和准确性；
[0056]
(3)根据对数据频率成分和结构的分析，构建了多步频率阻尼多尺度全波场反演方法，并在迭代中利用波场照明进行模型修正，实现海底硫化物矿体由轮廓到细节的高分辨成像。
附图说明
[0057]
图1是近底地震观测方式示意图，使用了低频海面震源、高频深海震源激发，通过深海垂直电缆获取宽频带数据；
[0058]
图2是通过图1所示观测系统获取的宽频带垂直电缆地震数据展示图；
[0059]
图3是本发明进行宽频带垂直电缆地震数据处理的总体流程示意图；
[0060]
图4是本发明压缩感知数据高进度重建结果展示图，缺失道集数据获得了较为准确的重建；
[0061]
图5是本发明多步频率阻尼多尺度波形反演使用的频率域阻尼因子的作用展示图。
[0062]
图6是本发明对深海硫化物矿体进行宽频带地震约束反演成像结果展示图；其中图(a)为真实模型，图(b)为基于多源数据建立的初始速度模型va，图(c)为第一步气枪震源数据全波形反演建立的速度模型(fb＝0hz，fe＝300hz，α＝0.080、0.008、0.003)，图(d)为在图(c)的基础上进行第二步深海震源数据迭代反演建立的速度模型vb(fb＝200hz，fe＝800hz，α＝0.080、0.008、0.003)。
具体实施方式
[0063]
下面结合具体实施方式对本发明做进一步阐述和说明。所述实施例仅是本公开内容的示范且不圈定限制范围。本发明中各个实施方式的技术特征在没有相互冲突的前提下，均可进行相应组合。
[0064]
如图1所示，深海海底硫化物勘探中使用水面气枪震源和深海拖曳换能器声源组成宽频带声源激发声波，前者主要激发150hz以内的低频地震波，而后者主要激发250～800hz的声波，并用靠近目标体的多通道垂直电缆地震电缆接收声学信号。这种观测方式的工作频带较宽，可以提高纵向分辨率；靠近海底激发和接收的工作方式还可以压缩第一菲涅尔带半径，从而提高横向分辨率，进而达到高分辨率勘探的目的。
[0065]
如图2所示，本实施例使用宽频带震源激发-近海底垂直电缆接收的地震数据进行后续分析和处理，该数据具有多方向和宽频带特征，对于实现海底精准成像具有重要意义。
[0066]
如图3所示，本实施例提供的海底硫化物矿体结构宽频带地震约束反演成像方法，包括以下步骤：
[0067]
步骤1：对硫化物矿体样品进行物性测试，获取硫化物矿体有限范围的声速、密度、孔隙度；
[0068]
以已有资料大致确定的硫化物矿体的展布范围为参考，选择在硫化物矿体上不同位置和不同深度的底质取样和钻探岩心样品，开展常压和高压环境的物性测试，获取硫化物矿体有限范围的包括声速、密度、孔隙度在内的物性特征。
[0069]
步骤2：对垂直电缆地震数据进行倾斜状态估计、定位信息融合、直达波压制、数据去噪、振幅恢复、数据反褶积处理；
[0070]
利用超短基线数据确定测线作业的震源和垂直电缆的水下位置，并基于姿态仪数据进行垂直电缆倾斜状态估计，获取各接收点的位置。对近底地震数据进行定位信息融合、直达波压制、数据去噪、振幅恢复、数据反褶积处理。
[0071]
步骤3：对垂直电缆地震数据进行压缩感知高精度重建处理；
[0072]
实际数据采集中，受仪器自身、其他搭载设备影响及海浪等环境条件的影响，气枪和深海声源数据中均有部分炮数据因噪声过大而无法使用。由于数据缺失具有一定的随机性，使用基于压缩感知的方法进行数据高精度重建，可以最大程度上克服数据缺失的问题，满足后续处理精度的要求。由于震源和接收两个方向上均可能有数据缺失，需要在两个方面进行数据的压缩感知重建，对缺失炮和道的数据进行恢复：
[0073]
a)地震采集的数据是成百上千炮数据，对于某些缺失或无效的炮数据，首先进行共接收点道集数据的重建，获取规则共接收点道集数据；
[0074]
b)将上一步重建的数据转换为共炮道集数据，并对每个共炮道集数据进行重建，获取规则的共炮道集数据。
[0075]
如图4所示，上图为含缺失道的垂直电缆地震数据，经过两个步骤就可以对缺失炮和道的数据进行恢复，下图为重建后的规则数据，可以有效提高后续反演处理的精度。
[0076]
所述的步骤3垂直电缆地震数据的压缩感知重建处理中，使用基于迭代阈值的压缩感知稀疏约束重建方法对垂直电缆地震数据进行压缩重建处理。
[0077]
根据压缩感知理论，数据缺失的过程可以表示为：
[0078]
b＝rf
[0079]
其中，f为规则网格完整的地震数据(设共有m道数据)，b为含缺失道的观测地震数据(设共有n道数据，缺失了m-n道数据)，r为缺失采样矩阵(设矩阵维度为n
×
m)。对于地震数据而言，其波形可以使用curvelet变换或者seislet变换进行稀疏表示，设变换矩阵为s，则：
[0080]
b＝ax，with a＝rs
*
[0081]
其中，x为数据f在s域的稀疏表示(维度为m)，a＝rs
*
为测量矩阵，上标*表示复数的共轭。这样数据的重建问题可以表示为：
[0082][0083]
其中，变量上面的波浪线～表示对变量的估计，λ为拉格朗日乘子。
[0084]
由于数据缺失具有随机性，r为随机采样矩阵，在满足稀疏度要求(kd log2(m/n)≤n，k为稀疏度，d为一个常数)的情况下，使用稀疏约束反演方法可以求解上式，获取最终数据重建结果
[0085]
使用迭代阈值方法进行稀疏约束反演，迭代式如下：
[0086][0087]
其中，i为迭代次数，表示阈值运算，γi为迭代步长，θ为使用的阈值。
[0088]
步骤4：多源数据约束初始速度建模；
[0089]
依据步骤1物性测试获取的围岩、底质和岩心的声速特征，结合已有其他方法的探测结果，构建一个平滑速度模型va，平滑速度模型包括海底以上的水体部分和海底以下的固体部分。两部分建模方法如下：
[0090]
1)平滑速度模型的固体部分，其背景速度为围岩速度，矿体范围以已有其他资料(如电法探测)勾勒的硫化物矿体轮廓为基础，块状硫化物使用岩心速度，网脉状硫化物的速度根据已有慢性扩张洋中脊硫化物区两层之间的关系和目标区实测块状硫化物的速度确定。
[0091]
2)平滑速度模型中水体部分，通过在探测区域开展若干站位的ctd探测作业获取
的ctd数据建立水体声速剖面。使用的方法为：以若干站位的ctd获取的一定深度的声速数据作为控制点，通过数据拟合和插值建立整个探测剖面的声速数据。
[0092]
多源数据约束下建立的初始速度模型是对实际速度模型的一种约束和逼近，可以有效提高后续反演的精度。
[0093]
步骤5：多步频率阻尼多尺度全波形反演；
[0094]
全波形反演是通过优化方法寻找目标函数的极小值，满足目标函数最小值的速度即为最终反演结果。利用近底地震多方向和宽频带的特点，分步使用阻尼波场全波形反演方法从低频到高频进行频率域全波形反演，实现由轮廓到细节的逐步重建，建立准确速度模型vb。
[0095]
对于声波方程：
[0096][0097]
其中，v为速度，t为时间，r为地下位置，rs为震源位置，u(r,t)为r处的波场，s(t)δ(r-rs)为震源函数。其全波形反演的目标函数可以表示为：
[0098][0099]
其中，d
obs
是观测数据的逆传播波场，u为根据声波方程计算获取的合成地震记录，t、s、g分别表示时间、震源和检波点位置，频率域的目标函数可以表示为：
[0100][0101]
全波形反演是通过一定的优化方法寻找目标函数e
t
或e
ω
的极小值，满足目标函数最小值的速度即为最终反演结果。对于本实施例中使用海面气枪震源和深海拖曳换能器声源组成宽频带声源的地震数据，多步频率阻尼多尺度全波形反演的过程如下：
[0102]
a)利用近底地震多方向(海面、深海震源多角度联合激发)和宽频带(气枪低频带，深海震源中高频带)的特点，通过对数据的频谱分析可知，海面气枪震源数据有效频带大致为3～150hz，而深海换能器声源的有效频带宽度大致为250～800hz。因此，根据上述认识，确定分两步进行全波形反演处理；
[0103]
b)根据步骤a)的认识，确定第1个频率带的起始频率fb＝0hz，终止频率fe＝300hz，确定第2个频率带的起始频率fb＝200hz，终止频率fe＝800hz；
[0104]
c)对地震数据进行频率域波场阻尼处理，获取每个频率分布带的阻尼波场，计算方法如下：
[0105][0106]
其中，u(f)为频率域波场，γd(f)为频率阻尼算子，f为频率，α为和反演尺度相关的阻尼因子，用于控制反演尺度，0＜α＜1表示波场按照指数曲线压制高频，α＝0表示波场未阻尼。阻尼算子γd如图5所示。
[0107]
d)分两步进行全波形反演，每步对应1个地震数据的频率分布带，且前者反演的结果作为下一步反演的初始速度模型。
[0108]
第1步：针对频率分布带1，对原始地震数据进行α＝0.003、0.008和0.080的阻尼处理获取对应的阻尼波场数据，进而进行反演处理。其中，α＝0.003的阻尼波场数据反演的初始速度为步骤4多源数据约束初始速度建模结果；三组数据全波形反演处理中，前者反演的结果作为后者反演的输入速度模型。阻尼算子γd的曲线如图5a所示。
[0109]
第2步：针对频率分布带2，对原始地震数据进行α＝0.003、0.008和0.080的阻尼处理获取对应的阻尼波场数据，进而进行反演处理。其中，α＝0.003的阻尼波场数据反演的初始速度为上一步全波形反演的结果；三组数据全波形反演处理中，前者反演的结果作为后者反演的输入速度模型。阻尼算子γd的曲线如图5b所示。
[0110]
上述两步反演中，第1步主要是针对气枪震源低频数据，第2步是针对换能器中高频数据，这主要对起始频率fb和终止频率fe和α进行调节，通过多步阻尼迭代反演可以达到预期效果。
[0111]
优选地，在本发明的一个具体实施例中，所述的步骤5多步频率阻尼多尺度全波形反演中，在进行波形反演的过程中增加波场照明，以实现对深部信息的高分辨成像。每次迭代计算一次总的波场照明矩阵，计算公式如下：
[0112][0113]
其中，imi为第i次迭代的波场照明矩阵，vi为第i次迭代生成的速度模型，为模型vi下的第j个震源在位置r处的正演波场，为模型vi下第k个检波点在位置r处的逆传播波场，ns为震源总数，ng为接收点总数。全波形反演中使用以下带有照明信息的模型修改方法：
[0114][0115]
其中，δvi为通过第i次迭代反演出来的模型修改量，模型v
i+1
为第i+1次迭代生成的速度模型。当迭代次数达到预定值或满足预设迭代终止条件时停止迭代，获取最终速度模型vb。预设终止条件为：
[0116][0117]
其中，e(i)为经过i次迭代后目标函数的最小二乘误差，η为一个不小于零的正数。
[0118]
如图6所示，图a为真实速度模型，图b为基于多源数据建立的初始速度模型va，图c为第一步气枪震源数据全波形反演建立的速度模型(fb＝0hz，fe＝300hz，α＝0.080、0.008、0.003)，图d为在图c的基础上进行第二步深海震源数据迭代反演建立的速度模型(fb＝200hz，fe＝800hz，α＝0.080、0.008、0.003)，也就是最终反演结果vb，可见具有较高的分辨率和准确度。
[0119]
以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。对于本领域的普通技术人员来说，在不脱离
本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。