CN110333548B - 一种基于归一化异常权函数的高分辨率密度反演方法 - Google Patents
一种基于归一化异常权函数的高分辨率密度反演方法 Download PDFInfo
- Publication number
- CN110333548B CN110333548B CN201910685270.3A CN201910685270A CN110333548B CN 110333548 B CN110333548 B CN 110333548B CN 201910685270 A CN201910685270 A CN 201910685270A CN 110333548 B CN110333548 B CN 110333548B
- Authority
- CN
- China
- Prior art keywords
- inversion
- function
- density
- underground
- weight
- 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.)
- Expired - Fee Related
Links
- 238000000034 method Methods 0.000 title claims abstract description 49
- 230000002159 abnormal effect Effects 0.000 title claims abstract description 25
- 230000005484 gravity Effects 0.000 claims abstract description 49
- 239000011159 matrix material Substances 0.000 claims abstract description 19
- 238000009826 distribution Methods 0.000 claims abstract description 14
- 230000008569 process Effects 0.000 claims abstract description 11
- 238000010606 normalization Methods 0.000 claims description 8
- 230000000694 effects Effects 0.000 claims description 7
- 238000012545 processing Methods 0.000 claims description 7
- 238000003825 pressing Methods 0.000 claims description 3
- 230000002547 anomalous effect Effects 0.000 claims 1
- 230000006870 function Effects 0.000 description 48
- XEEYBQQBJWHFJM-UHFFFAOYSA-N Iron Chemical compound [Fe] XEEYBQQBJWHFJM-UHFFFAOYSA-N 0.000 description 22
- 229910052742 iron Inorganic materials 0.000 description 11
- 238000009933 burial Methods 0.000 description 6
- 238000010586 diagram Methods 0.000 description 6
- 239000011435 rock Substances 0.000 description 6
- 230000005856 abnormality Effects 0.000 description 3
- BVKZGUZCCUSVTD-UHFFFAOYSA-L Carbonate Chemical compound [O-]C([O-])=O BVKZGUZCCUSVTD-UHFFFAOYSA-L 0.000 description 2
- 238000002939 conjugate gradient method Methods 0.000 description 2
- 238000007796 conventional method Methods 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 238000000926 separation method Methods 0.000 description 2
- 238000001228 spectrum Methods 0.000 description 2
- 230000001629 suppression Effects 0.000 description 2
- 229910052612 amphibole Inorganic materials 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 229910052500 inorganic mineral Inorganic materials 0.000 description 1
- 239000011707 mineral Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
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
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/15—Correlation function computation including computation of convolution operations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Data Mining & Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- Computational Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Computing Systems (AREA)
- Geophysics (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公开了一种基于归一化异常权函数的高分辨率密度反演方法,包括:具体反演过程将地下三维空间分为有限个密度一定的长方体单元,由地下密度分布所产生的异常近似为所有地下块体单元分别对地表观测点产生的异常的叠加,当存在n个观测点和m个长方体模型单元时,重力异常正演表达式用矩阵相乘的方式来表示。本发明提出归一化异常权和深度权相结合的密度反演方法,能有效地提高反演结果的水平和垂直分辨率,且由于在获取不同深度异常特征过程中主要依据向上延拓函数进行,因此反演过程更加稳定,能有效地压制噪声的干扰。
Description
技术领域
本发明涉及密度反演方法,具体涉及一种基于归一化异常权函数的高分辨率密度反演方法。
背景技术
目前,三维密度反演是重力数据反演方法中反演目标最为全面的方法。因为密度反演是一个求解欠定方程的问题,也就无法避免地球物理反演中的多解性问题,如何通过重力场的物理性质以及地质条件等密度分布必须遵守的规律来限制反演结果,从而减小反演的多解性,获得更加符合实际的密度分布特征至关重要。
发明内容
本发明的主要目的在于提供一种基于归一化异常权函数的高分辨率密度反演方法。
本发明采用的技术方案是:一种基于归一化异常权函数的高分辨率密度反演方法,包括:
具体反演过程将地下三维空间分为有限个密度一定的长方体单元,由地下密度分布所产生的异常近似为所有地下块体单元分别对地表观测点产生的异常的叠加,当存在n个观测点和m个长方体模型单元时,重力异常正演表达式用矩阵相乘的方式来表示:
核函数矩阵为:
通过以下形式来求解此欠定问题:
进一步地,基于归一化异常权函数的高分辨率密度反演方法还包括:
在公式(6)中将权函数中的原始异常g替换为异常的垂向导数,反演结果还会更加收敛:
其中,为相应i层深度对应高度的异常值,不同高度的值通过延拓的方式计算获得;为重力原始异常的最大值和最小值,作用是对整体数据进行归一化处理,避免不同的重力数据整体数值的大小对反演产生影响;a为常数,为了均衡异常权与其他权之间的数量级;的作用是将括号内的一维数组变为n×n的对角矩阵,只有对角线上元素与A中的元素一一对应,其余部分为0;k为地下块体纵向上分割的个数。
更进一步地,基于归一化异常权函数的高分辨率密度反演方法还包括:
(8)式的矩阵形式为:
本发明的优点:
本发明针对现有重力密度反演在同一平面上权函数相同所导致的结果分辨率低的问题,提出归一化异常权和深度权相结合的密度反演方法,新型的权函数是依据异常与深度的对应关系来建立,因此所提出方法能有效地提高反演结果的水平和垂直分辨率,且由于在获取不同深度异常特征过程中主要依据向上延拓函数进行,因此反演过程更加稳定,能有效地压制噪声的干扰。通过模型正演异常对所提出密度反演方法开展试验,反演结果表明该方法相对常规方法能获得更加准确和收敛的结果,且对于噪声具有良好的压制效果,并且利用重力垂向梯度异常进行加权时会取得更加高的分辨率。将该方法应用于山东地区的铁矿实测重力异常的反演,获得铁矿的空间分布特征,为下一步精细勘探提供了重要的支撑。
除了上面所描述的目的、特征和优点之外,本发明还有其它的目的、特征和优点。下面将参照图,对本发明作进一步详细的说明。
附图说明
构成本申请的一部分的附图用来提供对本发明的进一步理解,本发明的示意性实施例及其说明用于解释本发明,并不构成对本发明的不当限定。
图1是本发明的基于深度权的三维密度反演结果切片图;
图2是本发明的不同深度权函数的获得方式示意图;
图3是本发明的矿区实测重力异常图;
图4是本发明的总水平导数及断裂划分图;
图5是本发明的实测重力异常的平均对数功率谱图;
图6是本发明的位场分离后矿区重力异常图;
图7是本发明的三维密度反演结果高值三维成像图;
图8是本发明的估计铁矿范围解释图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白,以下结合附图及实施例,对本发明进行进一步详细说明。应当理解,此处所描述的具体实施例仅仅用以解释本发明,并不用于限定本发明。
本发明的具体反演过程是将地下三维空间分为有限个密度一定的长方体单元,由地下密度分布所产生的异常近似为所有地下块体单元分别对地表观测点产生的异常的叠加,当存在n个观测点和m个长方体模型单元时,重力异常正演表达式用矩阵相乘的方式来表示:
通过Holstein(1999)提出的长方体模型的正演计算知,核函数矩阵为:
图1为本发明基于深度权的三维密度反演结果切片图,(a)为重力异常剖面; (b)为深度权反演结果切片。
图1中为只用深度加权的三维密度反演结果的切片图,白色虚线为模型体的实际位置。从图中看出,由于核函数矩阵G随着深度的增加快速衰减,所以反演结果在深部具有较低的分辨率,引入下半空间的空间分布权函数,从而获得较高水平和垂直分辨率的结果。
图2为本发明不同深度权函数的获得方式示意图,(a)为原始异常向上延拓5m;(b)为原始异常向上延拓35m;(c)为原始异常向上延拓65m;(d)为原始异常向上延拓95m。
因为不同的重力异常的极值是不同的,所以为了去除异常极值大小带来的影响,将异常进行归一化处理,这样处理后的结果仅能体现异常高低的特征。又因为实测重力异常显示的是地下的综合结果,为了使权函数在深部时能主要体现深部的重力异常特征,利用向上延拓的方式来确定不同深度的权函数所利用的重力异常。如图2所示,
将向上延拓5m的异常数据的计算结果作为5m深处的权函数,按照这种方式形成一个和核函数相对应的对角矩阵。具体形式如下:采用不同高度的延拓异常进行归一化处理,从而反演结果将更加收敛;而不会像以往方法一样产生十分发散的结果,因此以往方法反演结果发散的区域异常权函数会使其无变化,接近零值,从而结果将更加收敛;
在公式(6)中将权函数中的原始异常g替换为异常的垂向导数,反演结果还会更加收敛:
其中,为相应i层深度对应高度的异常值,不同高度的值通过延拓的方式计算获得;为重力原始异常的最大值和最小值,作用是对整体数据进行归一化处理,避免不同的重力数据整体数值的大小对反演产生影响;a为常数,为了均衡异常权与其他权之间的数量级;的作用是将括号内的一维数组变为n×n的对角矩阵,只有对角线上元素与A中的元素一一对应,其余部分为0;k为地下块体纵向上分割的个数;将(5)式、(7)式结合后代入到(3)式中的权函数里,得到具体的反演函数为:
(8)式的矩阵形式为:
实际数据处理:
为了验证的方法在实际数据中的效果,将提出的方法应用于山东铁矿矿区的实测重力数据中。该地区位于华北板块,地处鲁西与鲁中隆起的西北部。地表存在大面积的第四系覆盖,并且隐伏碳酸盐岩分布广泛。附近已经发现的矿床资料显示,该地区主要为矽卡岩型铁矿床。矿体主要在闪长岩体与碳酸盐岩的接触带上,所以侵入岩的接触带为重要的找矿位置,为高密度体。重力数据主要是通过寻找重力异常高值的侵入岩来判断铁矿位置。将地下划分为40×30×20个长方形块体,每个长方体的大小为750×1500×500m。对数据重新进行网格化,获得观测数据个数为40×30个,观测点坐标对应长方体中心位置。实测异常如图3所示。我们对实测重力数据进行三维密度反演,并筛选出反演结果中的高密度值,来判断铁矿位置与范围。
因为断裂可以为岩浆上涌提供通道,矿体所在的岩浆侵入的位置经常会存在断裂构造,首先对原始重力异常进行总水平导数的计算来划分出当地的断裂构造,如图4所示,图中黑色虚线为断裂位置,即总水平导数的高值区域。
因为重力异常为地下所有密度的综合相应,所以对于实测异常中存在的背景场为干扰反演的异常。通过匹配滤波法对实测异常进行位场分离,图5为异常的平均对数功率谱,将波数较小的区域滤掉,获得浅部的剩余重力异常,如图5所示。
从图6的重力异常中我们可以看出,重力高值区大致分为六块,其中Ⅰ、Ⅲ、Ⅵ处的重力异常较高,说明这三处的高密度体范围较大,图6中西南角的高值不完全在勘探区域中,所以无法反演出准确的结果,所以本次处理并不考虑此高值异常。利用提出的方法对此异常进行反演,获得的反演结果如图7所示。
图7为反演结果中高密度块体的三维成像,表面覆盖的为原始重力异常,并且按照图4所示用罗马数字标出了六块高值异常的位置,很好的对应了反演结果中的六块地下高密度体,并且矿体倾向与异常显示完全吻合,从图7中可以看出,总共有六个区域是高密度块体汇聚的地区,这些区域的铁矿分布不是相互连通的,通过对高密度块体的坐标可以获得各个区域铁矿的埋深与范围,最终得到图8。
图8展示的是矿体的水平范围,即该地区的三维密度反演结果解释图。通过图7的三维结果可以获得各处矿体的顶面埋深,其中图8中Ⅰ和Ⅴ的顶面埋深为800m;Ⅲ和Ⅵ的顶面埋深为1000m;Ⅱ的顶面埋深为1200m;Ⅳ的顶面埋深为1400m。
本发明针对现有重力密度反演在同一平面上权函数相同所导致的结果分辨率低的问题,提出归一化异常权和深度权相结合的密度反演方法,新型的权函数是依据异常与深度的对应关系来建立,因此所提出方法能有效地提高反演结果的水平和垂直分辨率,且由于在获取不同深度异常特征过程中主要依据向上延拓函数进行,因此反演过程更加稳定,能有效地压制噪声的干扰。通过模型正演异常对所提出密度反演方法开展试验,反演结果表明该方法相对常规方法能获得更加准确和收敛的结果,且对于噪声具有良好的压制效果,并且利用重力垂向梯度异常进行加权时会取得更加高的分辨率。将该方法应用于山东地区的铁矿实测重力异常的反演,获得铁矿的空间分布特征,为下一步精细勘探提供了重要的支撑。
以上所述仅为本发明的较佳实施例,并不用以限制本发明,凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。
Claims (2)
1.一种基于归一化异常权函数的高分辨率密度反演方法,其特征在于,包括:
具体反演过程将地下三维空间分为有限个密度一定的长方体单元,由地下密度分布所产生的异常近似为所有地下块体单元分别对地表观测点产生的异常的叠加,当存在n个观测点和m个长方体模型单元时,重力异常正演表达式用矩阵相乘的方式来表示:
核函数矩阵为:
通过以下形式来求解此欠定问题:
所述的基于归一化异常权函数的高分辨率密度反演方法,还包括:
在公式(6)中将权函数中的原始异常g替换为异常的垂向导数,反演结果还会更加收敛:
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910685270.3A CN110333548B (zh) | 2019-07-27 | 2019-07-27 | 一种基于归一化异常权函数的高分辨率密度反演方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910685270.3A CN110333548B (zh) | 2019-07-27 | 2019-07-27 | 一种基于归一化异常权函数的高分辨率密度反演方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN110333548A CN110333548A (zh) | 2019-10-15 |
CN110333548B true CN110333548B (zh) | 2021-01-29 |
Family
ID=68147681
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201910685270.3A Expired - Fee Related CN110333548B (zh) | 2019-07-27 | 2019-07-27 | 一种基于归一化异常权函数的高分辨率密度反演方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN110333548B (zh) |
Families Citing this family (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112666607B (zh) * | 2019-10-16 | 2024-06-25 | 中国石油天然气集团有限公司 | 重力反演黄土层厚度分布的方法及装置 |
CN113514899A (zh) * | 2021-07-12 | 2021-10-19 | 吉林大学 | 一种重力梯度的自适应剖分反演方法 |
Family Cites Families (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2007102973A2 (en) * | 2006-03-08 | 2007-09-13 | Exxonmobil Upstream Research Company | Efficient computation method for electromagnetic modeling |
CN105549106B (zh) * | 2016-01-07 | 2018-06-08 | 中国科学院地质与地球物理研究所 | 一种重力多界面反演方法 |
CN109507749A (zh) * | 2019-01-14 | 2019-03-22 | 中国地质调查局成都地质调查中心 | 一种重磁自约束三维反演与联合解释方法 |
-
2019
- 2019-07-27 CN CN201910685270.3A patent/CN110333548B/zh not_active Expired - Fee Related
Also Published As
Publication number | Publication date |
---|---|
CN110333548A (zh) | 2019-10-15 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN101285896B (zh) | 一种地球物理勘探中的重磁数据处理方法 | |
CN110333548B (zh) | 一种基于归一化异常权函数的高分辨率密度反演方法 | |
Xia et al. | Resolution of high-frequency Rayleigh-wave data | |
TSELENTIS et al. | A spectral approach to Moho depths estimation from gravity measurements in Epirus (NW Greece) | |
CN110907995B (zh) | 井中vsp地震数据的逆时偏移方法及装置 | |
CN112465191A (zh) | 隧道突涌水灾害预测方法、装置、电子设备及存储介质 | |
Yin et al. | Joint inversion of Rayleigh and Love wave dispersion curves for improving the accuracy of near-surface S-wave velocities | |
Malagnini et al. | Estimating absolute site effects | |
Milano et al. | The deep crust beneath the Trans‐European Suture Zone from a multiscale magnetic model | |
CN109902315B (zh) | 一种圈定隐伏花岗岩岩体深部边界的方法 | |
Zhu et al. | Seismic sedimentology of sand-gravel bodies on the steep slope of rift basins—A case study of the Shahejie Formation, Dongying Sag, Eastern China | |
CN110361788B (zh) | 一种空-地联合三维重力数据特征分析及密度反演方法 | |
AU2011286367A1 (en) | Obtaining a response based on differencing of outputs of sensors | |
CN101609163B (zh) | 基于波动理论的多尺度地震资料联合成像方法 | |
Czarny et al. | The application of seismic interferometry for estimating a 1D S-wave velocity model with the use of mining induced seismicity | |
CN108919351A (zh) | 基于逆时聚焦原理进行观测系统双向聚焦性的评价方法 | |
Šumanovac | Mapping of thin sandy aquifers by using high resolution reflection seismics and 2-D electrical tomography | |
Li‐Min et al. | Crustal magnetic anomalies and geological structure in the Yunnan region | |
EP2430481A1 (en) | Geophysical data processing systems | |
CN108375794B (zh) | 基于对称观测的vsp缝洞绕射成像技术方法 | |
CN107179548B (zh) | 一种基于真地表的叠前地震成像方法 | |
Sarris | Processing and Analysing Geophysical Data | |
CN105842731B (zh) | 基于波场延拓原理的起伏地表组合震源波场定向方法 | |
Hatanaka et al. | Three-dimensional modeling and inversion of the mise-a-la-masse data using a steel-casing borehole | |
dos S Filho et al. | The architecture of the Eopaleozoic Cococi Basin, northeastern Brazil: 2D geological modelling from magnetic and gravimetric data |
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 | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20210129 Termination date: 20210727 |
|
CF01 | Termination of patent right due to non-payment of annual fee |