双重约束下的密度三维反演

刘展，于会臻，陈挺

(1.中国石油大学地球科学与技术学院，山东青岛266555;2.四川省地矿局物探队，四川成都610072)

双重约束下的密度三维反演

刘展1，于会臻1，陈挺2

(1.中国石油大学地球科学与技术学院，山东青岛266555;2.四川省地矿局物探队，四川成都610072)

为增加密度三维反演过程中的约束及避免“上漂”现象，从定性与定量相结合的角度提出一种密度相关概率成像与基于钻井的位场特征约束相结合的双重约束机制，重点提出以核函数均值差作为扫描函数，并给出依据剖分网格密度属性与位场分布相关性差异进行密度值双重约束的具体方案。理论模型试验结果表明，利用所提方法在先得到地下剖分单元的密度定性分布特征的情况下，可给出定量反演过程中的剖分网格密度取值范围及其搜索顺序，尤其在有钻井资料的区域中反演结果的纵向分辨率得到了较大提高。实际资料处理结果进一步验证了方法的实用性。

地球物理勘探;重磁勘探;密度反演;成像技术;双重约束;黄金分割;约束扩展;赋存空间

重磁反演模型可以大致分为两类:形态模型和物性模型［1］。密度反演属于物性反演，具有严重的多解性问题。为此，地球物理工作者做了大量工作。张贵宾等［2］建立了位场广义线性综合反演系统，在重磁异常线性反演中将该理论与吉洪诺夫正则化方法相结合求解地下密度源分布及质心位置;Bear等［3］使用线性反演方法，在源所在范围内，使用最小二乘法反演求解，根据连续的Household变换进行矩阵正交分解，在迭代过程中使用奇异值分解方法，从而得到重力资料的三维密度分布;Oleg Portniaguine［4］基于Tikhonov正则化方法，使用最小稳定方程得到聚集的反演成像，迭代反演过程使用权重模型参数构建，从而使方法收敛并稳健;管志宁等［5］把神经网络与重磁异常反演理论相结合，提出了一种用于重磁反演的拟BP神经网络方法;还有许多学者［6-15］也提出了许多不同的方法。上述文献中许多密度反演的计算方法注重了求解过程中搜索、迭代优化方法的改进，但对优化场源体密度约束的研究并不多，已知信息利用率较低，反演结果不稳定。为此，笔者从约束范围扩展角度提出双重约束机制。

1 方法原理

通常密度反演方法首先将地下场源区域规则地划分成若干小长方体网格，如图1(a)所示，将场源可能分布区域Ω剖分为若干个规则网格单元，设XYZ方向的剖分间距分别为Δξ、Δγ、Δζ，网格任意节点坐标为(ξi，ηj，ζk)，其中ξi=iΔξ，ηj=jΔη，ζk=kΔζ(i=－N，N;j=－M，M;k=－K，K)，图1(b)为第ijk个剖分单元示意图，设其密度值为σijk;之后，通过调整各单元的密度值并计算正演场数据来拟合实测数据;最终反演确定各网格的密度最优解。

本文中密度定量反演方法采用文献［16］中提出的以单相关系数为基础构建目标函数的多维黄金分割法。其关键技术是引入多元线性回归分析中的相关系数R。R的物理意义是模型理论场与观测场之间的误差方差与观测场本身方差的逼近程度，其目的是求取模型与观测数据的拟合组成的目标函数的极小来获得密度值的解。R的取值为0～1，越接近1说明观测场与核函数的相关性越好，反之则越差，使得R最接近1的密度模型就是反演最优解。算法在解的逼近问题上采用n维黄金分割最优化方法，密度反演中模型修改更合理、速度更快，但仍存在以下约束问题:

(1)剖分网格密度初始约束区间的确定问题。可知，约束区间范围越精确，反演速度越快，解越稳定。但是由于已知资料往往较少，通常方法仅对钻遇网格给出较小的密度约束区间，多数未知网格往往只根据地球物理资料统一取一个比较大的约束范围，致使密度值搜索速度慢且不稳定。

(2)网格反演顺序的确定问题。通常网格反演计算顺序按由浅到深逐层进行，这样容易产生上层剖分网格占据大部分异常值的“上漂”现象。

事实上，由于地下地质体赋存空间内密度值必然存在着一定的连续性且在所引起的重力场相关特征中有所体现，以往简单的钻井点约束可利用该特征对附近网格进行密度约束空间的扩展。同时，剖分网格与观测场之间有一定的映射关系，找到这种关系可确定任一网格对于整个场的贡献值，按照贡献值大小进行反演可使存在场源体可能性大的区域先分配场值从而避免“上漂”。因此，按上述思想，本文中提出以下双重约束机制:

(1)利用密度概率成像约束方法计算每个网格的密度分布概率值，根据概率值的大小判断网格存在场源的可能性以确定反演顺序，同时对概率值进行阈值处理并对大于阈值的网格搜索得到场源体成像。

(2)利用基于钻井的位场特征进一步确定水平方向上密度分布情况，结合步骤(1)的场源体成像结果，分区域将钻遇网格的密度值扩展到邻近剖分网格中，以增加反映约束信息，把钻井所在位置的串状点约束变为“糖葫芦”状体约束。

1.1 密度相关概率成像

密度相关概率成像技术研究的是场源体之间的相对赋存关系，扫描函数的选取是概率成像技术的研究重点。由于某些理论上成像效果较好的扫描函数可能会引入较大的计算误差，实际应用效果并不理想，如导数扫描函数等。为避免该问题，选用的扫描函数应具有计算简单且可在背景场中突出目标异常源的特点。因此，针对重磁数据的特征，本文中提出利用核函数均值差作为扫描函数并给出了成像方案。

概率成像原理可见文献［17］～［21］，现以图1为例求改进的重力场扫描函数。

设Δg(r)为在r=(x，y，z)处的实测重力异常，其与观测异常均值之差为

式中，G为万有引力常数;N为观测点数。定义扫描窗口S，对式(1)在S内做平方积分得

取第q部分，并利用施瓦茨不等式，得

定义改进的重力场扫描函数为

定义密度相关概率成像函数为

其中

现以图2(箭头为钻井位置)为例介绍场源体成像方案。

(1)计算密度相关概率值。图2(a)为利用核函数均值差作为扫描函数得到的概率值，该值反映了剖分网格相对实测异常的贡献，可确定最终密度定量反演时剖分网格的计算顺序。

(2)阈值处理。假设根据地质情况选阈值为0.6，概率值大于0.6的赋1，如图2(b)所示。

图2 成像概率值搜索结果网格示意图Fig.2 Grid schemes of imaging probability value search results

(3)场源体二维成像。以第一个网格a1b1作为扫描起点，逐行找到第一个为1的网格，如网格a4b2。之后采用八方位邻域搜索算法搜索相邻为1的网格，圈定如图2(b)中的4个地质体边界。

(4)场源体三维成像。将相邻两个剖面的地质体单体成像概率值突变点作为边界，确定场源体三维空间分布。

1.2 基于钻井的位场特征约束

实际地质体密度值在局部空间内具有一定连续性且场源体密度横向分布与观测场具有一定的相关性，因此可利用基于钻井的位场特征约束对井周围网格密度分布区间进一步划分。

设任一钻遇场源点(ξ，γ，ζ)在计算点(x，y，z)处的位场统一写成如下形式:

式中，ΔE为位场异常;C为不依赖于xyz的常数。当场源与观测面足够远或者场源足够小时，近似为点源，N0=2。

分别对式(7)中x、y、z求导，得

表1是通过模型试验研究得到的重力异常Δg、水平一阶导数Δgx、垂向一次导数Δgz和垂向二次导数Δgzz的零值水平区间距L与场源体横向分布L0之间的关系。

表1 异常特征值水平范围与异常体宽度的关系Table 1 Relationship between abnormal eigenvalue level and abnormal body width

由表1可知，L/L0的值越接近1，说明位场特征对异常体水平边界刻画越准确。因此，同一埋深时，利用重力异常垂向二次导数的零值范围圈定场源体横向赋存空间最准确。按照该结论可对钻井附近的剖分网格密度值约束区间进一步扩展，视位场特征边界范围内的网格为密度值接近的场源体，为下一步对剖分网格的区域划分提供依据。

1.3 密度约束范围扩展

以场源体概率成像为主，以基于钻井的位场特征约束范围为辅，分4种情形对钻遇网格的密度值进行定量约束扩展，以图2模型为例，图3为其位场特征三维示意图，设钻遇网格a5b2的密度取值范围为［α，β］。

图3 钻井定量约束扩展示意图Fig.3 Sketch map of drilling quantitative constraint expansion

(1)钻遇概率成像圈定的场源体。此时假设圈定的地质体内密度值相同，采用紧约束。图2中网格a4b2、a6b2、a4b3、a5b3、a6b3的密度值都与a5b2相同，为［α，β］。

(2)在钻井的位场特征范围之内但并非概率成像圈定的剖分网格。假设密度在位场特征范围内渐进变化，采用渐进约束，即以图2(a)所示网格的初始概率值的比值关系对未知网格赋初值，所赋值均值与已知网格密度均值呈正比关系，上下限与已知网格范围相同，则a3b2取值范围为［0.2((β+α)/2)/0.9－(β－α)/2，0.2((β+α)/2)/0.9+(β－α)/2］;a3b3取值范围为［0.3((β+α)/2)/0.9－(β－α)/2，0.3((β+α)/2)/0.9+(β－α)/2］。若其下限小于0，则赋为0。

(3)未钻遇但处于钻井垂向延长线上的场源体剖分网格。设密度的深度分布呈指数函数形式，即σ=σ0eAΔh，Δh为深度差，A为常数，则图2中a5b5、a6b5、a7b5的密度值均为［(β+α)/2exp(AΔh)－(β－α)/2，(β+α)exp(AΔh)/2+(β－α)/2］，此处Δh为a5b5与a5b2的深度差。

(4)其他未知网格。此时与常规方法采用同样的约束方法，即按照研究区内地球物理密度参数对其赋初值。

2 理论模型试验

构建两个模型对上述方法原理进行验证，测网为2500 m×2000 m，测点距50 m，按100 m×100 m×100 m进行网格剖分。

(1)模型a模拟两个水平位置不同的密度体横向叠加异常，两个地质体中心横向间距为600 m。钻井位于(350，0)处，深300 m处密度值为0.5 g/cm3。模型参数如表2所示。反演中概率成像阈值选为0.6，反演结果剖面等值线如图4所示。

表2 模型a、b参数Table 2 Parameters of model a and b

图4 模型a反演成像结果剖面等值线Fig.4 Section contour line of model a inversion imaging results

(2)模型b模拟两个深度不同的地质体纵向叠加异常，模型中心纵向间距为300 m。钻井位于(0，0)处，在200、500 m处密度值为0.5 g/cm3。模型参数如表2所示。反演中概率成像阈值选为0.6，反演结果剖面等值线如图5所示。

图5 模型b反演成像结果剖面等值线Fig.5 Section contour line of model b inversion imaging results

模型反演结果精度如表3所示。

表3 模型反演结果相对误差Table 3 Relative error of model inversion results%

将模型(b-2)在x方向上变宽至－200～200 m时，反演精度降为3.13%，尤其远离钻井的网格密度反演结果精度下降明显，这是由于离钻遇网格越远密度约束值定量扩展所得结果就越不准确。因此，误差会随着场源体体积的增大而变大，但相比常规反演，该反演方法仍具有更高的可信度。

模型试验表明:就密度横向分布来说，本文中反演结果边界更加清晰，误差较小;对于纵向分布，在有钻井的情况下，叠加的地质体得到有效分离，纵向埋深归位准确，无钻井时反演结果与真实值之间仍有一定误差，但相对常规方法有较大改善。

3 实测资料处理

为验证方法实用性，选某研究区为试验区。该区火成岩比较发育，主要为侵入性的辉绿岩、喷发性的玄武岩。已知研究区内有6口井Y27、Y30、Y101、Y1、Y20、Y29钻遇火成岩，选取其中3口井(Y27、Y30和Y101)的资料作已知先验约束，另外3口井(Y1、Y20和Y29)的资料作为结果验证。采用插值切割方法计算火成岩体引起剩余重力异常(图6(a));采用半径为10 km的四阶趋势面分析计算，利用原始化极异常减去趋势面分析场计算所得的磁异常即为火成岩体引起的磁力异常(图6(b))。

图6 火成岩引起的重磁异常等值线Fig.6 Gravity and magnetic anomaly contour caused by igneous rock

依据本文方法反演得到研究区内密度值分布(图7)，其中概率阈值选为0.6。由得到重磁三维定量反演的视密度和视磁化强度数据体确定火成岩分布情况(图8)，其中火成岩圈定的物性值下限视密度为0.1 g/cm3，视磁化强度为0.5 A/m。

将反演结果与Y1、Y20、Y29三口井资料进行多深度综合对比，统计所得视密度相对误差分别为4.02%、2.53%、3.54%，火成岩深度相对误差为3.48%、4.11%、3.01%。可见双重约束方法求得视密度三维空间分布精度较高。由图8看出，研究区内火成岩较发育，多发育在沙三段、沙四段。侵入火成岩体多发育于Y1北部、Y14、Y27和Y29等井区附近，在Y29与Y1井之间，发育了多套火成岩体。由实际资料处理可知，双重约束方法圈定出了火成岩体三维空间位置，展示了其空间赋存状况。

图7 不同深度切片火成岩密度等值线Fig.7 Density contour lines for different depth igneous rock petrophysics slice

图8 不同深度火成岩圈定图Fig.8 Delineating different depth location of igneous rocks

4 结束语

针对实际应用中重磁物性三维定量出现的约束不充分的问题，从约束范围扩展角度提出了双重约束机制。关键采用了两种手段:一是采用核函数均值差作为扫描函数，其更加适合密度反演成像，同时可为后续定量反演提供网格计算顺序，在一定程度上改善了“上漂”现象;二是将基于钻井的位场特征研究与成像结果结合起来，并分4种有效情况对所有剖分网格进行定量约束。利用该方法进行密度反演在一定程度上改善了常规方法中的已知约束利用率低的问题，提高了反演计算速度和合理性。但是，

还需注意:概率阈值应在对比反演结果与已知资料后进行最优选择;若能对不同深度范围内地质体引起的异常进行更有效地分离，则会使反演结果更加合理;在研究区内无钻井区域内需通过结合其他物探资料来提高反演精度。

［1］姚长利，郝天珧，管志宁.重磁反演约束条件及三维物性反演技术策略［J］.物探与化探，2002，26(1):253-256.YAO Chang-li，HAO Tian-yao，GUAN Zhi-ning.Restrictions in gravity and magnetic inversions and technical strategy of 3D properties inversion［J］.Geophysical＆Geochemical Exploration，2002，26(1):253-256.

［2］张贵宾，申宁华，王喜臣，等.位场广义线性综合反演系统的建立［J］.长春地质学院学报，1993，23(2):197-204.ZHANG Gui-bin，SHEN Ning-hua，WANG Xi-chen，et al.Location-field generalized linear inversion system establishment［J］.Journal of Changchun University of Earth Science，1993，23(2):197-204.

［3］BEAR G W，AL-SHUKRI H J，RUDMAN A J.Linear inversion of gravity data for 3-D density distributions［J］.Geophysics，1995，60(5):1354-1364.

［4］PORTNIAGUINE O，ZHDANOV M S.3-D magnetic inversion with data compression and image focusing［J］.Geophysics，2002，67(5):1532-1541.

［5］管志宁，侯俊胜，黄临平，等.重磁异常反演的拟BP神经网络方法及其应用［J］.地球物理学报，1998，41(2):242-251.GUAN Zhi-ning，HOU Jun-sheng，HUANG Lin-ping，et al.Inversion of gravity and magnetic anomalies using pseduo-BP neural network method and its application［J］.Chinese Journal of Geophysics，1998，41(2):242-251.

［6］刘展，王万银.曲面位场的正则化线性规划直接反演技术［J］.地质与勘探，1999，35(3):62-66.LIU Zhan，WANG Wan-yin.A method of regularized linear programming applying to directly inversion of potential fields on A C［J］.Geology and Exploration，1999，35(3):62-66.

［7］王妙月，底青云，许琨，等.磁化强度矢量反演方程及二维模型正反演研究［J］.地球物理学报，2004，47(3):528-534.WANG Miao-yue，DI Qing-yun，XU Kun，et al.Magnetization vector inversion equations and 2D forward and inversed model study［J］.Chinese Journal of Geophysics，2004，47(3):528-534.

［8］杨辉.重力、地震联合反演基岩密度及综合解释［J］.石油地球物理勘探，1998，33(1):496-502.YANG Hui.Basement density inversion using gravimetricand seismic data and the integrative interpretation［J］.Oil Geophysical Prospecting，1998，33(1):496-502.

［9］柯小平，王勇，许厚泽.用遗传算法反演地壳的变密度模型［J］.武汉大学学报，2004，29(11):981-984.KE Xiao-ping，WANG Yong，XU Hou-ze.Inversion of variable density model of crust from genetic algorithms［J］.Geomatics and Information Science of Wuhan University，2004，29(11):981-984.

［10］MONTESIONS F G，ARNOSO J，VIEIRA R.Using a genetic algorithm for 3-D inversion of gravity data in Fuerteventura［J］.Int J Sci(Geol Rundsch)，2005，94(2):301-316.

［11］柯小平，王勇，许厚泽.三维密度分布的遗传算法反演［J］.大地测量与地球动力学，2009，29(1):41-45.KE Xiao-ping，WANG Yong，XU Hou-ze.3D density inversion with genetic algorithm［J］.Journal of Geodesy and Geodynamics，2009，29(1):41-45.

［12］王德人.非线性方程组解法与最优化方法［M］.北京:人民教育出版社，1980:13-17.

［13］陈东敬，张新兵.带模拟退火的拟BP神经网络在伊朗某地区重力资料反演中的应用［J］.勘探地球物理进展，2005，28(3):215-218.CHEN Dong-jing，ZHANG Xin-bing.Application of SABPNN in inversion of gravity data from Iran［J］.Progress in Exploration Geophysics，2005，28(3):215-218.

［14］裴正林，余钦范，牟永光.小波多尺度井间地震层析成像方法［J］.地球学报，2002，23(1):383-386.PEI Zheng-lin，YU Qin-fan，MU Yong-guang.Crosshole seismic tomography by wavelet multiscale［J］.Acta Geosicientia Sinica，2002，23(1):383-386.

［15］姚长利，郑元满.重磁异常三维物性反演随机子域法方法技术［J］.地球物理学报，2007，50(2):1576-1583.YAO Chang-li，ZHENG Yuan-man.3-D gravity and magnetic inversion for physical properties using stochastic subspaces［J］.Chinese Journal of Geophysics，2007，50(2):1576-1583.

［16］刘展，班丽，魏巍，等.济阳坳陷花沟地区火成岩重磁成像解释方法［J］.中国石油大学学报:自然科学版，2007，31(1):30-34.LIU Zhan，BAN Li，WEI Wei，et al.A method of inversing igneous rock by gravity and magnetic imaging in Huagou area of Jiyang depression［J］.Journal of China University of Petroleum(Edition of Natural Science)，2007，31(1):30-34.

［17］许令周，关继腾，房文静.高次导数的概率成像原理［J］.青岛大学学报，2003，16(1):32-36.XU Ling-zhou，GUAN Ji-teng，FANG Wen-jing.Theory of probability tomography about second derivative formula［J］.Journal of Qingdao University，2003，16(1):32-36.

［18］徐世浙，余海龙，李海侠，等.基于位场分离与延拓的视密度反演［J］.地球物理学报，2009，52(6):1592-1598.XU Shi-zhe，YU Hai-long，LI Hai-xia，et al.The inversion of apparent density based on the separation and continuation of potential field［J］.Chinese Journal of Geophysics，2009，52(6):1592-1598.

［19］PATELLA D.Introduction to ground surface selfpotential tomography［J］.Geophysics Prospecting，1997，45(4):653-681.

［20］IULIANO T，MAURIELLO P，PATELLA D.Looking inside Mount Vesuvius by potential fields integrated probability tomographies［J］.Journal of Volcanology and Geothermal Research，2002，113(3/4):363-378.

［21］毛立峰，王绪本，高永才.大地电磁概率成像的效果评价［J］.地球物理学报，2005，48(2):429-433.MAO Li-feng，WANG Xu-ben，GAO Yong-cai.Appraisement on the MT probability tomography［J］.Chinese Journal of Geophysics，2005，48(2):429-433.

(编辑 修荣荣)

3D inversion method of density based on double constraint

LIU Zhan1，YU Hui-zhen1，CHEN Ting2

(1.School of Geosciences in China University of Petroleum，Qingdao 266555，China;2.Prospecting Team in Sichuan Department of Geology and Mineral，Chengdu 610072，China)

To increase the density of the restriction of 3D inversion process and avoid"resting on"phenomenon，from a combination of qualitative and quantitative point of view，a dual constraints mechanism was proposed using the combination of the probability density associated with the imaging and field characteristics based drilling bit constraint method.In the mechanism，the mean difference of core function was taken as scanning function，and specific programs for the density's dual constraints based on the relevant difference between subdivision grid density distribution of property and field were given.The theoretical model test results show that if the qualitative distribution characteristics are given，the value range and the search order of decomposing mesh in the process of quantitative inversion can be obtained by using this method，especially in the area of drilling data of seismic data，the inversion vertical resolution is greatly improved.Finally，the actual data processing results further prove that this method is of practicability.

geophysical prospecting;gravity and magnetic prospecting;density inversion;imaging techniques;double constraint;golden section;constraint expansion;occurrence space

P 631.1

A

10.3969/j.issn.1673-5005.2011.06.007

2011－01－10

国家重大专项课题(2008ZX05020－006)

刘展(1957－)，男(汉族)，四川内江人，教授，博士，博士生导师，从事物探资料处理与解释研究工作。

1673-5005(2011)06-0043-08