CN112147709B - 一种基于部分光滑约束的重力梯度数据三维反演方法 - Google Patents
一种基于部分光滑约束的重力梯度数据三维反演方法 Download PDFInfo
- Publication number
- CN112147709B CN112147709B CN202010766212.6A CN202010766212A CN112147709B CN 112147709 B CN112147709 B CN 112147709B CN 202010766212 A CN202010766212 A CN 202010766212A CN 112147709 B CN112147709 B CN 112147709B
- Authority
- CN
- China
- Prior art keywords
- matrix
- model
- gravity gradient
- roughness
- 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.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V7/00—Measuring gravitational fields or waves; Gravimetric prospecting or detecting
Landscapes
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明涉及一种基于部分光滑约束的重力梯度数据三维反演方法,首先计算长方体单元引力位在三维空间x,y,z三个方向上的二阶偏导数得到长方体元的正演公式,给定模型后根据正演公式正演计算得到观测数据;再依据一定的先验信息对粗糙度约束矩阵进行人为赋值,得到粗糙度矩阵,并采用积分灵敏度的加权方式克服“趋肤效应”;最后写出含有积分灵敏度矩阵,粗糙度矩阵,正则化参数,观测数据等参数的目标方程,并采用共轭梯度法解目标方程,得到反演结果。该技术适用于在有一定地质或地球物理先验条件下的三维重力和重力梯度数据反演。
Description
技术领域
本发明涉及地球物理勘探领域中三维重力梯度数据反演问题,具体说,涉及一种基于部分光滑约束的重力梯度数据三维反演方法。
背景技术
如何快速准确的确定地质***置和物性是找矿地质勘探中最重要的问题之一,而重力或重力梯度的物性反演技术是解决该问题的关键技术环节。在实际工作中,因为 采集的有效数据有限,测量精度受各种因素干扰大,特别是场源的等效性使得三维物 性反演存在严重的多解性,故反演问题一般是欠定问题。解决这种问题目前最常用的 办法就是加入地质和其他地球物理方法提供的先验信息或者合理适当的约束条件来处 理。Li和Oldenburg(1998)提出的粗糙度矩阵给出了一种反演的光滑模型,并通过 增加深度加权矩阵压制了“趋肤效应”,Li(2000)进一步将光滑反演方法推广到重力 梯度全张量数据的反演中,Martinez和Li(2013)采用该方法对航空重力梯度张量数 据的不同组合进行反演,证明了全张量数据联合反演结果能够更好的反演地质体结构。 Silva等(2001)讨论了现有物性反演方法中的约束条件,认为应根据实际的地质问题 选择约束条件。光滑反演一般采取的是全局光滑,得到的反演结果边界模糊且范围扩 大,模型分辨率较低。
发明内容
本发明所要解决的问题是提高反演结果的分辨率,进而较真实反映地下地质体的位置和物性。
本发明提供一种基于部分光滑约束的重力梯度数据三维反演方法,包括如下步骤:
1)分别计算引力位在三维空间x,y,z三个方向上的二阶偏导数,得到重力梯 度正演公式,再利用上述公式,计算给定模型正演重力梯度值作为反演的观测数据。
2)依据先验信息对粗糙度约束矩阵进行人为赋值,得到粗糙度矩阵;
3)采用积分灵敏度的加权方式形成积分灵敏度矩阵。
4)基于步骤1)中观测数据,步骤2)中的粗糙度矩阵,步骤3)中的积分灵敏 度矩阵,选取正则化参数,得到目标方程,
5)采用共轭梯度法解目标方程,得到反演结果。
所述步骤1)具体为:以(ξ,η,ζ)为长方体内体积微元的坐标,长方体单元坐标值最小和最大的两个点的坐标分别为(ξ1,η1,ζ1),(ξ2,η2,ζ2),点(x,y,z)为观测点,舒晴 等人推导得到的无解析奇点的重力梯度全张量正演计算公式如下所示:
其中,G为万有引力常量,r=[(x-ξ)2+(y-η)2+(z-ζ)2]1/2,ρ为长方体的剩余 密度,s=(-1)i+j+k。
所述步骤1)还包括:建立一个三维地质理论模型,将地下空间剖分为J行K列 L层整齐排列的长方体单元,对于不同的正演模型,各长方体单元除密度值外完全相 同,以靠近地面那一层每个单元的顶部中心位置为观测点,根据正演公式计算观测点 重力梯度值。
所述步骤2)中,因为重力及其梯度数据对于异常体深度信息不敏感,故一般需 要其它地质或地球物理方法提供异常体的深度信息,采用深度信息作为先验信息进行 约束。同时,异常体倾向也可以作为先验信息添加到粗糙度矩阵中。
上述先验信息可以由其它地质或地球物理信息提供,如异常体埋深,异常体倾向等,先验信息越多则反演结果越准确。
所述步骤2)中的粗糙度矩阵为:
式中M为光滑约束的条件数,N为地下模型单元的个数。
所述步骤3)具体为:观测数据d与模型m变化的关系可表示为
δdi=Fikδmk
其中,Fik为观测数据关于参数的灵敏度(导数)。数据灵敏度对模型mk的积分 可表示为
故积分灵敏度矩阵S为对角矩阵;
S=diag(FTF)1/2
其中,F为观测数据关于参数的灵敏度,FT为F的转置,diag表示取对角矩阵。
所述步骤4)中的目标方程为:φ(m)=(d-Gm)T(d-Gm)+αmTWm TWmm;式中,α为 正则化参数,Wm=RS,S为积分灵敏度对角矩阵,R为粗糙度矩阵,m为要求解的模 型单元密度值。
所述步骤5)中,目标函数可以写成以下增广矩阵的形式:
采用共轭梯度法解上述目标函数,得到矩阵[ρ]的解。
本发明的有益效果在于:相对于现有技术中针对当进行全局光滑约束时,反演结果模型边界模糊且范围扩大,模型分辨率较低,本发明首次提出了根据已知先验信息 控制光滑范围的部分光滑约束矩阵,将部分光滑约束矩阵加入到重力梯度异常数据三 维反演中。与现有技术相比,能够显著提高反演结果的分辨率,提高了模型反演结果 的准确性。
附图说明
图1三维重力梯度反演方法流程图;
图2模型一模型示意图;
图3模型一切片图
图4模型一全局光滑反演结果图;
图5模型一部分光滑反演结果图;
图6模型二切片图;
图7模型二部分光滑反演结果图。
具体实施方式
以下结合附图对本发明的原理和特征进行描述,所举实例只用于解释本发明,并非用于限定本发明的范围。
本发明涉及一种基于部分光滑约束的重力梯度数据三维反演方法,包括以下步骤:
a、建立笛卡尔直角坐标系,具体为:以(ξ,η,ζ)为长方体内体积微元的坐标,长 方体单元坐标值最小和最大的两个点的坐标分别为(ξ1,η1,ζ1),(ξ2,η2,ζ2),以点(x,y,z) 为观测点,则一个剩余密度为ρ长方体的重力梯度正演公式如下:
其中,G为万有引力常量,r=[(x-ξ)2+(y-η)2+(z-ζ)2]12,s=(-1)i+j+k。
b、输入观测数据d,这个观测数据由理论模型给出。将地下空间剖分为一系列21×21×10 的长方体单位元,其中每个单位元的大小为10m×10m×20m,观测点为最上层每个单位 元的顶部中心位置。模型一是一个地下存在两个异常体的三维地质模型。地质***置如图2所示。其中左侧地质体剩余密度为0.8g/cm3,右侧地质体剩余密度为1g/cm3, 做三维切片图如图3所示,正演其重力梯度数据Uzz,得到观测数据。模型二为台阶模 型,剩余密度为0.8g/cm3,模型切片图如图6所示,正演其UZZ作为观测数据。
c、粗糙度矩阵按照以下示例的形式给出。
式中M为光滑约束的条件数,N为地下模型单元的个数。如上述矩阵的第一行 前两个元素分别为-1和1,该行其余元素全部为零,则表示第一个模型单元和第二 个模型单元之间进行了光滑约束,每多一条约束,则矩阵增加一行,最终形成粗糙 度矩阵。粗糙度矩阵的具体形式由其他地质地球物理信息给出,如在本示例模型一 中已知该异常体的埋深,模型二中除了埋深还把异常体倾向作为先验信息。
d、由于各个剖分单元与观测网格点之间几何关系,导致了距离观测点越远的 剖分单元的重力梯度张量导数越小,如果不对模型进行加权约束就会使反演结果出 现“趋肤效应”,即异常都集中在地表很难反演出真实构造。因此,本文引入积分灵 敏度矩阵S,对模型进行加权,矩阵S由于具有x,y,z三方向加权效果,使得反演 具有三方向的分辨能力,因此,反演算法的旁侧反演能力得到大大提高。
观测数据与模型变化的关系可表示为
δdi=Fikδmk
其中,Fik为观测数据关于参数的灵敏度(导数)。数据灵敏度对模型mk的积分 可表示为
故积分灵敏度矩阵S为对角矩阵;
S=diag(FTF)1/2
其中,F为观测数据关于参数的灵敏度,FT为F的转置,diag表示取对角矩阵。
e、写出含有约束矩阵Wm,正则化参数,观测数据等参数的目标方程,即 φ(m)=(d-Gm)T(d-Gm)+αmTWm TWmm。式中,α为正则化参数,Wm=RS,S为积分灵 敏度对角矩阵,R为粗糙度矩阵。本文正则化参数选取的方法是一种易于编程的方法, 即选用α的初值为一个较大的常数值,之后每次迭代过程中α总为之前的一半,从而 实现需要光滑约束的模型空间尽可能光滑,且数据拟合误差在可接受范围内。
f、则目标函数可以写成以下增广矩阵的形式:
采用共轭梯度法解上述目标函数,得到矩阵[ρ]的解。
按照上述过程模型一反演结果如图4、图5所示,全局光滑反演结果如图4所示, 部分光滑时,采用左半侧空间上五层仅z方向光滑,右侧空间上四层仅z方向光滑, 使该空间范围内单元体之间x,y方向无约束,其余空间位置仍然采用x,y,z三方向的全 局光滑,反演结果如图6所示。其中红色线框表示地质体实际位置。
对比图4和图5,两种光滑方式均可以确定地质体的水平位置,但部分光滑反演 得到的地质体边界位置更加清楚,这是因为在全局光滑中地质体与围岩之间存在光滑 约束,使地质体边界模糊。
对比图4和图5,部分光滑相比于全局光滑有明显的优势,在图4中右侧异常体 纵向位置不明显,而图5中右侧地质体纵向位置可以准确的的反演得到,同时对于左 侧的异常体,图5边界更加清晰且反演得到的剩余密度值更接近于所给定模型的剩余 密度值。
对于模型二,采用部分光滑约束为上五层没有x,y方向的光滑,若已知地下模型可能为倾斜构造及构造方向,则添加右倾对角单元之间的光滑关系,下五层仍然采用 x,y,z三方向的全局光滑,反演结果如图7所示。可以看出,在浅部可以很好的得到地 质体的水平位置及纵向位置,在较深部区域水平位置略有偏差,则部分光滑反演基本 可以得到台阶模型的轮廓位置。
以上所述仅为本发明的较佳实施例,并不用以限制本发明,凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。
Claims (3)
1.一种基于部分光滑约束的重力梯度数据三维反演方法,其特征在于,包括如下步骤:
1)分别计算引力位在三维空间x,y,z三个方向上的二阶偏导数,得到重力梯度正演公式,再利用上述公式,计算给定模型正演重力梯度值作为反演的观测数据;
2)依据先验信息对粗糙度约束矩阵进行人为赋值,得到粗糙度矩阵;
3)采用积分灵敏度的加权方式形成积分灵敏度矩阵;
4)基于步骤1)中观测数据,步骤2)中的粗糙度矩阵,步骤3)中的积分灵敏度矩阵,选取正则化参数,得到目标方程,
5)采用共轭梯度法解目标方程,得到反演结果;
所述步骤2)中的先验信息为通过地质或地球物理方法提供异常体的深度信息或者异常体倾向;
所述步骤2)中的粗糙度矩阵为:
式中M为光滑约束的条件数,N为地下模型单元的个数,上述矩阵的第一行前两个元素分别为-1和1,该行其余元素全部为零,则表示第一个模型单元和第二个模型单元之间进行了光滑约束,每多一条约束,则矩阵增加一行,最终形成粗糙度矩阵;
所述步骤4)中的目标方程为:φ(m)=(d-Gm)T(d-Gm)+αmTWm TWmm;式中,α为正则化参数,Wm=RS,S为积分灵敏度对角矩阵,R为粗糙度矩阵,d为观测数据,m为模型,G为万有引力常量;
所述步骤5)中,目标方程为以下增广矩阵的形式:
采用共轭梯度法解上述目标方程,得到矩阵[ρ]的解,ρ为剩余密度。
3.根据权利要求2所述的基于部分光滑约束的重力梯度数据三维反演方法,其特征在于,所述步骤1)还包括:建立一个三维地质理论模型,将地下空间剖分为J行K列L层整齐排列的长方体单元,对于不同的正演模型,各长方体单元除密度值外完全相同,以靠近地面那一层每个单元的顶部中心位置为观测点,根据正演公式计算观测点重力梯度值。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010766212.6A CN112147709B (zh) | 2020-08-03 | 2020-08-03 | 一种基于部分光滑约束的重力梯度数据三维反演方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010766212.6A CN112147709B (zh) | 2020-08-03 | 2020-08-03 | 一种基于部分光滑约束的重力梯度数据三维反演方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN112147709A CN112147709A (zh) | 2020-12-29 |
CN112147709B true CN112147709B (zh) | 2022-07-29 |
Family
ID=73887815
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202010766212.6A Active CN112147709B (zh) | 2020-08-03 | 2020-08-03 | 一种基于部分光滑约束的重力梯度数据三维反演方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN112147709B (zh) |
Families Citing this family (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN115407423B (zh) * | 2021-05-27 | 2024-02-06 | 中国自然资源航空物探遥感中心 | 重、磁测数据三维反演方法及装置 |
CN113486503B (zh) * | 2021-06-24 | 2023-05-23 | 西南交通大学 | 一种重力及梯度异常正演方法 |
CN113504575B (zh) * | 2021-07-09 | 2022-05-03 | 吉林大学 | 基于权相交及多次交叉梯度约束的联合反演方法 |
CN113514900B (zh) * | 2021-07-12 | 2022-05-17 | 吉林大学 | 基于密度约束的球坐标系重力和重力梯度联合反演方法 |
CN113591030B (zh) * | 2021-08-17 | 2024-01-30 | 东北大学 | 基于多gpu的重力梯度数据灵敏度矩阵压缩及调用方法 |
CN114966878B (zh) * | 2022-04-12 | 2023-04-14 | 成都理工大学 | 一种基于混合范数和互相关约束的三维重力场反演方法 |
CN115220119B (zh) * | 2022-06-21 | 2023-02-24 | 广州海洋地质调查局 | 一种适用于大规模数据的重力反演方法 |
Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108957554A (zh) * | 2018-08-09 | 2018-12-07 | 北京探创资源科技有限公司 | 一种地球物理勘探中地震反演方法 |
Family Cites Families (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2017048445A1 (en) * | 2015-09-15 | 2017-03-23 | Exxonmobil Upstream Research Company | Accelerated occam inversion using model remapping and jacobian matrix decomposition |
CN110398782B (zh) * | 2019-07-17 | 2021-04-23 | 广州海洋地质调查局 | 一种重力数据和重力梯度数据联合正则化反演方法 |
CN111399074B (zh) * | 2020-04-28 | 2022-08-12 | 中国自然资源航空物探遥感中心 | 一种重力和重力梯度模量联合三维反演方法 |
-
2020
- 2020-08-03 CN CN202010766212.6A patent/CN112147709B/zh active Active
Patent Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108957554A (zh) * | 2018-08-09 | 2018-12-07 | 北京探创资源科技有限公司 | 一种地球物理勘探中地震反演方法 |
Also Published As
Publication number | Publication date |
---|---|
CN112147709A (zh) | 2020-12-29 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN112147709B (zh) | 一种基于部分光滑约束的重力梯度数据三维反演方法 | |
CN105549106B (zh) | 一种重力多界面反演方法 | |
CN112363236B (zh) | 一种基于pde的重力场数据等效源延拓与数据类型转换方法 | |
AU2015101990B4 (en) | Conditioning of object or event based reservoir models using local multiple-point statistics simulations | |
CN105334542B (zh) | 任意密度分布复杂地质体重力场快速、高精度正演方法 | |
GB2375635A (en) | a method for updating a fine geologic model | |
CN111337993A (zh) | 一种基于变密度变深度约束的重力密度界面反演方法 | |
CN113552625B (zh) | 一种用于常规陆域地震数据的多尺度全波形反演方法 | |
CN111856599B (zh) | 一种基于pde的磁测数据等效源化极与类型转换方法 | |
CN111221035B (zh) | 一种地震反射波斜率和重力异常数据联合反演方法 | |
CN111856598A (zh) | 一种磁测数据多层等效源上延拓与下延拓方法 | |
CN114943178A (zh) | 一种三维地质模型建模方法、装置及计算机设备 | |
CN109444973B (zh) | 一种球坐标系下重力正演加速方法 | |
CN107402409A (zh) | 一种三维不规则地层起伏界面重力正演方法 | |
CN112346139B (zh) | 一种重力数据多层等效源延拓与数据转换方法 | |
CN111859251B (zh) | 一种基于pde的磁测数据等效源上延拓与下延拓方法 | |
CN112596106A (zh) | 一种球坐标系下重震联合反演密度界面分布的方法 | |
Waheed et al. | A holistic approach to computing first-arrival traveltimes using neural networks | |
CN107748834A (zh) | 一种计算起伏观测面磁场的快速、高精度数值模拟方法 | |
CN114200541B (zh) | 一种基于余弦点积梯度约束的三维重磁联合反演方法 | |
CN112748471B (zh) | 一种非结构化等效源的重磁数据延拓与转换方法 | |
CN114740540A (zh) | 一种基于方向约束的洋中脊区磁异常图构建方法及*** | |
CN117688785B (zh) | 一种基于种植思想的全张量重力梯度数据反演方法 | |
CN114966878B (zh) | 一种基于混合范数和互相关约束的三维重力场反演方法 | |
CN113468727B (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 |