CN116660974A - 一种基于结构耦合约束的体波和面波三维联合反演方法 - Google Patents
一种基于结构耦合约束的体波和面波三维联合反演方法 Download PDFInfo
- Publication number
- CN116660974A CN116660974A CN202310382018.1A CN202310382018A CN116660974A CN 116660974 A CN116660974 A CN 116660974A CN 202310382018 A CN202310382018 A CN 202310382018A CN 116660974 A CN116660974 A CN 116660974A
- Authority
- CN
- China
- Prior art keywords
- model
- wave
- inversion
- travel time
- constraint
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 116
- 230000008878 coupling Effects 0.000 title claims abstract description 17
- 238000010168 coupling process Methods 0.000 title claims abstract description 17
- 238000005859 coupling reaction Methods 0.000 title claims abstract description 17
- 239000006185 dispersion Substances 0.000 claims abstract description 88
- 239000011159 matrix material Substances 0.000 claims description 96
- 238000005259 measurement Methods 0.000 claims description 18
- 238000004364 calculation method Methods 0.000 claims description 14
- 238000009499 grossing Methods 0.000 claims description 11
- 238000010276 construction Methods 0.000 claims description 6
- 238000002224 dissection Methods 0.000 claims description 6
- 238000003384 imaging method Methods 0.000 abstract description 9
- 230000010365 information processing Effects 0.000 abstract description 2
- 230000008569 process Effects 0.000 description 13
- 230000009286 beneficial effect Effects 0.000 description 8
- 230000002159 abnormal effect Effects 0.000 description 6
- 238000001514 detection method Methods 0.000 description 6
- 230000000704 physical effect Effects 0.000 description 6
- 238000010586 diagram Methods 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 4
- 238000012545 processing Methods 0.000 description 4
- 238000011161 development Methods 0.000 description 3
- 238000001914 filtration Methods 0.000 description 3
- 238000011835 investigation Methods 0.000 description 3
- 238000013459 approach Methods 0.000 description 2
- 238000009933 burial Methods 0.000 description 2
- 238000004422 calculation algorithm Methods 0.000 description 2
- 238000005314 correlation function Methods 0.000 description 2
- 238000013461 design Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 229910052500 inorganic mineral Inorganic materials 0.000 description 2
- 230000007246 mechanism Effects 0.000 description 2
- 239000011707 mineral Substances 0.000 description 2
- 238000012544 monitoring process Methods 0.000 description 2
- 238000004613 tight binding model Methods 0.000 description 2
- 238000012952 Resampling Methods 0.000 description 1
- 238000003491 array Methods 0.000 description 1
- 230000000295 complement effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 230000004927 fusion Effects 0.000 description 1
- 230000001939 inductive effect Effects 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 238000000691 measurement method Methods 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 238000003908 quality control method Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 230000035939 shock Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 230000002087 whitening effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/282—Application of seismic models, synthetic seismograms
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/303—Analysis for determining velocity profiles or travel times
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公开了一种基于结构耦合约束的体波和面波三维联合反演方法，属于地震信息处理技术领域，该方法包括获取体波走时数据和面波频散数据；获取网格剖分后的剖分结果；根据剖分结果，分别利用第一初始模型和第二初始模型，构建体波单方法反演模型和面波频散单方法反演模型；分别将第一初始模型和第二初始模型更新为体波单方法反演模型和面波频散单方法反演模型；若达到初始迭代次数，分别根据第一初始模型和第二初始模型构建面波频散约束反演模型和体波约束反演模型，并基于体波约束反演模型和面波频散约束反演模型完成联合反演，否则，重复迭代。本发明克服了传统联合反演两种不同数据体权重因子选择困难的问题，提高了单方法反演成像分辨率。
Description
技术领域
本发明属于地震信息处理技术领域，尤其涉及一种基于结构耦合约束的体波和面波三维联合反演方法。
背景技术
地震成像是研究地球内部结构的最重要工具之一，通过地震成像获取的三维速度模型为油气矿产资源探测及开发、大地震的深部孕震及发震机理、地质灾害监测及预警、地球内部结构及深部动力学过程等领域提供关键信息。随着探测技术不断发展，野外地震观测也逐渐往多参数、智能化、高精度、超密集台阵探测等方向发展，海量的、高精度的地震数据不断积累，高精度成像技术也不断推进。体波和面波是地震记录上两种重要的信号，充分地利用这两种数据体可以有效的提高反演精度。
基于射线理论的反演方法是体波和面波成像的重要手段，因其数据资料处理简单、正反演技术成熟等优势，广泛应用于地球不同尺度介质非均匀性特征研究。目前体波和面波数据反演主要以单方法反演为主，然而单方法反演往往受限于方法本身的局限，成像分辨率受到限制。体波反演精度受限于地震和台站的分布，在地震发生方位不均匀或台站分布稀疏的情况下，成像分辨率有限。而面波反演方法的纵向分辨率随深度增加而降低、对介质P波速度的灵敏度低于S波速度等方面的约束，仅靠面波单方法反演难以获得高精度的速度结构。
已有的体波和面波联合思路是将体波和面波反演过程中的灵敏核函数构建成一个大型线性稀疏矩阵，通过选取两种地震数据的权重等参数，同步求解出P波和S波的速度结构。该联合反演策略是在同一个反演系统中，设定在同一迭代终止条件，联合反演体波走时和面波频散数据，同步获取地下介质的P波和S波速度结构。然而该方法需要选取合适的权重因子，才能保证在同一的评价体系下中获得最优的P波和S波速度结构，确保其联合反演结果优于其单方法反演结果。权重因子的选择不仅仅包含有不同数据体质量控制权重因子、不同数据体之间权重因子、正则化项权重因子等，如何选取合适的权重因子使得联合反演获取的速度结构优于单方法反演结果给该方法带来了很大的困难。
发明内容
针对现有技术中的上述不足，本发明提供的一种基于结构耦合约束的体波和面波三维联合反演方法克服了传统联合反演两种不同数据体权重因子选择困难的问题，提高了单方法反演成像分辨率。
为了达到上述发明目的，本发明采用的技术方案为：一种基于结构耦合约束的体波和面波三维联合反演方法，包括以下步骤：
S1、获取体波走时数据和面波频散数据；
S2、进行网格剖分，得到剖分结果；
S3、构建第一初始模型和第二初始模型；
S4、根据剖分结果对第一初始模型进行体波走时数据反演，得到体波单方法反演模型；
S5、根据剖分结果对第二初始模型进行面波频散数据反演，得到面波频散单方法反演模型；
S6、将第一初始模型更新为体波单方法反演模型，以及将第二初始模型更新为面波频散单方法反演模型，并判断是否达到初始迭代次数，若是，则进入步骤S7，否则，返回步骤S4；
S7、根据第一初始模型构建面波频散约束反演模型，以及根据第二初始模型构建体波约束反演模型，并基于体波约束反演模型和面波频散约束反演模型完成联合反演。
本发明的有益效果为：本发明充分利用了地震数据，并加入双差走时提高了震源区域速度结构精度；利用交叉梯度函数构建体波约束反演模型和面波频散约束反演模型，并在面波频散单方法反演模型的目标函数的构建中带入体波约束反演模型，在体波单方法反演模型的目标函数的构建中带入面波频散约束反演模型，实现两种方法的约束反演，采用交替约束策略联合反演，最终获取在结构上越来越相似的两种速度模型。交替约束联合反演策略可以在联合反演过程中不断更新约束模型，相比直接用固定模型做约束反演，避免了先验约束模型精度对反演结果的影响。
进一步地，所述步骤S2网格剖分具体为中间目标区域网格间距为0.1°，深度网格在26km内网格间距为2km，边界网格采用不均匀剖分。
上述进一步方案的有益效果为：网格划分考虑了探测目标体的规模以及数据覆盖范围，平衡了正演和反演的计算精度。
进一步地，所述步骤S4具体为：
S401、根据剖分结果对第一初始模型进行正演计算，得到体波理论模型正演数据：
其中，为体波理论模型正演数据；τi为地震事件i的发震时刻；i为地震事件编号；k为地震台站编号；u为体波慢度；ds为射线路径微元量；
S402、将体波理论模型正演数据和体波走时数据作差，得到体波实测走时数据与理论模型走时数据的差：
其中，为体波实测走时数据与理论模型走时数据的差，表示第i个地震事件和第k个地震台站之间的体波走时差；q为震源发震位置参数编号；/>为震源参数扰动矩阵；xq为震源发震位置参数；Δ为拉普拉斯算子；/>为偏导数符号；Cik为地震事件i到地震台站k的传播路径；N为网格总量；l为网格编号；Δul为第l个网格的慢度扰动；
S403、根据体波实测走时数与理论模型走时数据的差，获取所有地震事件及其对应的接收地震台站的体波走时差，得到数据走时残差矩阵：
ΔT＝AΔX+CΔU
其中，ΔT为数据走时残差矩阵；A为所有震源参数的扰动组合成的大型矩阵；ΔX为所有震源的空间扰动位置组成的一维矩阵；C为所有地震事件传播路径矩阵；ΔU为所有网格的慢度扰动组成的一维矩阵；
S404、在数据走时残差矩阵中加入双差矩阵和圆滑因子，得到体波单方法反演模型，所述体波单方法反演模型的目标函数为：
其中，L1为体波单方法反演模型的目标函数，表示体波单方法所有数据走时差；α和β均为权重因子；Wm为圆滑因子；η为圆滑权重系数；QDD为双差矩阵；0为零矩阵。
上述进一步方案的有益效果为：将通过不断修正模型，使得体波理论模型正演数据与体波走时数据越来越逼近，获取满足要求的速度结构，并加入双差走时，提高震源区域速度结构的精度。
进一步地，所述双差矩阵具体为双差走时组成的矩阵；所述双差走时为两个位置相近的震源到同一观测台的实测走时和理论模型走时的时差做差。
上述进一步方案的有益效果为：提高了震源区域速度结构的精度。
进一步地，所述步骤S5具体为：
S501、根据剖分结果对第二初始模型进行正演计算，得到面波理论模型正演数据：
其中，为面波理论模型正演数据；τi为地震事件i的发震时刻；i为地震事件编号；k为地震台站编号；v为面波慢度；ds为射线路径微元量；
S502、将面波理论模型正演数据和面波频散数据作差，得到面波实测频散走时数据和理论模型正演走时数据的差：
其中，δtj(ω)为面波实测频散走时数据和理论模型正演走时数据的差，即面波走时差；j为射线编号；ω为频率；ΔLb为射线在第b个网格的长度；b为网格编号；Cb(ω)为地表第b个网格频散速度；为地表第b个网格在深度za处纵波速度的偏导数；za为地表第a个网格的z深度处；p为P波速度；a为网格编号；L为网格个数；Rp为P波速度的系数矩阵；Rρ为密度的系数矩阵；δρb(za)为第b个网格在深度za处纵波速度的偏导数；δβb(za)为第b个网格在深度za处横波速度的扰动；Δ为拉普拉斯算子；
S503、根据面波实测频散走时数据和理论模型正演走时数据的差，获取所有地震事件及其对应的接收地震台站的面波走时差，得到面波频散走时差矩阵：
ΔdSW＝GSWΔms
其中，ΔdSW为面波频散走时差矩阵；GSW为灵敏核函数矩阵；Δms为待求解物性模型；
S504、在面波频散走时差矩阵中加入模型光滑项，得到面波频散单方法反演模型，所述面波频散单方法反演模型的目标函数为：
其中，L2为面波频散单方法反演模型的目标函数，表示面波频散单方法所有数据走时差；WSW为模型光滑项。
上述进一步方案的有益效果为：将通过不断修正模型，使得面波理论模型正演数据与面波频散数据越来越逼近，获取满足要求的速度结构，在走时差矩阵中加入模型光滑项，提高反演结果的稳定性。
进一步地，所述步骤S7具体为：
S701、构建并根据交叉梯度函数，分别利用第一初始模型和第二初始模型作为约束构建面波频散约束反演目标函数和体波约束反演目标函数；
S702、根据面波频散约束反演目标函数和体波约束反演目标函数分别得到面波频散约束反演模型和体波约束反演模型；
S703、采用交替迭代策略进行联合约束反演，判断是否满足迭代终止条件，若是，则得到P波速度模型和S波速度模型，否则，将第一初始模型更新为体波约束反演模型，以及将第二初始模型更新为面波频散约束反演模型，并返回步骤S701。
上述进一步方案的有益效果为：利用体波约束反演模型和面波频散约束反演模型，在体波单方法反演模型的目标函数中带入面波频散约束反演模型，实现体波约束反演，以及在面波频散单方法反演模型的目标函数中带入体波约束反演模型，实现面波频散约束反演，采用交替约束反演策略联合反演体波走时数据和面波频散数据，交替约束联合反演后获取的两种模型在结构上也越来越相似，相比直接用固定模型做约束反演对约束模型的依赖性小一些，避免了先验约束模型精度对反演结果的影响。
进一步地，所述步骤S701具体为：
S7011、构建交叉梯度函数：
其中，D为交叉梯度函数；Wcg,x为交叉梯度函数x分量的偏导数矩阵；Wcg,y为交叉梯度函数y分量的偏导数矩阵；Wcg,z为交叉梯度函数z分量的偏导数矩阵；ΔM为待求解的模型；t0,x为交叉梯度函数x分量；t0,y为交叉梯度函数y分量；t0,z为交叉梯度函数z分量；ΔU为所有网格的慢度扰动组成的一维矩阵；U为慢度；
S7012、根据体波单方法反演模型的目标函数和交叉梯度函数，以第二初始模型为约束构建体波约束反演目标函数：
其中，L3为体波约束反演模型目标函数；A为所有震源参数的扰动组合成的大型矩阵；QDD为双差矩阵；ζ、α和β均为权重因子；η为圆滑权重系数；C为所有地震事件传播路径矩阵；0为零矩阵；ΔT为数据走时残差矩阵；ΔX为所有震源的空间扰动位置组成的一维矩阵；
S7013、根据面波频散单方法反演模型和交叉梯度函数，以第一初始模型为约束构建面波频散约束反演目标函数：
其中，L4为面波频散约束反演目标函数；ΔdSW为走时差矩阵；μ为交叉梯度约束的权重因子；GSW为灵敏核函数矩阵；WSW为模型光滑项；Δms为待求解物性模型。
上述进一步方案的有益效果为：利用交叉梯度函数构建体波约束反演目标函数和面波频散约束反演目标函数，并在体波单方法反演模型的目标函数中带入面波频散约束反演目标函数，在面波频散单方法反演模型的目标函数中带入体波约束反演目标函数，采用交替约束策略联合反演体波走时数据和面波频散数据，经交替约束联合反演后获取的两种模型在结构上也越来越相似，相比直接用固定模型做约束反演对约束模型的依赖性小一些，避免了先验约束模型精度对反演结果的影响。
进一步地，所述步骤S704中迭代终止条件为满足数据拟合差的最小值，或相邻两次迭代的第一初始模型和第二初始模型之间的值差满足设定的阈值，或满足设定的最大迭代次数。
上述进一步方案的有益效果为：设定迭代终止条件，为联合反演提供结束出口，保证最终P波速度模型和S波速度模型的获取。
附图说明
图1为本发明的方法流程图。
图2为本发明提供的三维正、反演网格剖分示意图。
图3为本发明实施例提供的三维理论模型示意图。
图4为本发明实施例提供的联合反演后P波速度和S波速度交叉梯度分布图。
具体实施方式
下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。
如图1所示，在本发明的一个实施例中，一种基于结构耦合约束的体波和面波三维联合反演方法，包括以下步骤：
S1、获取体波走时数据和面波频散数据；
S2、进行网格剖分，得到剖分结果；
S3、构建第一初始模型和第二初始模型；
S4、根据剖分结果对第一初始模型进行体波走时数据反演，得到体波单方法反演模型；
S5、根据剖分结果对第二初始模型进行面波频散数据反演，得到面波频散单方法反演模型；
S6、将第一初始模型更新为体波单方法反演模型，以及将第二初始模型更新为面波频散单方法反演模型，并判断是否达到初始迭代次数，若是，则进入步骤S7，否则，返回步骤S4；
S7、根据第一初始模型构建面波频散约束反演模型，以及根据第二初始模型构建体波约束反演模型，并基于体波约束反演模型和面波频散约束反演模型完成联合反演。
本实施例中，体波走时数据可以从区域和地方地震台网、流动密集台阵的连续地震波形数据中获取。这两种数据来源对体波走时的提取方式不一样。对于区域和地方地震台网观测数据，选取观测报告中震级大等于1.0且走时残差小等于0.5s的震相数据。对于密集流动地震台阵观测数据，采用统软件拾取波形初至。将所有走时数据画时距曲线，剔除走时曲线中离散较大的震相，确保每个地震至少有5个台站记录到P波初至到时。
本实施例中，面波处理所需的连续地震波形资料来源与体波走时资料相同，但处理手段不一样。具体为：对连续数据重采样、去仪器响应、去均值、去趋势、滤波；对时间序列进行谱白化和时间归一化处理；将连续信号按天或者小时分割，对相同时间段数据进行互相关计算，获取该时段台站对之间的互相关函数。然后将同一台站对不同时间段的互相关函数进行叠加，获取该台站对的经验格林函数；对经验格林函数进行滤波、归一化、时间窗选取等处理。采用多窗口测量法获取理论格林函数和经验格林函数不同频段的走时残差。
所述步骤S2网格剖分具体为中间目标区域网格间距为0.1°，深度网格在26km内网格间距为2km，边界网格采用不均匀剖分。
本实施例中，网格剖分是正、反演计算精度的关键因素之一。密集网格可以提高正演精度，但会增加反演的多解性。网格稀疏会导致正演精度降低，从而影响反演精度。同时网格划分还需要重点考虑探测目标体的规模以及数据覆盖范围。图2为本发明提供的三维正、反演网格剖分示意图，同时考虑了台站、地震分布、分辨率等需求进行网格剖分，在网格剖分方面中间目标区网格间距为0.1°，深度网格在26km内网格间距为2km，边界网格采用不均匀剖分。
所述步骤S4具体为：
S401、根据剖分结果对第一初始模型进行正演计算，得到体波理论模型正演数据：
其中，为体波理论模型正演数据；τi为地震事件i的发震时刻；i为地震事件编号；k为地震台站编号；u为体波慢度；ds为射线路径微元量；
S402、将体波理论模型正演数据和体波走时数据作差，得到体波实测走时数据与理论模型走时数据的差：
其中，为体波实测走时数据与理论模型走时数据的差，表示第i个地震事件和第k个地震台站之间的体波走时差；q为震源发震位置参数编号；/>为震源参数扰动矩阵；xq为震源发震位置参数；Δ为拉普拉斯算子；/>为偏导数符号；Cik为地震事件i到地震台站k的传播路径；N为网格总量；l为网格编号；Δul为第l个网格的慢度扰动；
S403、根据体波实测走时数与理论模型走时数据的差，获取所有地震事件及其对应的接收地震台站的体波走时差，得到数据走时残差矩阵：
ΔT＝AΔX+CΔU
其中，ΔT为数据走时残差矩阵；A为所有震源参数的扰动组合成的大型矩阵；ΔX为所有震源的空间扰动位置组成的一维矩阵；C为所有地震事件传播路径矩阵；ΔU为所有网格的慢度扰动组成的一维矩阵；
S404、在数据走时残差矩阵中加入双差矩阵和圆滑因子，得到体波单方法反演模型，所述体波单方法反演模型的目标函数为：
其中，L1为体波单方法反演模型的目标函数，表示体波单方法所有数据走时差；α和β均为权重因子；Wm为圆滑因子；η为圆滑权重系数；QDD为双差矩阵；0为零矩阵。
本实施例中，正演计算是本项目开发的算法的重要过程之一。在理论模型测试中，需要通过正演计算获取模型的正演数据，并将理论模型正演数据与观测数据做差，带入反演计算流程。在反演计算中需要通过计算每次反演模型的正演数据，计算理论数据和观测数据的走时差，通过走时差用于评估反演效果。以下简要介绍体波和面波正演计算过程。基于射线理论的走时成像正演计算理论可简化成射线从源(虚拟源)到地震台站的旅行事件。体波和面波走时的计算区别在于慢度的计算。体波慢度的获取是介质的速度的倒数，通过直接获取介质的速度物性即可。而面波慢度的计算需要对面波进行分频段滤波，分别对每个频段计算地表网格面波传播速度，对地表网格速度求倒数即可获得面波慢度。
本实施例中，采用的地震体波主要针对的是近震体波，即震源在研究区域之内。对于震源在研究区外的事件，可以通过射线追踪到研究区域内的边界网格上，将边界网格交点作为源，然后当成区域内的震源进行处理。
本实施例中，基于射线理论的体波反演问题跟正演问题中的已知量和未知量刚好相反，在反演过程中，采用一阶平滑、二阶平滑或高斯平滑等方式进行模型光滑可以提高反演的稳定性和反演结果的可靠性；对于每个地震信息的ΔX(所有震源的空间扰动位置组成的一维矩阵)是相互独立的，不能进行光滑处理，此处处理方式为补0。
所述双差矩阵具体为双差走时组成的矩阵；所述双差走时为两个位置相近的震源到同一观测台的实测走时和理论模型走时的时差做差。
所述步骤S5具体为：
S501、根据剖分结果对第二初始模型进行正演计算，得到面波理论模型正演数据：
其中，为面波理论模型正演数据；τi为地震事件i的发震时刻；i为地震事件编号；k为地震台站编号；v为面波慢度；ds为射线路径微元量；
S502、将面波理论模型正演数据和面波频散数据作差，得到面波实测频散走时数据和理论模型正演走时数据的差：
其中，δtj(ω)为面波实测频散走时数据和理论模型正演走时数据的差，即面波走时差；j为射线编号；ω为频率；ΔLb为射线在第b个网格的长度；b为网格编号；Cb(ω)为地表第b个网格频散速度；为地表第b个网格在深度za处纵波速度的偏导数；za为地表第a个网格的z深度处；p为P波速度；a为网格编号；L为网格个数；Rp为P波速度的系数矩阵；Rρ为密度的系数矩阵；δρb(za)为第b个网格在深度za处纵波速度的偏导数；δβb(za)为第b个网格在深度za处横波速度的扰动；Δ为拉普拉斯算子；
S503、根据面波实测频散走时数据和理论模型正演走时数据的差，获取所有地震事件及其对应的接收地震台站的面波走时差，得到面波频散走时差矩阵：
ΔdSW＝GSWΔms
其中，ΔdSW为面波频散走时差矩阵；GSW为灵敏核函数矩阵；Δms为待求解物性模型；
S504、在面波频散走时差矩阵中加入模型光滑项，得到面波频散单方法反演模型，所述面波频散单方法反演模型的目标函数为：
其中，L2为面波频散单方法反演模型的目标函数，表示面波频散单方法所有数据走时差；WSW为模型光滑项。
本实施例中，模型光滑项可以采用一阶、二阶、三阶、高斯等平滑公式，模型光滑项可以提高反演的稳定性和反演结果的可靠性。
所述步骤S7具体为：
S701、构建并根据交叉梯度函数，分别利用第一初始模型和第二初始模型作为约束构建面波频散约束反演目标函数和体波约束反演目标函数；
S702、根据面波频散约束反演目标函数和体波约束反演目标函数分别得到面波频散约束反演模型和体波约束反演模型；
S703、采用交替迭代策略进行联合约束反演，判断是否满足迭代终止条件，若是，则得到P波速度模型和S波速度模型，否则，将第一初始模型更新为体波约束反演模型，以及将第二初始模型更新为面波频散约束反演模型，并返回步骤S701。
所述步骤S701具体为：
S7011、构建交叉梯度函数：
其中，D为交叉梯度函数；Wcg,x为交叉梯度函数x分量的偏导数矩阵；Wcg,y为交叉梯度函数y分量的偏导数矩阵；Wcg,z为交叉梯度函数z分量的偏导数矩阵；ΔM为待求解的模型；t0,x为交叉梯度函数x分量；t0,y为交叉梯度函数y分量；t0,z为交叉梯度函数z分量；ΔU为所有网格的慢度扰动组成的一维矩阵；U为慢度；
S7012、根据体波单方法反演模型的目标函数和交叉梯度函数，以第二初始模型为约束构建体波约束反演目标函数：
其中，L3为体波约束反演模型目标函数；A为所有震源参数的扰动组合成的大型矩阵；QDD为双差矩阵；ζ、α和β均为权重因子；η为圆滑权重系数；C为所有地震事件传播路径矩阵；0为零矩阵；ΔT为数据走时残差矩阵；ΔX为所有震源的空间扰动位置组成的一维矩阵；
S7013、根据面波频散单方法反演模型和交叉梯度函数，以第一初始模型为约束构建面波频散约束反演目标函数：
其中，L4为面波频散约束反演目标函数；ΔdSW为走时差矩阵；μ为交叉梯度约束的权重因子；GSW为灵敏核函数矩阵；WSW为模型光滑项；Δms为待求解物性模型。
本实施例中，从Gallardo和Meju(2003)提出的不同物性参数的结构耦合关系式出发，所述结构耦合关系式为：
其中，ms和mp分别为S波速度和P波速度模型，为这两种物性的交叉梯度向量，tx为交叉梯度函数的x分量，ty为交叉梯度函数的y分量，tz为交叉梯度函数的z分量，m0为第一个初始模型和第二个初始模型。将结构耦合关系式在m0处进行泰勒级数展开取一阶偏导数，得到公式(1)：
公式(1)中m可以是ms或mp，根据交叉梯度函数离散化公式，可以将公式(1)写成矩阵形式，得到公式(2)：
其中，Wcg为交叉梯度函数系数矩阵，ΔM为m-m0。将公式(2)写成三份量形式，得到公式(3)：
公式(3)即为交叉梯度函数离散化的线性方程组。本发明将公式(3)加入体波和面波单方法反演目标函数中，即可实现约束反演的目标函数构建.
所述步骤S704中迭代终止条件为满足数据拟合差的最小值，或相邻两次迭代的第一初始模型和第二初始模型之间的值差满足设定的阈值，或满足设定的最大迭代次数。
本实施例中，本发明采用交替迭代的策略实现体波走时和面波频散的联合反演。该策略是以约束反演目标函数为主线，在两种方法迭代过程中，将两种方法反演迭代过程中获取的更新模型作为约束模型去约束另一种方法，直至满足迭代终止条件。本发明迭代终止条件主要采用三种判断：(1)满足数据拟合差的最小值；(2)相邻两次迭代的模型值差满足满足设定的阈值；(3)满足设定的最大迭代次数。在迭代过程中，随着每种方法获取的模型越来越接近于真实模型，带入另一种方法的约束模型也越来越可靠，交替约束联合反演后获取的两种模型在结构上也越来越相似，相比直接用固定模型做约束反演对约束模型的依赖性小一些，避免先验约束模型精度对反演结果的影响。
本实施例中，为了考察本发明的算法对不同规模、埋深、异常强度的异常体的分辨率能力，本项目设计多种理论物性模型(包括规模不同、异常强度不同、埋深不同的单个或多个异常体)，通过理论模型合成数据进行反演试算，理论模型即为第一初始模型和第二初始模型的原始状态，所述理论模型极为地震波的模拟模型，第一初始模型和第二初始模型未开始反演前均为理论模型。图3为本发明实施例提供的三维理论模型示意图。该模型设计思路是为了充分发挥体波和面波的优势，浅部同一深度设计了4个规模一样的低速异常体(异常值低于背景速度3％)，用于考察两种方法横向分辨率的差异。在四个异常体下方设计了一个规模较大的高速异常体(异常值高于背景速度3％)，用于考察两种方法的垂向分辨率的差异。
本实施例中，由于反演迭代开始时，P波速度和S波速度模型跟初始模型相近，交叉梯度约束的作用小，因此刚开始迭代的时候，先进行地震体波走时和面波频散的单方法反演，当经过一定次数迭代后，启动两种方法的约束反演模块，将单方法反演的结果作为约束，约束另一种方法，在迭代更新过程中，两种方法的约束模型交替更新，实现两种方法的约束联合反演。
本实施例中，图4为联合反演后P波速度模型和S波速度模型的交叉梯度分布图，图中显示了交叉梯度值主要集中在这两种物性边界区域，这也间接说明在约束反演过程中，边界信息会通过交叉梯度项影响目标函数，从而影响反演结果，获得两种物性结构上更相似的模型。在一些结构不相似或物性没有变换的区域交叉梯度值几乎为零，在反演过程中，就不会影响目标函数的求解，从而保留了单方法的优势。这也是基于交叉梯度实现联合反演的目的，在联合反演过程中寻求两种物性相似的区域，使其结构更为耦合，最终获取形态上更为相似的多种物理模型。
本发明通过提供了一种地震体波走时和面波频散的联合反演方法，该方法不仅可以充分利用已有的地震数据，发挥不同地震数据体的优势，获取更高精度的速度结构，推进多种地震数据融合的高精度反演技术。同时其获取的高精度P波和S波速度模型还可以为油气矿产资源探测及开发、大地震的深部孕震及发震机理、地质灾害监测及预警、地球内部结构及深部动力学过程等领域提供关键信息。
Claims (8)
1.一种基于结构耦合约束的体波和面波三维联合反演方法，其特征在于，包括以下步骤：
S1、获取体波走时数据和面波频散数据；
S2、进行网格剖分，得到剖分结果；
S3、构建第一初始模型和第二初始模型；
S4、根据剖分结果对第一初始模型进行体波走时数据反演，得到体波单方法反演模型；
S5、根据剖分结果对第二初始模型进行面波频散数据反演，得到面波频散单方法反演模型；
S6、将第一初始模型更新为体波单方法反演模型，以及将第二初始模型更新为面波频散单方法反演模型，并判断是否达到初始迭代次数，若是，则进入步骤S7，否则，返回步骤S4；
S7、根据第一初始模型构建面波频散约束反演模型，以及根据第二初始模型构建体波约束反演模型，并基于体波约束反演模型和面波频散约束反演模型完成联合反演。
2.根据权利要求1所述基于结构耦合约束的体波和面波三维联合反演方法，其特征在于，所述步骤S2网格剖分具体为中间目标区域网格间距为0.1°，深度网格在26km内网格间距为2km，边界网格采用不均匀剖分。
3.根据权利要求1所述基于结构耦合约束的体波和面波三维联合反演方法，其特征在于，所述步骤S4具体为：
S401、根据剖分结果对第一初始模型进行正演计算，得到体波理论模型正演数据：
其中，为体波理论模型正演数据；τi为地震事件i的发震时刻；i为地震事件编号；k为地震台站编号；u为体波慢度；ds为射线路径微元量；
S402、将体波理论模型正演数据和体波走时数据作差，得到体波实测走时数据与理论模型走时数据的差：
其中，为体波实测走时数据与理论模型走时数据的差，表示第i个地震事件和第k个地震台站之间的体波走时差；q为震源发震位置参数编号；/>为震源参数扰动矩阵；xq为震源发震位置参数；Δ为拉普拉斯算子；/>为偏导数符号；Cik为地震事件i到地震台站k的传播路径；N为网格总量；l为网格编号；Δul为第l个网格的慢度扰动；
S403、根据体波实测走时数与理论模型走时数据的差，获取所有地震事件及其对应的接收地震台站的体波走时差，得到数据走时残差矩阵：
ΔT＝AΔX+CΔU
其中，ΔT为数据走时残差矩阵；A为所有震源参数的扰动组合成的大型矩阵；ΔX为所有震源的空间扰动位置组成的一维矩阵；C为所有地震事件传播路径矩阵；ΔU为所有网格的慢度扰动组成的一维矩阵；
S404、在数据走时残差矩阵中加入双差矩阵和圆滑因子，得到体波单方法反演模型，所述体波单方法反演模型的目标函数为：
其中，L1为体波单方法反演模型的目标函数，表示体波单方法所有数据走时差；α和β均为权重因子；Wm为圆滑因子；η为圆滑权重系数；QDD为双差矩阵；0为零矩阵。
4.根据权利要求3所述基于结构耦合约束的体波和面波三维联合反演方法，其特征在于，所述双差矩阵具体为双差走时组成的矩阵；所述双差走时为两个位置相近的震源到同一观测台的实测走时和理论模型走时的时差做差。
5.根据权利要求1所述基于结构耦合约束的体波和面波三维联合反演方法，其特征在于，所述步骤S5具体为：
S501、根据剖分结果对第二初始模型进行正演计算，得到面波理论模型正演数据；
其中，为面波理论模型正演数据；τi为地震事件i的发震时刻；i为地震事件编号；k为地震台站编号；v为面波慢度；ds为射线路径微元量；
S502、将面波理论模型正演数据和面波频散数据作差，得到面波实测频散走时数据和理论模型正演走时数据的差：
其中，δtj(ω)为面波实测频散走时数据和理论模型正演走时数据的差，即面波走时差；j为射线编号；ω为频率；ΔLb为射线在第b个网格的长度；b为网格编号；Cb(ω)为地表第b个网格频散速度；为地表第b个网格在深度za处纵波速度的偏导数；za为地表第a个网格的z深度处；p为P波速度；a为网格编号；L为网格个数；Rp为P波速度的系数矩阵；Rρ为密度的系数矩阵；δρb(za)为第b个网格在深度za处纵波速度的偏导数；δβb(za)为第b个网格在深度za处横波速度的扰动；Δ为拉普拉斯算子；
S503、根据面波实测频散走时数据和理论模型正演走时数据的差，获取所有地震事件及其对应的接收地震台站的面波走时差，得到面波频散走时差矩阵：
ΔdSW＝GSWΔms
其中，ΔdSW为面波频散走时差矩阵；GSW为灵敏核函数矩阵；Δms为待求解物性模型；
S504、在面波频散走时差矩阵中加入模型光滑项，得到面波频散单方法反演模型，所述面波频散单方法反演模型的目标函数为：
其中，L2为面波频散单方法反演模型的目标函数，表示面波频散单方法所有数据走时差；WSW为模型光滑项。
6.根据权利要求1所述基于结构耦合约束的体波走时和面波频散的三维联合反演方法，其特征在于，所述步骤S7具体为：
S701、构建并根据交叉梯度函数，分别利用第一初始模型和第二初始模型作为约束构建面波频散约束反演目标函数和体波约束反演目标函数；
S702、根据面波频散约束反演目标函数和体波约束反演目标函数分别得到面波频散约束反演模型和体波约束反演模型；
S703、采用交替迭代策略进行联合约束反演，判断是否满足迭代终止条件，若是，则得到P波速度模型和S波速度模型，否则，将第一初始模型更新为体波约束反演模型，以及将第二初始模型更新为面波频散约束反演模型，并返回步骤S701。
7.根据权利要求6所述基于结构耦合约束的体波走时和面波频散三维联合反演方法，其特征在于，所述步骤S701具体为：
S7011、构建交叉梯度函数：
其中，D为交叉梯度函数；Wcg,x为交叉梯度函数x分量的偏导数矩阵；Wcg,y为交叉梯度函数y分量的偏导数矩阵；Wcg,z为交叉梯度函数z分量的偏导数矩阵；ΔM为待求解的模型；t0,x为交叉梯度函数x分量；t0,y为交叉梯度函数y分量；t0,z为交叉梯度函数z分量；ΔU为所有网格的慢度扰动组成的一维矩阵；U为慢度；
S7012、根据体波单方法反演模型的目标函数和交叉梯度函数，以第二初始模型为约束构建体波约束反演目标函数：
其中，L3为体波约束反演模型目标函数；A为所有震源参数的扰动组合成的大型矩阵；QDD为双差矩阵；ζ、α和β均为权重因子；η为圆滑权重系数；C为所有地震事件传播路径矩阵；0为零矩阵；ΔT为数据走时残差矩阵；ΔX为所有震源的空间扰动位置组成的一维矩阵；
S7013、根据面波频散单方法反演模型和交叉梯度函数，以第一初始模型为约束构建面波频散约束反演目标函数：
其中，L4为面波频散约束反演目标函数；ΔdSW为走时差矩阵；μ为交叉梯度约束的权重因子；GSW为灵敏核函数矩阵；WSW为模型光滑项；Δms为待求解物性模型。
8.根据权利要求6所述基于结构耦合约束的体波走时和面波频散三维联合反演方法，其特征在于，所述步骤S704中迭代终止条件为满足数据拟合差的最小值，或相邻两次迭代的第一初始模型和第二初始模型之间的值差满足设定的阈值，或满足设定的最大迭代次数。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310382018.1A CN116660974B (zh) | 2023-04-11 | 2023-04-11 | 一种基于结构耦合约束的体波和面波三维联合反演方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310382018.1A CN116660974B (zh) | 2023-04-11 | 2023-04-11 | 一种基于结构耦合约束的体波和面波三维联合反演方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116660974A true CN116660974A (zh) | 2023-08-29 |
| CN116660974B CN116660974B (zh) | 2024-02-23 |
Family
ID=87723097
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310382018.1A Active CN116660974B (zh) | 2023-04-11 | 2023-04-11 | 一种基于结构耦合约束的体波和面波三维联合反演方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116660974B (zh) |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20110134722A1 (en) * | 2009-12-07 | 2011-06-09 | Mediwound Ltd. | Simultaneous Joint Inversion of Surface Wave and Refraction Data |
| US20140136170A1 (en) * | 2011-07-21 | 2014-05-15 | Exxonmobile Upstream Research Company | Adaptive Weighting of Geophysical Data Types in Joint Inversion |
| CN106772575A (zh) * | 2016-11-28 | 2017-05-31 | 安徽理工大学 | 一种基于折射波与面波联合反演剩余煤层厚度的方法 |
| CN109061731A (zh) * | 2018-09-18 | 2018-12-21 | 中国地震局地壳应力研究所 | 面波频散与体波谱比联合反演浅层速度的全局优化方法 |
| US20210141113A1 (en) * | 2018-08-06 | 2021-05-13 | Southern University Of Science And Technology | Active source surface wave prospecting method, surface wave exploration device and computer-readable storage medium |
| US20220326403A1 (en) * | 2020-05-12 | 2022-10-13 | Shandong University | Multi-wavefield seismic detection method and system based on construction noise of shield machine |
| CN115508908A (zh) * | 2022-09-27 | 2022-12-23 | 成都理工大学 | 一种地震面波走时和重力异常联合反演方法与系统 |
-
2023
- 2023-04-11 CN CN202310382018.1A patent/CN116660974B/zh active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20110134722A1 (en) * | 2009-12-07 | 2011-06-09 | Mediwound Ltd. | Simultaneous Joint Inversion of Surface Wave and Refraction Data |
| US20140136170A1 (en) * | 2011-07-21 | 2014-05-15 | Exxonmobile Upstream Research Company | Adaptive Weighting of Geophysical Data Types in Joint Inversion |
| CN106772575A (zh) * | 2016-11-28 | 2017-05-31 | 安徽理工大学 | 一种基于折射波与面波联合反演剩余煤层厚度的方法 |
| US20210141113A1 (en) * | 2018-08-06 | 2021-05-13 | Southern University Of Science And Technology | Active source surface wave prospecting method, surface wave exploration device and computer-readable storage medium |
| CN109061731A (zh) * | 2018-09-18 | 2018-12-21 | 中国地震局地壳应力研究所 | 面波频散与体波谱比联合反演浅层速度的全局优化方法 |
| US20220326403A1 (en) * | 2020-05-12 | 2022-10-13 | Shandong University | Multi-wavefield seismic detection method and system based on construction noise of shield machine |
| CN115508908A (zh) * | 2022-09-27 | 2022-12-23 | 成都理工大学 | 一种地震面波走时和重力异常联合反演方法与系统 |
Non-Patent Citations (2)
| Title |
|---|
| E. M. GOLOS ET AL.: "Shear Wave Tomography Beneath the United States Using a Joint Inversion of Surface and Body Waves", JOURNAL OF GEOPHYSICAL RESEARCH: SOLID EARTH * |
| 吴萍萍 等: "基于交叉梯度约束的电阻率法和背景噪声法三维联合反演研究", 地球物理学报, vol. 63, no. 10 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116660974B (zh) | 2024-02-23 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107505654B (zh) | 基于地震记录积分的全波形反演方法 | |
| CN104597490A (zh) | 基于精确Zoeppritz方程的多波AVO储层弹性参数反演方法 | |
| CN104614763A (zh) | 基于反射率法的多波avo储层弹性参数反演方法及系统 | |
| CN104459784B (zh) | 基于单台、双台和双事件数据二维Lg波Q值层析成像方法 | |
| CN104122585A (zh) | 基于弹性波场矢量分解与低秩分解的地震正演模拟方法 | |
| CN105549068A (zh) | 一种三维各向异性微地震干涉逆时定位方法及系统 | |
| CN109100792B (zh) | 基于台站与三维地震联合采集资料的速度反演方法 | |
| CN109459787B (zh) | 基于地震槽波全波形反演的煤矿井下构造成像方法及系统 | |
| CN105242328B (zh) | 古热岩石圈厚度的确定方法及装置 | |
| Zhang et al. | Static corrections in mountainous areas using Fresnel-wavepath tomography | |
| CN104360396A (zh) | 一种海上井间tti介质三种初至波走时层析成像方法 | |
| US10132945B2 (en) | Method for obtaining estimates of a model parameter so as to characterise the evolution of a subsurface volume | |
| CN113671570B (zh) | 一种地震面波走时和重力异常联合反演方法与系统 | |
| CN113917560B (zh) | 一种三维重磁电震多参数协同反演方法 | |
| CN111665556B (zh) | 地层声波传播速度模型构建方法 | |
| CN116466402B (zh) | 一种基于地质信息和电磁数据联合驱动的电磁反演方法 | |
| Waheed et al. | A holistic approach to computing first-arrival traveltimes using neural networks | |
| CN116660974B (zh) | 一种基于结构耦合约束的体波和面波三维联合反演方法 | |
| CN103217715B (zh) | 多尺度规则网格层析反演静校正方法 | |
| CN111025388B (zh) | 一种多波联合的叠前波形反演方法 | |
| CN112100906A (zh) | 数据驱动的大尺度密度建模方法、计算设备及存储介质 | |
| Zunino et al. | Integrating gradient information with probabilistic traveltime tomography using the Hamiltonian Monte Carlo algorithm | |
| CN115373022A (zh) | 一种基于振幅相位校正的弹性波场Helmholtz分解方法 | |
| Wang et al. | Rock Fracture Monitoring Based on High‐Precision Microseismic Event Location Using 3D Multiscale Waveform Inversion | |
| CN110857999B (zh) | 一种基于全波形反演的高精度波阻抗反演方法及系统 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |