CN111158059A - 基于三次b样条函数的重力反演方法 - Google Patents
基于三次b样条函数的重力反演方法 Download PDFInfo
- Publication number
- CN111158059A CN111158059A CN202010017136.9A CN202010017136A CN111158059A CN 111158059 A CN111158059 A CN 111158059A CN 202010017136 A CN202010017136 A CN 202010017136A CN 111158059 A CN111158059 A CN 111158059A
- Authority
- CN
- China
- Prior art keywords
- density
- gravity
- inversion
- nodes
- cubic
- 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
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
- G01V7/02—Details
- G01V7/06—Analysis or interpretation of gravimetric records
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
本发明涉及一种基于三次B样条函数的重力反演方法,属于地球物理重力勘探技术领域。其包括如下步骤:确定地下反演区域范围;对反演区进行密度模型参数化,即设置密度节点;计算在某重力观测点处各个节点引起的重力核;循环上一步依次计算完毕所有观测点处各节点引起的重力核,并生成重力核函数矩阵;建立目标函数,求解获得各节点密度值;利用三次B样条函数插值得到整个反演区域的平滑密度场。本发明优点是能够计算平滑密度场重力响应的解析表达,避免了常密度单元重力场叠加的数值解法,解决了现有方法不能够反演空间二阶连续可导的密度场的问题,反演过程由密度节点控制,不需要额外的空间平滑约束,反演效果优于常用的L2模约束反演。
Description
技术领域
本发明涉及一种基于三次B样条函数的重力反演方法,属于地球物理重力勘探技术领域。
背景技术
重力反演分为界面模型反演和密度模型反演两大类,前者是在密度分布规律已知的情况下,反演起伏界面的形态;后者是直接反演地下区域的空间密度分布。在重力的密度模型反演中,人们普遍采用常密度网格剖分的方式,每个网格单元的密度为一常数,将所有常密度单元产生的重力异常累加起来作为某个观测点的重力响应。这种方法的缺点是当网格间距较大时,相邻网格单元的密度会产生突变,无法模拟密度随空间连续可导变化的情况,进而造成模拟重力异常与真实密度场的重力响应存在偏差;若网格间距过小,则造成计算成本和反演多解性增加的问题。
为此,学者们研究了变密度体的重力正反演问题。将密度分布表示为多项式函数的形式,计算密度随空间呈多项式函数变化时产生的重力异常的解析表达式是一种变密度体正演的有效手段。目前,已经发展了密度随深度呈直线、双曲线、抛物线、三阶函数乃至更高阶函数变化时重力异常的正演解析表达形式,但变密度体重力异常反演的发展远落后于变密度体重力异常正演研究。现有的变密度体反演集中在界面模型反演,包括三阶及以下函数形式和三次B样条函数形式下的界面模型反演。现有的变密度体密度模型反演仅限于密度随深度变化的情况,没有考虑密度在水平方向呈多项式连续变化的情况。以上重力反演方法均不能够反演空间二阶连续可导的光滑密度场模型。
发明内容
针对现有技术存在的上述缺陷,本发明提出了一种新的基于三次B样条函数的重力反演方法,能够反演获得空间二阶连续可导的光滑密度场。
本发明是采用以下的技术方案实现的:本发明所述一种基于三次B样条函数的重力反演方法,包括如下步骤:
步骤一:根据重力观测点位置和目标埋深确定地下反演区域范围;
步骤二:对反演区进行密度模型参数化,即设置密度节点:
其中,k表示归一化相对坐标;
步骤三:计算在某重力观测点处各个密度节点引起的重力核:
步骤四:循环步骤三依次计算完毕所有观测点处各节点引起的重力核,并生成重力核函数矩阵:
S6:构建如下的目标函数:
其中,第一项为数据残差项;第二项为先验模型约束项,为重力核函数矩阵,表示先验密度模型,在先验密度差模型信息未知时,该项用零模型约束;为聚焦约束矩阵,为深度加权矩阵,为权重系数;为数据权重矩阵,且为第i个测点观测数据的标准偏差,和分别为密度差的上下限;
S7:目标函数式(5)采用约束最小二乘算法求解,得到L个节点的密度值;
步骤六:利用三次B样条函数和求得的L个节点的密度值插值得到整个反演区域的平滑密度场。
本发明的有益效果是:采用本发明所述的基于三次B样条函数的重力反演方法,通过三次B样条节点对模型进行参数化,能够得到连续密度模型重力异常的精确解析解,避免了常密度单元网格离散造成的重力总场的计算误差,解决了现有方法不能够反演空间二阶连续可导的光滑密度场的问题,反演过程不需要额外的空间平滑约束,反演结果可靠。
附图说明
图1为本发明的流程图;
图2为模型1真实密度模型及其产生的重力异常图;
图3为模型1下常密度单元L2模约束下反演的密度模型图;
图4为模型1下本发明反演的密度模型图;
图5为模型2真实密度模型及其产生的重力异常图;
图6为模型2下常密度单元L2模约束下反演的密度模型图;
图7为模型2下本发明反演的密度模型图。
具体实施方式
为了使本发明目的、技术方案更加清楚明白,下面结合附图,对本发明作进一步详细说明。本发明的流程图,如图1所示,包括如下步骤:
步骤一:根据重力观测点位置和目标埋深确定地下反演区域范围;
步骤二:对反演区进行密度模型参数化,即设置密度节点:
其中,k表示归一化相对坐标;
步骤三:计算在某重力观测点处各个密度节点引起的重力核:
步骤四:循环步骤三依次计算完毕所有观测点处各节点引起的重力核,并生成重力核函数矩阵:
S6:构建如下的目标函数:
其中,第一项为数据残差项;第二项为先验模型约束项,为重力核函数矩阵,表示先验密度模型,在先验密度差模型信息未知时,该项用零模型约束;为聚焦约束矩阵,为深度加权矩阵,为权重系数;为数据权重矩阵,且为第i个测点观测数据的标准偏差,和分别为密度差的上下限;
S7:目标函数式(5)采用约束最小二乘算法求解,得到L个节点的密度值;
步骤六:利用三次B样条函数和求得的L个节点的密度值插值得到整个反演区域的平滑密度场。
下面结合具体实施方式,对于本发明的模型测试进行解释和说明。
实施例一:
为了说明本方法的实现思路及实现过程并证明方法的有效性,用密度异常一正一负的异常体(模型1)进行测试,并和常规常密度单元L2模约束反演的结果进行比较。
S1:将如图2所示的模型1作为真实密度模型。真实密度模型由一正一负的密度异常体组成,两个异常体截面均为边长为1.1 km的正方形,密度差分别为0.5 g/cm3和-0.5g/cm3。
S2:沿水平地表布设重力观测点,测点间距为100 m,共计95个观测点。
S3:计算模型1产生的重力异常,如图2曲线所示。
S4: 将计算的重力异常数据作为观测数据,设置三次B样条节点间距为400 m,共得到19×51=969个节点。
S5:根据公式(3)计算95个重力观测点处各个密度节点引起的重力核,生成大小为95×969的重力核函数矩阵。
S6:建立反演目标函数,利用约束最小二乘算法求解得到969个节点的密度值。
S7:根据969个节点的密度值插值得到整个研究区域连续平滑的密度场模型,如图4所示。
为了说明本发明方法反演的效果,将本发明反演结果与常规常密度单元L2模约束反演结果对比。图3是常规常密度单元L2模约束方法的反演结果。从图中可以明显看出,该方法反演密度幅值偏小,深度越深,密度分布越发散。图4是本发明的反演结果,可以看出本发明反演的异常体的位置,规模大小和密度幅值与真实模型更加接近。
实施例二:
为了进一步说明本方法的实用性,用规模大小各不相同的三个密度异常体(模型2)进行测试,并和常规常密度单元L2模约束反演的结果进行比较。模型2由三个中心埋深和大小不同的正密度异常体组成,三个异常体截面分别为边长1.1 km, 0.6 km, 0.3 km的正方形,密度差均为0.5 g/cm3,如图5所示。计算模型2产生的重力异常如图5曲线所示,实例二反演具体步骤同实例一的S4~S7。图6是常规常密度单元L2模约束方法的反演结果。从图中可以明显看出,该方法反演密度幅值偏小,深度越深,密度分布越发散,异常体边界区分不明显。图7是本发明的反演结果,可以看出本发明反演的异常体的位置,规模大小和密度幅值与真实模型更加接近,异常体边界更清晰。本发明方法优于常规常密度单元L2模约束反演方法。
以上所述仅为本发明的较佳实施例而己,并不以本发明为限制,凡在本发明的精神和原则之内所作的均等修改、等同替换和改进等,均应包含在本发明的专利涵盖范围内。
Claims (5)
1.一种基于三次B样条函数的重力反演方法,其特征在于,包括如下步骤:
步骤一:根据重力观测点位置和目标埋深确定地下反演区域范围;
步骤二:对反演区进行密度模型参数化,即设置密度节点:
其中,k表示归一化相对坐标;
步骤三:计算在某重力观测点处各个密度节点引起的重力核:
(3)
步骤四:循环步骤三依次计算完毕所有观测点处各节点引起的重力核,并生成重力核函数矩阵:
S6:构建如下的目标函数:
其中,第一项为数据残差项;第二项为先验模型约束项,为重力核函数矩阵,表示先验密度模型,在先验密度差模型信息未知时,该项用零模型约束;为聚焦约束矩阵,为深度加权矩阵,为权重系数;为数据权重矩阵,且为第i个测点观测数据的标准偏差,和分别为密度差的上下限;
S7:目标函数式(5)采用约束最小二乘算法求解,得到L个节点的密度值;
步骤六:利用三次B样条函数和求得的L个节点的密度值插值得到整个反演区域的平滑密度场。
2.根据权利要求1所述的基于三次B样条函数的重力反演方法,其特征在于,所述步骤二中,地下任一点的密度用该点周围的16个密度节点通过三次B样条函数插值表示。
3.根据权利要求1所述的基于三次B样条函数的重力反演方法,其特征在于,所述步骤三中,任一重力观测点处某个密度节点引起的重力核利用公式(2)解析计算获得。
4.根据权利要求1所述的基于三次B样条函数的重力反演方法,其特征在于,所述步骤五中的目标函数不需要施加空间平滑约束。
5.根据权利要求1所述的基于三次B样条函数的重力反演方法,其特征在于,所述步骤六中,利用三次B样条函数和节点密度插值可以获得空间二阶连续可导的平滑密度场。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010017136.9A CN111158059B (zh) | 2020-01-08 | 2020-01-08 | 基于三次b样条函数的重力反演方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010017136.9A CN111158059B (zh) | 2020-01-08 | 2020-01-08 | 基于三次b样条函数的重力反演方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN111158059A true CN111158059A (zh) | 2020-05-15 |
CN111158059B CN111158059B (zh) | 2021-04-27 |
Family
ID=70561920
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202010017136.9A Active CN111158059B (zh) | 2020-01-08 | 2020-01-08 | 基于三次b样条函数的重力反演方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN111158059B (zh) |
Cited By (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113468727A (zh) * | 2021-06-11 | 2021-10-01 | 中国海洋大学 | 基于先验结构和已知点双重约束的层状密度建模方法 |
CN113504575A (zh) * | 2021-07-09 | 2021-10-15 | 吉林大学 | 基于权相交及多次交叉梯度约束的联合反演方法 |
CN118094075A (zh) * | 2024-04-28 | 2024-05-28 | 中国地质调查局成都地质调查中心(西南地质科技创新中心) | 基于动态更新加权矩阵的密度模型求取方法 |
CN118094075B (zh) * | 2024-04-28 | 2024-07-09 | 中国地质调查局成都地质调查中心(西南地质科技创新中心) | 基于动态更新加权矩阵的密度模型求取方法 |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105043390A (zh) * | 2015-06-29 | 2015-11-11 | 中国船舶重工集团公司第七0七研究所 | 基于泛克里金法的重力场插值方法 |
US20160139299A1 (en) * | 2014-11-18 | 2016-05-19 | IFP Energies Nouvelles | Method of constructing a geological model |
CN107491411A (zh) * | 2017-06-23 | 2017-12-19 | 中国海洋大学 | 基于n阶多项式密度函数的重力异常反演方法 |
CN108304618A (zh) * | 2018-01-05 | 2018-07-20 | 台州创兴环保科技有限公司 | 一种重力数据与大地电磁数据联合反演方法 |
WO2019058294A1 (en) * | 2017-09-21 | 2019-03-28 | Chevron U.S.A. Inc. | ENHANCED COMPLETE WAVEFORM INVERSION SYSTEM AND METHOD |
-
2020
- 2020-01-08 CN CN202010017136.9A patent/CN111158059B/zh active Active
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20160139299A1 (en) * | 2014-11-18 | 2016-05-19 | IFP Energies Nouvelles | Method of constructing a geological model |
CN105043390A (zh) * | 2015-06-29 | 2015-11-11 | 中国船舶重工集团公司第七0七研究所 | 基于泛克里金法的重力场插值方法 |
CN107491411A (zh) * | 2017-06-23 | 2017-12-19 | 中国海洋大学 | 基于n阶多项式密度函数的重力异常反演方法 |
WO2019058294A1 (en) * | 2017-09-21 | 2019-03-28 | Chevron U.S.A. Inc. | ENHANCED COMPLETE WAVEFORM INVERSION SYSTEM AND METHOD |
CN108304618A (zh) * | 2018-01-05 | 2018-07-20 | 台州创兴环保科技有限公司 | 一种重力数据与大地电磁数据联合反演方法 |
Non-Patent Citations (3)
Title |
---|
PAULA BERKEL · VOLKER MICHEL: "On Mathematical Aspects of a Combined Inversion of Gravity and Normal Mode Variations by a Spline Method", 《MATH GEOSCI》 * |
牛源源等: "近地表密度估计的重力贝叶斯分析方法及在云南地区的应用", 《地球物理学报》 * |
王硕儒等: "二维物性界面深度重磁反演的一种非线性迭代解法-样条函数法", 《海洋与湖沼》 * |
Cited By (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113468727A (zh) * | 2021-06-11 | 2021-10-01 | 中国海洋大学 | 基于先验结构和已知点双重约束的层状密度建模方法 |
CN113468727B (zh) * | 2021-06-11 | 2024-02-06 | 中国海洋大学 | 基于先验结构和已知点双重约束的层状密度建模方法 |
CN113504575A (zh) * | 2021-07-09 | 2021-10-15 | 吉林大学 | 基于权相交及多次交叉梯度约束的联合反演方法 |
CN113504575B (zh) * | 2021-07-09 | 2022-05-03 | 吉林大学 | 基于权相交及多次交叉梯度约束的联合反演方法 |
CN118094075A (zh) * | 2024-04-28 | 2024-05-28 | 中国地质调查局成都地质调查中心(西南地质科技创新中心) | 基于动态更新加权矩阵的密度模型求取方法 |
CN118094075B (zh) * | 2024-04-28 | 2024-07-09 | 中国地质调查局成都地质调查中心(西南地质科技创新中心) | 基于动态更新加权矩阵的密度模型求取方法 |
Also Published As
Publication number | Publication date |
---|---|
CN111158059B (zh) | 2021-04-27 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Korn | Formulation of an unstructured grid model for global ocean dynamics | |
Bodin et al. | Seismic tomography with the reversible jump algorithm | |
CN106199742B (zh) | 一种频率域航空电磁法2.5维带地形反演方法 | |
CN108563837B (zh) | 一种冲积河流水沙模型的模型参数实时校正方法和系统 | |
CN106483559B (zh) | 一种地下速度模型的构建方法 | |
CN110398782A (zh) | 一种重力数据和重力梯度数据联合正则化反演方法 | |
Minniakhmetov et al. | High-order spatial simulation using Legendre-like orthogonal splines | |
Reimond et al. | Spheroidal and ellipsoidal harmonic expansions of the gravitational potential of small Solar System bodies. Case study: Comet 67P/Churyumov‐Gerasimenko | |
CN111158059B (zh) | 基于三次b样条函数的重力反演方法 | |
CN108984939A (zh) | 基于3D Gauss-FFT的三维重力场正演方法 | |
CN109636912A (zh) | 应用于三维声呐图像重构的四面体剖分有限元插值方法 | |
CN113311494A (zh) | 一种卫星重力场反演方法 | |
CN111597753B (zh) | 数据深度变化特征自适应的二维电阻率反演方法及系统 | |
CN114943178A (zh) | 一种三维地质模型建模方法、装置及计算机设备 | |
Mazzia et al. | A comparison of numerical integration rules for the meshless local Petrov–Galerkin method | |
CN113109883B (zh) | 基于等参变换全球离散网格球坐标下卫星重力场正演方法 | |
CN113671570B (zh) | 一种地震面波走时和重力异常联合反演方法与系统 | |
Chen et al. | A high speed method of SMTS | |
CN109239776B (zh) | 一种地震波传播正演模拟方法和装置 | |
CN114139335A (zh) | 基于单松弛时间格子玻尔兹曼模型的粘滞声波模拟方法 | |
Zhang et al. | Bayesian spatial modelling for high dimensional seismic inverse problems | |
CN112748471B (zh) | 一种非结构化等效源的重磁数据延拓与转换方法 | |
CN114155354A (zh) | 基于图卷积网络的电容层析成像重建方法与装置 | |
CN110794469B (zh) | 基于最小地质特征单元约束的重力反演方法 | |
Agoshkov et al. | Recovery of boundary functions on external and internal open boundaries in an open sea hydrodynamic problem |
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 |