CN117555025B - 一种基于重力数据的多层地壳结构反演方法 - Google Patents

一种基于重力数据的多层地壳结构反演方法 Download PDF

Info

Publication number
CN117555025B
CN117555025B CN202410039732.5A CN202410039732A CN117555025B CN 117555025 B CN117555025 B CN 117555025B CN 202410039732 A CN202410039732 A CN 202410039732A CN 117555025 B CN117555025 B CN 117555025B
Authority
CN
China
Prior art keywords
crust
depth
inversion
data
gravity
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.)
Active
Application number
CN202410039732.5A
Other languages
English (en)
Other versions
CN117555025A (zh
Inventor
田甜
姜文亮
何仲太
龚丽霞
李强
张景发
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
National Institute of Natural Hazards
Original Assignee
National Institute of Natural Hazards
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by National Institute of Natural Hazards filed Critical National Institute of Natural Hazards
Priority to CN202410039732.5A priority Critical patent/CN117555025B/zh
Publication of CN117555025A publication Critical patent/CN117555025A/zh
Application granted granted Critical
Publication of CN117555025B publication Critical patent/CN117555025B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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. analysis, for interpretation, for correction
    • G01V1/30—Analysis
    • G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
    • 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. analysis, for interpretation, for correction
    • 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. analysis, for interpretation, for correction
    • G01V1/30—Analysis
    • G01V1/301—Analysis for determining seismic cross-sections or geostructures
    • G01V1/302—Analysis for determining seismic cross-sections or geostructures in 3D data cubes
    • 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. analysis, for interpretation, for correction
    • G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Abstract

本发明公开了一种基于重力数据的多层地壳结构反演方法，涉及地壳结构反演领域，本发明基于布格重力异常数据，通过多尺度小波分解，将重力异常分解为不同深度的信号，寻找对应深度，进行层面反演，利用反演结果建立三维地壳分层初始模型；为初始模型的分层进行密度参数赋值，观察残差，调整密度参数，固定密度赋值；在固定密度参数赋值的基础上，对模型进行正演，根据正演残差迭代调整各层深度，实现地壳多层面几何形态反演。

Description

一种基于重力数据的多层地壳结构反演方法
技术领域
本发明涉及地壳结构反演领域，具体涉及一种基于重力数据的多层地壳结构反演方法。
背景技术
地壳一般分为上、中、下三层（包括Pg层、P1层、P2层、P3层、Pm层，其中Pg层为沉积层的底界，P1、P2、P3层为地壳的内部分层，Pm为莫霍面）。在基于重力数据的地壳结构反演中，通常假设地壳密度均匀，反演地壳的底面，获得莫霍面埋深，讨论莫霍面起伏。其中常用的Parker法是一种频率域位场方法，将空间域位场数据变换到频率域进行数据处理，可以充分发挥位场数据的横向分辨率优势，使得空间域复杂的正演计算公式在频率域中大幅度简化。引入泰勒公式，通过一系列不同阶次的导数来逼近起伏界面，配合相应的反演算法，可以为界面起伏类的位场异常解释提供快速解决方案。除反演莫霍面外，Parker法还可以用于反演任意深度的密度界面起伏，只要输入相应的界面平均深度以及界面密度差参数。
但是在地壳结构的研究中，除了莫霍面的深度，还需要关注上、中、下地壳之间的分界面。多层界面反演中除了单层界面反演所需要考虑的内容外，还需要考虑首先对多个层面间不同分层面引起的重力异常分配进行识别。
发明内容
针对现有技术中的上述不足，本发明提供的一种基于重力数据的多层地壳结构反演方法根据不同深度密度异常引起的重力异常频谱差异，分别提取不同地壳分界面起伏引起的重力异常，据此实现地壳多层面几何形态反演。
为了达到上述发明目的，本发明采用的技术方案为：
提供一种基于重力数据的多层地壳结构反演方法，其包括以下步骤：
S1、获取布格重力异常数据，并对其进行多尺度小波分解，将重力异常分解为不同深度的信号；
S2、在不同深度的信号中寻找对应深度并进行层面反演，得到第一反演结果；
S3、基于第一反演结果建立三维地壳分层初始模型；
S4、调整三维地壳分层初始模型的分层的密度参数，获取对应的正演残差，通过残差值固定密度参数值；
S5、对固定了密度参数值的三维地壳分层初始模型进行正演，根据正演残差迭代调整各层深度，得到第二反演结果，即基于重力数据的多层地壳结构反演结果。
本发明的有益效果为：本发明基于布格重力异常数据，通过多尺度小波分解，将重力异常分解为不同深度的信号，寻找对应深度，进行层面反演，利用反演结果建立三维地壳分层初始模型；为初始模型的分层进行密度参数赋值，观察残差，调整密度参数，固定密度赋值；在固定密度参数赋值的基础上，对模型进行正演，根据正演残差迭代调整各层深度，实现地壳多层面几何形态反演。
附图说明
图1为本方法的流程示意图；
图2为实施例中Pg层埋深统计示意图；
图3为实施例中P1层埋深统计示意图；
图4为实施例中P2层埋深统计示意图；
图5为实施例中P3层埋深统计示意图；
图6为实施例中莫霍面埋深统计示意图；
图7为实施例中的重力异常小波多尺度分解结果；
图8为实施例中的重力异常小波多尺度分解结果；
图9为实施例中的重力异常小波多尺度分解结果；
图10为实施例中的重力异常小波多尺度分解结果；
图11为实施例中的重力异常小波多尺度分解结果；
图12为实施例中莫霍面深度对比示意图；
图13为实施例中中地壳底界深度对比示意图；
图14为实施例中上地壳底界深度对比示意图；
图15为实施例中三维地壳分层初始模型示意图；
图16为实施例中A4重力观测数据；
图17为实施例中莫霍面第3次迭代正演计算获得的模拟重力场；
图18为实施例中莫霍面第3次迭代模拟重力场对观测值的模拟残差；
图19为实施例中A3重力观测数据；
图20为实施例中中地壳底界第6次迭代后模型正演计算获得的模拟重力场；
图21为实施例中中地壳底界第6次迭代模拟重力场对观测值的模拟残差；
图22为实施例中A2重力观测数据；
图23为实施例中上地壳底界第25次迭代后模型正演计算获得的模拟重力场；
图24为实施例中上地壳底界第25次迭代模拟重力场对观测值的模拟残差；
图25为实施例中模型残差标准差变化曲线；
图26为实施例中25次迭代后的地壳分层界面三维形态。
具体实施方式
下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。
如图1所示，该基于重力数据的多层地壳结构反演方法包括以下步骤：
S1、获取布格重力异常数据，并对其进行多尺度小波分解，将重力异常分解为不同深度的信号；
S2、在不同深度的信号中寻找对应深度并进行层面反演，得到第一反演结果；
S3、基于第一反演结果建立三维地壳分层初始模型；
S4、调整三维地壳分层初始模型的分层的密度参数，获取对应的正演残差，通过残差值固定密度参数值；
S5、对固定了密度参数值的三维地壳分层初始模型进行正演，根据正演残差迭代调整各层深度，得到第二反演结果，即基于重力数据的多层地壳结构反演结果。
步骤S1的具体方法包括以下子步骤：
S1-1、获取布格重力异常数据，并根据公式：
进行多尺度小波分解；其中为布格重力异常，/>是小波变换第i阶逼近部分，/>是小波变换第/>阶高频细节部分，/>；/>表示区域、深部重力异常值；/>表示局部、浅部重力异常值；
S1-2、分别对小波变换第阶逼近部分和小波变换第/>阶高频细节部分进行场源深度估计，得到小波分解结果似场源深度估计值，并得到不同深度的信号。
步骤S2中的具体方法为：
对照地壳分层大致深度范围约束信息与小波分解结果似场源深度估计值选择界面反演依据并进行层面反演，得到第一反演结果：
将第阶小波分解的特征深度范围记为/>；其中/>为小波细节结果/>的似场源深度值，/>为小波逼近结果/>的似场源深度估计值；将莫霍面约束数据的深度范围记做，将中地壳底面约束数据的深度范围记做/>，将上地壳底面约束数据的深度范围记做/>；
使用完成地形、中间层改正的布格重力异常数据，选择特征深度范围与分布范围最接近的阶数，使得其特征深度范围/>最接近莫霍面的约束深度范围，以小波逼近结果/>为重力观测数据，以0海拔平面为上界面，使用Parker法下界面反演方法，得到整个地壳的底面，也即莫霍面；
选择特征深度范围与中地壳底界深度范围/>最相似的阶数/>，以/>为重力观测数据，以莫霍面为下界面，使用Parker法上界面反演方法，得到下地壳的顶界，也即中地壳底界；
选择特征深度范围与上地壳底界深度范围/>最相似的阶数，以/>为重力观测数据，以中地壳底界为下界面，使用Parker法上界面反演方法，得到中地壳的顶界，也即上地壳底界。
在具体实施过程中，使用完成地形、中间层改正的布格重力异常数据，选择特征深度范围与分布范围最接近的阶数（其特征深度范围为30km～63km，接近莫霍面的约束深度范围28km～46km，）以小波逼近结果/>为重力观测数据，以0海拔平面为上界面，使用Parker法下界面反演方法，得到整个地壳的底面，也即莫霍面；
选择特征深度范围最接近的阶数3（其特征深度范围为16km～37km，接近中地壳底界深度范围20km～33km），以为重力观测数据，以莫霍面为下界面，使用Parker法上界面反演方法，得到下地壳的顶界，也即中地壳底界；
选择特征深度范围最接近的阶数2（其特征深度范围为10km～20km，接近上地壳底界深度范围15km～24km），以为重力观测数据，以中地壳底界为下界面，使用Parker法上界面反演方法，得到中地壳的顶界，也即上地壳底界。
地壳分层大致深度范围约束信息包括：
莫霍面的深度在28km～46km之间；
下地壳与中地壳分界深度范围在20km～33km之间；
中地壳与上地壳分界深度范围在15km～24km之间。
步骤S3的具体方法包括以下子步骤：
S3-1、对第一反演结果中不同的界面进行插值，生成三角网；
S3-2、根据三角网的网格数据，计算界面上各网点的三维坐标、倾向信息和倾角信息；
S3-3、基于Geomodeller软件，根据建模范围、重力数据精度和模型分辨率信息，规划X轴方向上的剖面，生成控制剖面网；
S3-4、根据剖面位置，设置在剖面上的取样点，并获取取样点上的界面名称、三维坐标、倾向信息和倾角信息；
S3-5、根据剖面与取样点的三维坐标、倾向信息和倾角信息，计算取样点在剖面上的局部二维坐标（u，v），其中u为沿剖面走向坐标，v为沿剖面深度坐标；
S3-6、基于剖面与取样点在剖面上的局部二维坐标（u，v）进行建模，得到三维地壳分层初始模型。
步骤S4的具体方法包括以下子步骤：
S4-1、根据全球地壳模型crust1.0的分层进行初始密度参数值设置：设定莫霍面以下的上地幔顶部密度为2.9g/cm3，下地壳密度为2.6 g/cm3，中地壳密度为2.44 g/cm3，上地壳密度为2.21 g/cm3；
S4-2、基于当前密度值，对三维地壳分层初始模型进行正演，获取正演残差；
S4-3、判断当前的正演残差的标准差是否趋于平稳，若是则输出当前的密度参数值；否则分层调整密度参数值并返回步骤S4-2。
步骤S5的具体方法包括以下子步骤：
S5-1、对固定了密度参数值的三维地壳分层初始模型进行正演，分别以、/>和作为观测值，获取正演残差/>、/>和/>；
S5-2、判断当前正演残差的标准差是否趋于平稳，若是则输出当前三维地壳分层模型，得到第二反演结果；否则进入步骤S5-3；
S5-3、根据公式：
获取上地壳底界深度调整值、中地壳底界深度调整值/>和莫霍面深度调整值/>；其中/>为调整步长因子；
S5-4、根据当前的深度调整值计算各层的倾向信息和倾角信息；
S5-5、基于当前的倾向信息和倾角信息重新计算剖面上界面中采样点数据，建立新的三维地壳分层模型，维持密度参数值不变，返回步骤S5-1。
在本发明的一个实施例中，在地质体形态的反演中，初始模型的合理构建是解决多解性的重要途径。一般可用于构建初始模型的参考资料包括地壳尺度的深地震反射剖面、全球地壳分层模型Crust等。地壳内部的界面往往比较复杂，不同的研究对地震剖面可能有不同的解释，不同资料中对同一地壳内部界面的深度认识也有差异，由于勘探手段的局限性，地下的情况难以精确感知，因此无法判断多种资料中哪一种最可靠。因此需要尽可能完整的收集不同学者、不同手段获得的地壳界面深度信息，再对这些数据进行统计，制作箱图，得出地壳不同分层的大致深度范围。
一般深地震反射剖面相关文献会同时提供反射剖面线的平面图、解译的剖面图。根据反射剖面线的平面图，数字化剖面线。结合剖面线位置、剖面图上的坐标值等，对剖面图进行二维-三维投影，将剖面图回复到实际的三维空间中，建立该剖面上坐标系统与实际三维坐标系统之间的坐标转换关系。
在深地震反射剖面的解译图上，一般会识别出Pg、P1、P2、P3、Pm等层面。其中的Pg表示沉积层底部，一般在5km以内，P1-P3为地壳内部分层，除了上、中、下地壳外，还将下地壳分为上下两部分。Pm为莫霍面，是地壳和地幔之间的分层。数字化后的分层线上的节点根据二-三维坐标转换关系，转换为三维空间中的分界点，具有x、y、z坐标值。
深地震反射剖面往往长度较长，或仅有部分落在研究区内。选择落在研究区内的部分点。区内地壳分层面存在高低起伏，因此相应的地壳分层的深度数据往往不是单个值，而是一个范围。不同剖面上，或相同剖面的不同解译，也往往存在不一致。因此要对研究区内数字化之后的层面点的深度进行统计，制作深度分布箱图，对比不同的结果，对已有的不同分界面的深度范围与认识差异有一个大致的认识。如在某地区，选择了区内前人发表的三个剖面，分别为孟连-马龙剖面（简写为MM），玉溪-临沧剖面（简写为YL），思茅-中甸剖面（简写为SZ）。三个剖面上都对地壳的分层进行了描绘，包括Pg层，P1层，P2层，P3层、Pm层。其中Pg层位沉积层的底界，P1、P2、P3层为地壳的内部分层，Pm为莫霍面。通过上述处理方式，对分层上的点的深度进行统计，制作箱图，不同剖面上的分层点深度对比如图2、图3、图4、图5和图6所示。
crust 1.0是一个常用的全球地壳模型，模型精度为1°×1°；网格点定义在单元的中心，模型分为8层（也可以说是9层）；对于每个1°×1°的单元，给出每层的深度（crust1.bnds）、Vp（crust1.vp）、Vs（crust1.vs）、ρ（crust1.rho）。从上到下每层对应水层、冰层、上沉积层、中沉积层、下沉积层、上地壳、中地壳、下地壳。
根据模型中属性点的x、y值，筛选研究区内的点，并对区内点的界面深度、波速、密度进行统计，制作箱图，作为初始界面深度、初始密度属性等的参考。
在本方法的具体实施过程中，地表测量的重力异常，是地下不同深度范围内所有重力异常场源的叠加结果。为了从中分辨出不同深度的异常，通常使用小波多尺度分解的方法，根据异常空间频率的区别，将异常分解成表达不同深度场源的重力异常，其中的高频部分代表浅部、局部异常，低频部分代表深部、区域异常。某区域的小波分解结果似场源深度估计值如表1所示。、/>、/>、/>和/>的重力异常小波多尺度分解结果分别如图7、图8、图9、图10和图11所示。
表1：小波分解结果似场源深度
在本实施例中，D4似场源深度约30km，A4似场源深度为63km，假设存在一界面将地壳分为上下两部分，其上部分产生的异常为D4+D3+D2+D1，由D4计算的似场源深度结果会略浅于该界面深度。其下部分产生的异常为A4，由A4计算的似场源深度结果会略深于该界面。也即将地壳分成A4和D4+ D3+ D2+ D1的界面其深度在30-63km之间。这与约束数据统计显示莫霍面的深度在28～46km之间莫霍面的深度相似，因此基于A4结果进行莫霍面的反演较为合理。
下地壳与中地壳分界，其深度范围在20～33km之间，而D3深度为16km，A3深度为37km，那么以莫霍面为底界，基于D4进行Parker法上界面反演，获得下地壳顶，也即中地壳底，其模拟的重力异常为A3，A3G=A4G +D4G。
中地壳与上地壳分界，其深度范围在15～24km之间，而D2深度为10km，A2深度为20km，那么以P2为底界，基于D3进行Parker法上界面反演，获得中地壳顶界，也即上地壳底界，其模拟的重力异常为A2，A2G =A3G +D3G。
在使用Parker法进行密度界面反演时，平均深度与界面密度两个参数初始参数选择约束数据的中值。需要根据约束数据输入平均深度与界面密度两个参数，需要根据反演结果与点约束的对比，来调整阶数、迭代次数、正则化因子等参数。将反演结果与单点深度约束进行对比，计算残差，并对残差分布进行分析。若残差均值偏离0较大，说明平均深度的选取存在偏差，对平均深度参数进行调整，调整范围在约束数据最大最小值范围内。若反演结果深度起伏过大，对界面密度进行调整，界面密度值越大，界面起伏越小，最终使得反演结果的最大最小值与约束数据最大最小值相当。若残差标准差过大，或层面起伏形态与输入重力异常之间的空间特征差异过大，则调整阶数、迭代次数、正则化因子等参数。阶数越大、迭代次数越大、正则化因子越大，则界面起伏增大，形态复杂化。
所有的统计针对研究区内的点进行，图12中crustPm代表与Crust1.0模型之间的莫霍面深度差异，MMPm代表与孟连-马龙剖面间的莫霍面深度差异，SZPm代表与思茅-中甸剖面的莫霍面深度差异，YLPm代表与玉溪-临沧剖面见的莫霍面深度差异；图13中crustP3代表与Crust1.0模型之间的中地壳底界深度差异，MMP3代表与孟连-马龙剖面间的中地壳底界深度差异，SZP3代表与思茅-中甸剖面的中地壳底界深度差异，YLP3代表与玉溪-临沧剖面间的中地壳底界深度差异；图14中的crustP2代表与Crust1.0模型之间的上地壳底界深度差异，MMP2代表与孟连-马龙剖面间的上地壳底界深度差异，SZP2代表与思茅-中甸剖面的上地壳底界深度差异，YLP2代表与玉溪-临沧剖面间的上地壳底界深度差异。
从图12至图14可以看出，主要以Crust1.0给出的深度值作为参考进行反演，对比结果显示反演结果与Crust1.0的逐点偏差的中值接近0，多数偏差在5km以内，最大不超过8km。
在本发明的一个实施例中，重力数据比例尺为10万，换算成分辨率为10km，因此模型分辨率不应大于10km。在两侧给出微量（不大于1km）边界的情况下，从下边界开始，每隔10km规划一条剖面，直至超出数据覆盖范围。在两侧给出微量（不大于1km）边界的情况下，从左边界开始，每隔10km规划一个取样点，直至超出数据覆盖范围。随着控制剖面网和采样点的增加，插值后的模型逐步逼近反演的界面，建模结果如图15所示，椎体为采样点，浅色点为界面反演结果，灰色面为界面三维模型。
在本发明的一个实施例中，参考crust1.0的分层密度设置，设定莫霍面以下的上地幔顶部密度为2.9，下地壳密度为2.6/>，中地壳密度为2.44/>，上地壳密度为2.21/>。以A2为整体的拟合目标，直接用这一套参数，正演残差较大，范围在-20～130之间。模拟场能够刻画出这一区域的重力的区域特性，细节刻画不足。
通过分层密度属性的优化，残差可以减少到-35～44之间，分层密度参数从上到下分别为2.31、2.54/>、2.7/>、2.8/>。密度的绝对值大小意义不大，其与模型的建模深度范围关系密切，该模型的建模深度为60km以上。对密度值的优化，可减小残差的数值，但基本不改变残差的形态。残差的进一步缩小需要通过改变界面形态来实现。
在本发明的具体实施过程中，A4重力观测数据如图16所示，莫霍面第3次迭代正演计算获得的模拟重力场如图17所示，莫霍面第3次迭代模拟重力场对观测值的模拟残差如图18所示，A3重力观测数据如图19所示，中地壳底界第6次迭代后模型正演计算获得的模拟重力场如图20所示，中地壳底界第6次迭代模拟重力场对观测值的模拟残差如图21所示，A2重力观测数据如图22所示，上地壳底界第25次迭代后模型正演计算获得的模拟重力场如图23所示，上地壳底界第25次迭代模拟重力场对观测值的模拟残差如图24所示，模型正演残差的标准差变化曲线如图25所示，25次迭代后的地壳分层界面三维形态如图26所示。可以看出，采用本方法调整各层深度，残差收敛较好，对重力观测值的拟合更好，模型对A2拟合的残差标准差约为1mGal，最终得到的反演结果细节饱满。
综上所述，本发明基于布格重力异常数据，通过多尺度小波分解，将重力异常分解为不同深度的信号，寻找对应深度，进行层面反演，利用反演结果建立三维地壳分层初始模型；为初始模型的分层进行密度参数赋值，观察残差，调整密度参数，固定密度赋值；在固定密度参数赋值的基础上，对模型进行正演，根据正演残差迭代调整各层深度，实现地壳多层面几何形态反演。

Claims (3)

1.一种基于重力数据的多层地壳结构反演方法，其特征在于，包括以下步骤：
S1、获取布格重力异常数据，并对其进行多尺度小波分解，将重力异常分解为不同深度的信号；
S2、在不同深度的信号中寻找对应深度并进行层面反演，得到第一反演结果；
S3、基于第一反演结果建立三维地壳分层初始模型；
S4、调整三维地壳分层初始模型的分层的密度参数，获取对应的正演残差，通过残差值固定密度参数值；
S5、对固定了密度参数值的三维地壳分层初始模型进行正演，根据正演残差迭代调整各层深度，得到第二反演结果，即基于重力数据的多层地壳结构反演结果；
步骤S1的具体方法包括以下子步骤：
S1-1、获取布格重力异常数据，并根据公式：
1阶小波分解
2阶小波分解
3阶小波分解
4阶小波分解；
进行多尺度小波分解；其中为布格重力异常，/>是小波变换第i阶逼近部分，是小波变换第i阶高频细节部分，i=1,2,3,4；/>表示区域、深部重力异常值；/>表示局部、浅部重力异常值；
S1-2、分别对小波变换第i阶逼近部分和小波变换第i阶高频细节部分进行场源深度估计，得到小波分解结果似场源深度估计值，并得到不同深度的信号；
步骤S2中的具体方法为：
对照地壳分层大致深度范围约束信息与小波分解结果似场源深度估计值选择界面反演依据并进行层面反演，得到第一反演结果：
将第阶小波分解的特征深度范围记为/>；其中/>为小波细节结果/>的似场源深度值，/>为小波逼近结果/>的似场源深度估计值；将莫霍面约束数据的深度范围记做，将中地壳底面约束数据的深度范围记做/>，将上地壳底面约束数据的深度范围记做/>；
使用完成地形、中间层改正的布格重力异常数据，选择特征深度范围与分布范围最接近的阶数，使得其特征深度范围/>最接近莫霍面的约束深度范围/>，以小波逼近结果/>为重力观测数据，以0海拔平面为上界面，使用Parker法下界面反演方法，得到整个地壳的底面，也即莫霍面；
选择特征深度范围与中地壳底界深度范围/>最相似的阶数/>，以为重力观测数据，以莫霍面为下界面，使用Parker法上界面反演方法，得到下地壳的顶界，也即中地壳底界；
选择特征深度范围与上地壳底界深度范围/>最相似的阶数/>，以为重力观测数据，以中地壳底界为下界面，使用Parker法上界面反演方法，得到中地壳的顶界，也即上地壳底界；
步骤S3的具体方法包括以下子步骤：
S3-1、对第一反演结果中不同的界面进行插值，生成三角网；
S3-2、根据三角网的网格数据，计算界面上各网点的三维坐标、倾向信息和倾角信息；
S3-3、基于Geomodeller软件，根据建模范围、重力数据精度和模型分辨率信息，规划X轴方向上的剖面，生成控制剖面网；
S3-4、根据剖面位置，设置在剖面上的取样点，并获取取样点上的界面名称、三维坐标、倾向信息和倾角信息；
S3-5、根据剖面与取样点的三维坐标、倾向信息和倾角信息，计算取样点在剖面上的局部二维坐标（u，v），其中u为沿剖面走向坐标，v为沿剖面深度坐标；
S3-6、基于剖面与取样点在剖面上的局部二维坐标（u，v）进行建模，得到三维地壳分层初始模型；
步骤S5的具体方法包括以下子步骤：
S5-1、对固定了密度参数值的三维地壳分层初始模型进行正演，分别以、/>和/>作为观测值，获取正演残差/>、/>和/>；
S5-2、判断当前正演残差的标准差是否趋于平稳，若是则输出当前三维地壳分层模型，得到第二反演结果；否则进入步骤S5-3；
S5-3、根据公式：
；
；
；
获取上地壳底界深度调整值、中地壳底界深度调整值/>和莫霍面深度调整值；其中/>为调整步长因子；
S5-4、根据当前的深度调整值计算各层的倾向信息和倾角信息；
S5-5、基于当前的倾向信息和倾角信息重新计算剖面上界面中采样点数据，建立新的三维地壳分层模型，维持密度参数值不变，返回步骤S5-1。
2.根据权利要求1所述的基于重力数据的多层地壳结构反演方法，其特征在于，地壳分层大致深度范围约束信息包括：
莫霍面的深度在28km～46km之间；
下地壳与中地壳分界深度范围在20km～33km之间；
中地壳与上地壳分界深度范围在15km～24km之间。
3.根据权利要求1所述的基于重力数据的多层地壳结构反演方法，其特征在于，步骤S4的具体方法包括以下子步骤：
S4-1、根据全球地壳模型crust1.0的分层进行初始密度参数值设置：设定莫霍面以下的上地幔顶部密度为2.9g/cm3，下地壳密度为2.6 g/cm3，中地壳密度为2.44 g/cm3，上地壳密度为2.21 g/cm3；
S4-2、基于当前密度值，对三维地壳分层初始模型进行正演，获取正演残差；
S4-3、判断当前的正演残差的标准差是否趋于平稳，若是则输出当前的密度参数值；否则分层调整密度参数值并返回步骤S4-2。
CN202410039732.5A 2024-01-11 2024-01-11 一种基于重力数据的多层地壳结构反演方法 Active CN117555025B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202410039732.5A CN117555025B (zh) 2024-01-11 2024-01-11 一种基于重力数据的多层地壳结构反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202410039732.5A CN117555025B (zh) 2024-01-11 2024-01-11 一种基于重力数据的多层地壳结构反演方法

Publications (2)

Publication Number Publication Date
CN117555025A CN117555025A (zh) 2024-02-13
CN117555025B true CN117555025B (zh) 2024-04-02

Family

ID=89820876

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202410039732.5A Active CN117555025B (zh) 2024-01-11 2024-01-11 一种基于重力数据的多层地壳结构反演方法

Country Status (1)

Country Link
CN (1) CN117555025B (zh)

Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104459795A (zh) * 2014-12-08 2015-03-25 中国科学院南海海洋研究所 一种深度变密度的地壳伸展系数热校正重力异常反演方法
JP2017068089A (ja) * 2015-09-30 2017-04-06 キヤノン株式会社 結像光学系及びそれを備える画像読取装置
EP3221724A1 (en) * 2012-05-22 2017-09-27 Nxt Energy Solutions, Inc. Gravity transducer
CN110221344A (zh) * 2019-06-17 2019-09-10 中国地质大学（北京） 一种地壳三维密度结构的地震全波形与重力联合反演方法
CN110244352A (zh) * 2019-06-11 2019-09-17 西安石油大学 一种基于变密度的地壳厚度重力反演方法
CN111650641A (zh) * 2020-06-05 2020-09-11 河南工业大学 一种地壳三维结构模型融合方法及装置
CN112883454A (zh) * 2020-12-29 2021-06-01 煤炭科学研究总院 一种快速评估地震作用下非饱和边坡永久变形的方法和装置
CN113985490A (zh) * 2021-09-22 2022-01-28 中国人民解放军战略支援部队信息工程大学 利用地形和地壳密度数据进行地表重力仿真的方法及装置
CN115437027A (zh) * 2022-09-21 2022-12-06 中国地质调查局西安地质调查中心(西北地质科技创新中心) 利用地质信息变密度正演计算布格重力异常的方法及装置
CN116187168A (zh) * 2022-12-30 2023-05-30 中国航天科技创新研究院 基于神经网络-重力信息小波分解提高水深反演精度方法
WO2023098441A1 (zh) * 2022-08-09 2023-06-08 中国科学院南海海洋研究所 基于地层记录沉降反演被动陆缘地壳结构的方法及装置

Patent Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP3221724A1 (en) * 2012-05-22 2017-09-27 Nxt Energy Solutions, Inc. Gravity transducer
CN104459795A (zh) * 2014-12-08 2015-03-25 中国科学院南海海洋研究所 一种深度变密度的地壳伸展系数热校正重力异常反演方法
JP2017068089A (ja) * 2015-09-30 2017-04-06 キヤノン株式会社 結像光学系及びそれを備える画像読取装置
CN110244352A (zh) * 2019-06-11 2019-09-17 西安石油大学 一种基于变密度的地壳厚度重力反演方法
CN110221344A (zh) * 2019-06-17 2019-09-10 中国地质大学（北京） 一种地壳三维密度结构的地震全波形与重力联合反演方法
CN111650641A (zh) * 2020-06-05 2020-09-11 河南工业大学 一种地壳三维结构模型融合方法及装置
CN112883454A (zh) * 2020-12-29 2021-06-01 煤炭科学研究总院 一种快速评估地震作用下非饱和边坡永久变形的方法和装置
CN113985490A (zh) * 2021-09-22 2022-01-28 中国人民解放军战略支援部队信息工程大学 利用地形和地壳密度数据进行地表重力仿真的方法及装置
WO2023098441A1 (zh) * 2022-08-09 2023-06-08 中国科学院南海海洋研究所 基于地层记录沉降反演被动陆缘地壳结构的方法及装置
CN115437027A (zh) * 2022-09-21 2022-12-06 中国地质调查局西安地质调查中心(西北地质科技创新中心) 利用地质信息变密度正演计算布格重力异常的方法及装置
CN116187168A (zh) * 2022-12-30 2023-05-30 中国航天科技创新研究院 基于神经网络-重力信息小波分解提高水深反演精度方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
利用小波变换提取川滇地区流动重力异常特征――以2014年鲁甸M_S6.5、景谷M_S6.6地震为例;张双喜;陈兆辉;王青华;刘金钊;张品;;大地测量与地球动力学;20200115(01);第91-97页 *

Also Published As

Publication number Publication date
CN117555025A (zh) 2024-02-13

Similar Documents

Publication Publication Date Title
AU2016305571B2 (en) System and method for gravity and/or gravity gradient terrain corrections
US8437960B2 (en) Gravity survey data processing
AU2006276744B2 (en) Method for tomographic inversion by matrix transformation
CN112363236B (zh) 一种基于pde的重力场数据等效源延拓与数据类型转换方法
CN112147709B (zh) 一种基于部分光滑约束的重力梯度数据三维反演方法
CN111340723B (zh) 一种地形自适应的机载LiDAR点云正则化薄板样条插值滤波方法
CN111305834B (zh) 基于多探测模式电阻率测井的三维反演初始模型构建方法
WO2016036411A1 (en) Sediment transport simulation with parameterized templates for depth profiling
CN114943178A (zh) 一种三维地质模型建模方法、装置及计算机设备
CN111221035A (zh) 一种地震反射波斜率和重力异常数据联合反演方法
CN114896564A (zh) 采用自适应泰森多边形参数化的瞬变电磁二维贝叶斯反演方法
Mastellone et al. Volume Continuation of potential fields from the minimum-length solution: An optimal tool for continuation through general surfaces
CN117555025B (zh) 一种基于重力数据的多层地壳结构反演方法
RU2431873C2 (ru) Обработка данных гравиметрической разведки
CN113109793A (zh) 自适应分辨率水深曲面滤波方法及装置
CN110989021B (zh) 水深反演方法、装置和计算机可读存储介质
Xing et al. Bathymetry inversion using the modified gravity-geologic method: application of the rectangular prism model and Tikhonov regularization
CN116187168A (zh) 基于神经网络-重力信息小波分解提高水深反演精度方法
Soycan et al. Digital elevation model production from scanned topographic contour maps via thin plate spline interpolation
Chen et al. The improved Kriging interpolation algorithm for local underwater terrain based on fractal compensation
Castiello et al. Enhanced methods for interpreting microgravity anomalies in urban areas
EP4073554B1 (en) Layering for geomodeling
US20230057978A1 (en) Painting For Geomodeling
Liu et al. An Extension of Multiquadric Method Based on Trend Analysis for Surface Construction
CN114779357A (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