地球物理学报  2016, Vol. 59 Issue (5): 1849-1860   PDF    
基于协同克里金方法的重力梯度全张量三维约束反演
耿美霞1,2, 黄大年2, 于平2, 杨庆节3    
1. 中国地质大学(武汉)地球物理与空间信息学院, 武汉 430074;
2. 吉林大学地球探测科学与技术学院, 长春 130026;
3. 中国科学院测量与地球物理研究所, 武汉 430077
摘要: 随着重力梯度全张量测量技术的日趋成熟,重力梯度全张量数据的三维反演技术日益受到重视与关注.全张量数据反演与重力数据反演一样仍然面临着严重的多解性问题.本文将基于地质统计学的协同克里金方法应用于重力梯度全张量数据三维反演,建立了密度约束下的多变量协同克里金联合反演方程,以降低反演的多解性.模型试验表明密度信息的加入能够有效降低反演的多解性,提高反演结果的分辨率,尤其是纵向分辨率能够得到显著提高.最后对美国德克萨斯州一个岩丘区所获得实际资料的应用表明了本文方法的实用性.
关键词: 重力梯度全张量     三维约束反演     协同克里金方法     地质统计学    
Three-dimensional constrained inversion of full tensor gradiometer data based on cokriging method
GENG Mei-Xia1,2, HUANG Da-Nian2, YU Ping2, YANG Qing-Jie3    
1. Institute of Geophysics and Geomatics, China University of Geosciences, Wuhan 430074, China;
2. Geo-Exploration Science and Technology Institute, Jilin University, Changchun 130026, China;
3. Chinese Academy of Sciences, Institute of Geodesy and Geophysics, Wuhan 430077, China
Abstract: As full tensor gradiometer (FTG) measurement techniques has increasingly become mature, three-dimensional inversion of FTG data is focused on more and more. Like inversion of gravity, inversion of FTG data confronts with multi-solution problem. In this paper, cokriging inversion equation is framed with several components with known densities as constrains. The method proposed can easily conclude the prior information, for example the known densities, into the inversion equation, so that non-uniqueness of inversion can be reduced. This method is tested on sythetic models and the inverted results show that the resolution, especially the vertical resolution of the resulted model can be significantly improved with the densities as constrains. The application of actual data obtained from a salt dome located in Texas proved the validity of the proposed method.
Key words: Full tensor gradiometer     3-D constrained inversion     Cokriging method     Geostatistics    
1 引言

重力梯度测量能够观测重力位的二阶导数,因而,与传统的重力测量相比,具有更高的分辨率,能够更好地反映浅部异常的边界(曾华霖,1999).在移动平台中(例如船和飞机),梯度仪不像重力仪那样易受共模噪声和大的加速的影响,因此与重力数据相比,重力梯度数据具有更高的信噪比(Martinez et al.,2013).随着重力梯度全张量(FTG)仪器的问世,FTG数据能够提供更加丰富的多分量信息,每个梯度分量包含的信息量不尽相同,综合利用各个张量信息有助于提高地质解释的准确性.

随着重力梯度测量的商业化和计算机技术的快速发展,重力梯度数据的三维反演技术成为国内外学者研究的热点.重磁及其梯度三维反演问题在数学上呈病态.求解病态问题可以得出多个解,这些解都满足数据本身的要求,在地球物理反演中,被称为多解性问题.引起多解性问题主要有两个原因:一方面,实际观测数据是离散且有限的,并伴有来自仪器和周围环境的噪声干扰,另一方面是场源等效性.通过增加野外观测数据和提高观测精度的办法,可以在一定程度上降低反演的多解性,然而针对实际有限信息的观测数据,更有效的方法是在反演中加入先验信息进行约束反演(姚长利等,2002).

多年来,约束条件的使用越来越受到重视,许多地球物理学者提出了各种不同的约束反演方法.Li和Oldenburg(1998)采用Tikhonov正则化方法对重力数据进行三维反演,通过在反演目标函数中引入粗糙度矩阵得到光滑反演模型,该方法也被称为光滑约束反演.Li(2001)将该方法应用于重力梯度数据三维反演.Portniaguine和Zhdanov(1999)提出了聚焦反演方法并将其应用于磁场数据、重力及重力梯度数据反演(Zhdanov et al.,2004),反演结果表明该方法对具有尖锐边界的地质体有较好的反演效果.在国内,姚长利等(2007)针对三维物性反演中存在的多解性和内存消耗巨大的问题,提出了随机子域反演技术.该方法通过概率方式控制子域生成的分布,实现三维反演约束新机制.郭良辉等(2009)提出了基于重力梯度数据异常分离的三维相关成像方法.孟小红等(2012)提出基于剩余异常相关成像的重磁物性反演方法,通过对物性模型的正演和实测结果的残差进行相关成像并迭代更新物性模型,实现了相关成像方法对物性参数的反演.

约束条件的使用在地球物理反演中至关重要,但约束条件的提取以及同反演计算的结合仍然存在很多困难,例如如何将一些地质信息转化成具体数学形式的约束条件还存在着困难,另外,有些约束条件还难以同反演计算过程融合,从而制约了约束条件作用的发挥(姚长利等,2002Li,2012).因此,发展新的容易与先验信息相结合的反演方法将具有重要实用价值(Sun and Li,2000).近年来,地质统计学方法被应用到地球物理反演中.地质统计学方法在地球物理反演中应用的好处,除了对区域地质特征规律信息的应用之外,还体现在能够灵活地结合地质学家对地质情况的认识上.Asli等(2000)建立了以待反演密度作为主变量,重力数据作为次级变量的两变量协同克里金反演方程;Gloaguen等(20052007)将协同克里金反演方法应用于井中雷达速度反演和层析成像领域;Shamsipour等(2010)采用协同克里金方法对重力数据进行反演;Geng等(2014)采用协同克里金方法对重力梯度全张量数据进行反演,反演结果表明重力梯度全张量数据联合反演有助于提高反演结果的准确性与分辨率,但反演结果在纵向上仍然存在一定程度的发散.

为了进一步提高协同克里金反演方法结果的分辨率,提高该方法与先验信息有效融合的能力,本文在Geng等人(2014)工作基础上,提出密度约束下的协同克里金反演方法.为检验方法的有效性,研究密度信息的加入对反演结果的影响,对直立长方体组合模型和倾斜岩脉组合模型的重力梯度全张量数据进行了反演计算.最后将本文研究方法应用于美国德克萨斯州的一个研究构造上所获得的数据,证明了本文方法的有效性和应用前景.

2 重力梯度反演的基本原理2.1 重力梯度张量与正演

在笛卡尔坐标系下,对于地下任意形状的地质体已知其密度分布,其重力梯度张量,即重力场三分量沿3个坐标系方向的梯度可表示为

由于在地球外部,引力位满足拉普拉斯方程Txx+Tyy+Tzz=0,以及重力场的无旋性:Txy=Tyx,Txz=Tzx,Tyz=Tzy,因此,上式(1)中重力梯度张量的9个分量中只有5个是独立的.将待反演的地下空间离散成一系列大小相等紧密排列的直立长方体单元,并赋予每个单元固定的密度值,地表重力梯度张量各分量观测值与地下密度分布可以表示为如下线性关:

其中m为网格单元的剩余密度,d为观测点上的重力梯度张量分量的理论观测值,Gαβ(α,β=x,y,z)为该张量的正演核函数,它是一个离散单元相对观测点位置所决定的常数矩阵(郭志宏等,2004).2.2 协同克里金反演

协同克里金方法(Myers,1982;Chilès and Delfiner,2012)是一种基于地质统计学的估计方法,它通过利用主变量和次级变量的空间相关性使主变量的估计方差达到最小,从而对主变量进行估计.在重力梯度张量数据反演中,主变量为密度值ρ,其估计值为ρ*,次级变量为重力梯度张量数据Tαβ(α,β=x,y,z).假设密度和重力梯度数据的期望值具有空间一致性,且它们的期望值都为0.我们可以由下式得到主变量ρ的估计方差(Myers,1982Shamsipour et al.,2010):

其中,Λ是加权系数矩阵,Cρρ是密度协方差矩阵,CTαβ,ρ是重力梯度数据和密度的互协方差矩阵,CTαβ,Tαβ是重力梯度数据的协方差矩阵(Asli et al.,2000).通过使主变量ρ估计方差最小,可以得到单分量反演的协同克里金方程:

根据式(4),我们可以得到最优加权系数矩阵Λ,然后利用最优加权系数矩阵,我们可以得到密度估计值:

协同克里金方差向量可由下式得到:

估计方差位于式(6)的主对角线上,非对角元素是估计误差的互协方差.由式(2),我们知道密度和重力梯度数据间是线性关系,因此密度和重力梯度数据互协方差矩阵也是线性相关的(Myers,1982):

其中Q0是重力梯度观测误差协方差矩阵.观测数据经过一系列校正和处理后,可以认为数据中的噪声是不相关的,则矩阵Q0=σ2I是观测数据的标准偏差,I为单位矩阵.

重力梯度数据和密度的互协方差与密度协方差之间存在如下的关系(Asli et al.,2000):

不考虑重力梯度数据中的观测误差,即Q0=0,可以证明下式成立:
式(9)证明了采用协同克里金方法反演得到的模型的正演数据与观测异常在不考虑噪声影响的情况下是相等,从而保证了反演结果的正演数据与实际观测数据能够在噪声水平内得到较好的拟合.

如果已知先验信息,例如岩石的密度信息,用ρF表示,则建立密度数据约束下的重力梯度全张量数据反演的协同克里金方程如下式所示:

其中ΛΓΩγHΦ分别为分量TxzTyzTzz TxxTxyρF对应的最优加权矩阵.左侧矩阵中对角元素为次级变量的协方差矩阵,非对角元素为不同次级变量之间的互协方差矩阵.对式(10)求解可以得到最优加权矩阵,从而可以得到6个次级变量联合反演的密度估计值:
协同克里金估计方差可由下式计算得到:
其中 .估计方差位于矩阵σck的主对角线上.

由式(11)可知,采用协同克里金方法反演密度必须先计算出密度协方差矩阵Cρρ.要想得到Cρρ,需要先计算变差函数参数,即块金值C0x、y、z三方向的变程ax、ay、az,以及基台值C+C0(C称为拱高).然后由计算得到的变差函数参数计算出密度协方差矩阵Cρρ(孙洪泉,1990).估计变差函数参数常用的方法有三种(Chasseriau and Chouteau,2003Shamsipour et al.,2010):(1)根据已有的地质信息,建立初步估计模型,然后利用地质统计学方法得到变差函数参数(附录A);(2)当来自地下密度数据(例如测井密度)充足时,可以直接根据密度数据通过地质统计学方法计算变差函数参数(附录A);(3)在先验地质信息缺乏的情况下,可采用Alis等(2000)提出的V-V绘制方法估计变差函数参数(附录B).

2.3 积分灵敏度矩阵

在位场反演中核函数随积分单元与观测位置之间的距离增大而迅速衰减的性质会导致反演结果的物性参数分布全部集中于地表,这种现象被称为反演的“趋肤效应”.为了克服这种效应,Li和Oldenburg(1996,1998)在三维磁化率和密度反演中引入深度加权函数w(z)=(z+z0)-β.但该深度加权函数仅考虑深度方向的情况,未能考虑反演的横向分辨率.此外,反演结果受参数β影响很大.因此本文在协同克里金方程中引入积分灵敏度矩阵.由于积分灵敏度矩阵能够对x,y,z三个方向加权,使得反演结果具有三个方向的分辨能力,从而使反演算法的横向分辨能力得到提高.

观测数据与模型变化的关系可表示为

其中,Gik为重力梯度数据的核函数.数据灵敏度对模型mk的积分可以表示为
可以看到,数据灵敏度只与核函数有关,不涉及其他参数,从而避免了因参数选择不同对反演结果的影响.在离散化的协同克里金反演方程中,积分灵敏度矩阵是一个对角矩阵:
向式(5)中引入积分灵敏度矩阵后,单分量协同克里金反演方程变为
3 理论模型试验

为了检验所建立的反演方程的正确性,并研究密度信息的加入对反演结果的影响,设计了两个理论模型进行反演试验.理论模型试验的具体过程为:首先根据设计的理论模型采用快速正演算法(姚长利等,2003)计算出重力梯度各个分量的正演数据,并向各个分量中加入5%的高斯白噪声;然后对重力梯度张量的五个独立分量进行联合反演,最后在已知的密度数据约束下对五个独立分量进行联合反演.

3.1 双直立长方体模型

首先设计由两个直立长方体组成的组合模型一.设置地下模型空间尺寸x,y,z为2400 m×2400 m×1500 m.将地下模型空间剖分为24×24×15个紧密排列的直立立方体单元,每个单元为100 m×100 m×100m.地面观测网格为24×24的网格,观测点位于每个网格中心.设地下模型空间中有两个400 m×400 m×500 m的长方体模型,剩余密度分别为-1 g·cm-3和1 g·cm-3.图 1a显示了合成模型的密度分布三维图.反演计算前,向该理论模型在地表产生的重力梯度张量的6个分量各加入最大幅值5%的高斯白噪声.加入噪声干扰后的Tzz分量如图 1b所示.

图 1 (a) 组合长方体理论模型; (b) 理论模型Tzz分量正演结果(含5%的高斯随机噪声) Fig. 1 (a) The synthesized prism model; (b) calculated Tzz component with 5% uncorrelated Gaussian noise of maximum of the datum magnitude

模型试验中,我们采用上述第一种方法计算变差函数参数值.采用高斯模型对估计模型(这里为真实模型)拟合结果如图 2(ac)所示,得到变差函数参数:C0=2000(kg·m-3)2C=24000(kg·m-3)2ax=450 m,ay=450 m,az=450 m.

图 2 高斯模型拟合变差函数
(a) x方向拟合结果; (b) y方向拟合结果; (c) z方向拟合结果.
Fig. 2 Fitting variogram function using Gaussian model
(a) x-direction variogram; (b) y-direction variogram; (c) z-direction variogram.

首先对张量的5个独立分量TxyTxzTyyTyzTzz进行联合协同克里金反演,反演结果在z=600 m横向切片和y=1300 m纵向切片如图 3a3b所示,其中黑色的方框表示真实模型的边界.可以看到,反演结果较准确地反映出真实异常体的空间位置和形态,但是剩余密度值的绝对值比真实值小,并且异常体底部密度值存在一定程度的发散.然后加入两口测井中的密度数据,采用新建立的式(11)对五个独立分量数据进行联合反演,其中两口测井的水平中心坐标分别为(0.8 km,1.3 km)和(1.8 km,1.3 km),即测井恰好垂直穿过两个异常体的中心,地表下方测井所穿过的网格单元的剩余密度为已知,模拟多井资料.图 3c3d分别显示了在密度约束下五个独立分量联合反演结果的横向和纵向切片.可以看到与无密度约束反演结果相比,加入密度信息约束反演得到的异常体边界更加清晰,尤其是纵向分辨率得到显著提高;此外,还可以看到测井穿过的单元附近的剖分单元也能够比较准确反演.加入密度数据前后反演结果与真实模型的均方根误差分别为0.10 g·cm-3和0.07 g·cm-3,从而从定量角度说明了加入密度数据进行约束反演能够提高反演结果的精度.

图 3 无密度约束下张量的五个独立分量联合反演结果:(a) z=600 m横向切片; (b) y=1200 m纵向切片.密度约束下张量的五个独立分量联合反演结果图: (c) z=600 m横向切片; (d) y=1200 m纵向切片 Fig. 3 The recovered model with five independent componentswithout constraints: (a) In the z=600 m depth section; (b) In the y=1200 m cross section. The recovered model with five independent components with constraints: (c) In the z=600 m depth section; (d) In the y=1200 m cross section
3.2 脉状体组合模型

下面设计一个较为复杂的倾斜脉状体组合模型,如图 4a所示,两个脉状体的倾角都为45°,剩余密度分别为0.8 g·cm-3和1 g·cm-3.我们知道,位场数据的纵向分辨率非常差,纵向上存在叠加的异常体很难被区分.这一组合模型用来研究密度约束在反演复杂模型中的作用.设置地下模型空间尺寸x,y,z为2400 m× 2400 m ×1500 m.将地下模型空间剖分为24×24×15个紧密排列的直立立方体单元,每个单元为100 m×100 m×100 m.地面观测网格为24×24的网格,观测点位于每个网格中心.反演计算前,向该理论模型在地表产生的重力梯度张量的六个分量各加入5%的高斯白噪声,带有噪声干扰的Tzz分量如图 4b所示.

图 4 (a) 组合倾斜脉状体理论模型; (b)理论模型Tzz分量正演结果(含5%的高斯随机噪声) Fig. 4 (a) The synthesized dike model; (b) Calculated Tzz component with 5% uncorrelated Gaussian noise of the datum magnitude

采用高斯模型对估计模型(这里为真实模型)进行拟合.拟合结果如图 5(ac)所示,得到变差函数参数:C0=2000(kg·m-3)2C=24000(kg·m-3)2ay=250 m,a45°=1000 m,a135°=380 m,其中a45°的方向为平行于脉状体倾向的方向,a135°的方向垂直于y轴和脉状体倾向方向所构成的平面.首先对5个独立分量TxyTxzTyyTyzTzz进行联合协同克里金反演,反演结果在z=600 m横向切片和y=1300 m纵向切片如图 6a6b所示,可以看到纵向上存在叠加的两个倾斜脉状体都得到了比较准确的恢复:两个脉状体的位置、倾角都与真实模型基本一致,但是脉状体底部密度值发散比较严重.然后加入一口测井中的密度数据作为约束条件再对五个独立分量数据进行联合反演,测井的水平中心坐标为(1.3 km,1.3 km),地表下方测井所穿过的网格剩余密度已知,模拟单口井资料.图 6c6d分别显示了五个独立分量在密度约束下联合反演结果的横向和纵向切片图,可以看到与无密度约束反演结果相比,加入密度信息约束后,反演得到的两个倾斜脉状体的边缘发散现象得到明显改善,脉状体边缘更加清晰,并且两个脉状体的剩余密度与各自的真实值也非常接近.加入密度数据前后反演结果与真实模型的均方根误差分别为0.09 g·cm-3和0.07 g·cm-3,再次证明了密度数据的加入对提高重力梯度全张量数据反演精度的能力.

图 5 高斯模型拟合变差函数
(a) y方向拟合结果; (b) 平行于脉状体倾向方向拟合结果; (c) 垂直于脉状体倾向方向拟合结果.
Fig. 5 Fitting variogram function using Gaussian model
(a) y-direction variogram; (b) Variogram along the inclination of dikes; (c) Variogram perpendicular to the inclination of dikes.

图 6 无密度约束下张量的五个独立分量联合反演结果:(a) z=600 m横向切片; (b) y=1200 m纵向切片.密度约束下张量的五个独立分量联合反演结果: (c) z=600 m横向切片; (d) y=1200 m纵向切片 Fig. 6 The recovered model with five independent componentswithout constraints: (a) The z=600 m depth section; (b) The y=1200 m cross section. The recovered model with five independent components with constraints: (c) The z=600 m depth section; (d) The y=1200 m cross section
4 实测数据应用

将协同克里金反演方法应用于美国德克萨斯州的一个研究构造上所获得的重力数据(Nettleton,1976)的解释中.利用FFT变换(Michus and Hinojosa,2001)计算得到异常的全张量数据.研究区内为沉积岩相地下发育盐丘构造.石油和天然气的沉积及储层与盐丘有密切关系,盐丘区附近常常拥有丰富的油气储量,因此对盐丘地区的研究具有重要价值(Zhou,2006Ennen and Hall,2011梁杰等,2010).由于盐丘下部构造往往非常复杂,地震波在盐丘下部迅速衰减,地震照明度低,其地震成像面临着非常大的挑战(Liu et al.,2011).在盐丘地区,对重力或重力梯度全张量数据进行三维反演,得到的密度模型可为该地区的研究工作提供依据;此外,得到的密度模型可以通过适当的密度-速度函数转化成速度模型,从而实现地震-重力资料的综合解释,提高解释结果的可靠性.

重力梯度反演选择重力梯度异常低值区域13 km×13 km,并将数据网格化,网格间距是500 m,网格化后的重力梯度张量的各个分量等值线如图 7所示.将反演目标区域地下空间划分为26×26×20个紧密排列的直立正方体单元,每个单元的尺寸为500 m×500 m×500 m.采用V-V绘制方法得到变差函数参数值:C0=2000(kg·m-3)2C=26000(kg·m-3)2a-40°=5000 m(岩丘延伸方向),a50°=4000 m(垂直岩丘延伸方向),az=3000 m.

图 7 重力梯度全张量异常
(a) Txx; (b) Txy; (c) Tzx; (d) Tyy; (e) Tyz; (f) Tzz.
Fig. 7 Gravity gradient tensor data

首先对五个独立分量TxxTxyTxzTyzTzz进行联合协同克里金反演,反演结果剩余密度三维分布如8a所示,反演结果在z=4500 m的横向切片如图 8b所示,岩丘延伸方向的纵向切片(图 8b上虚线l1所示的位置)和垂直岩丘延伸方向的纵向切片(图 8b虚线l2所示的位置)分别如图 8c8d所示.为了检验密度数据约束反演在提高反演结果分辨率方面的能力,设置了五口测井,已有地质调查和钻井数据显示(Said,1962)该岩丘的剩余密度在-0.3 g·cm-3左右,因此令这五口测井穿越过的岩石的剩余密度为-0.3 g·cm-3.在五口井中密度数据约束下,对重力梯度全张量数据进行协同克里金反演,反演结果剩余密度三维分布如图 8e所示,相对应的切片结果如图 8(fgh)所示(其中白色的线和点表示测井所在位置和穿越的单元格).对比图 8(ad)和8(eh)可以明显看到,在密度数据的约束下,反演结果的分辨率得到明显提高,测井穿过的单元附近的剖分单元剩余密度也非常接近-0.3 g·cm-3.从盐丘区的反演结果看,低密度值分布约在地下3~6.5 km的深度;岩丘中心埋深约为4.5 km,这一结果与以往研究结果(见表 1)相一致.本文采用的5组密度数据是根据已有调查和研究结果所模拟的数据,在实际数据应用中,在缺少井数据的情况下,可以将地质调查和/或其他物探方法得到的密度信息作为约束进行联合反演,同样可以起到降低反演多解性,提高反演结果的分辨率的作用.

图 8 重力梯度张量的五个独立分量联合反演结果: (a) 剩余密度分布三维图; (b) z=4500横向切片; (c) 岩丘走向方向纵向切片; (d) 垂直于走向方向纵向切片. 密度约束下五个独立分量联合反演结果: (e) 剩余密度分布三维图; (f) z=4500 m横向切片; (g) 岩丘走向方向纵向切片; (h) 垂直于走向方向纵向切片 Fig. 8 The recovered model with five independent components: (a) Recovered 3-D model; (b) The depth slice of 4500 m; (c) Cross section along the strike direction of salt dome; (d) Cross section perpendicular to the strike direction of salt dome. The recovered model with five independent components with constraints: (e) Recovered 3-D model; (f) The depth slice of 4500 m; (g) Cross section along the strike direction of salt dome; (h) Cross section perpendicular to the strike direction of salt dome

表 1 密度约束下协同克里金反演结果与先前计算岩丘中心埋深结果比较 Table 1 The comparison of the inversion result computed by the constrained cokriging method and previous results
5 结论

本文建立了密度约束下的重力梯度全张量数据联合协同克里金反演方程,使协同克里金反演方法不仅可以将地质学家对地质情况的认识,例如地质体走向、倾角、空间尺度等结合到反演方程,还能够将来自于测井或其他方法获得的密度数据结合到反演方程中,从而大大降低反演的多解性.理论模型试验表明加入密度数据进行约束反演能够显著提高反演结果的分辨率与精度.最后将其应用于美国德克萨斯州盐丘的重力梯度数据的解释,获得结果与前人解释结果相一致,从而证明了该方法的实用性.

致谢 特别感谢提供实际资料的同仁以及为评审本文所付出努力的专家们.附录A 变差函数的估计与拟合

在地质统计学中,变差函数是最基本也是最重要的模拟工具,它用于描述数据值的空间互相关性,将数据点的相关性随空间距离的增大而减小的现象通过数学函数模拟出来.变差函数可以用四个参数来描述,即变差函数类型、变程a、块金效应C0和基台值C+C0.

(1)变差函数的计算公式与估算

变差函数的定义是:区域化变量Z(x)和Z(x+h)两点之差的方差之半.Z(x)的变差函数数学定义如下(孙洪泉,1990):

其中Var(·)表示方差,h为滞后距.如果有了区域变量Z(x)的一部分采样样本,就可以估算该区域变化量Z(x)的变差函数,具体计算公式如下:
其中i为样本序号.当来自地下(例如通过测井方法得到)密度数据量充足时,变差函数可以直接利用已知的密度数据根据式(A2)计算得到;或解释人员根据已有的地质信息,将地下解释空间剖分成紧密排列的直立棱柱体单元,建立初步解释模型,然后再根据式(A2)计算变差函数.(2)变差函数的拟合

下面我们给出根据初步解释模型估算变差函数的示例.

设置地下模型空间尺寸x,y,z为2100 m× 2100 m ×1500 m.将地下模型空间剖分为21×21×15个紧密排列的直立立方体单元,每个单元尺寸为100 m×100 m×100 m.设在地下模型空间中有一个300 m×500 m×700 m的长方体模型,剩余密度为1 g·cm-3,如图A1所示.

图 A1 直立长方体模型 Fig. A1 Prism model

由于这一地质模型属性具有各向异性特征,应当计算它在不同方向的变差函数值.滞后距h是指向某方向移动的距离,计算时,在指定的方向上,对指定的h,搜索所有相距h的点对[Z(xi),[Z(xi+h)],并统计点对个数N,根据式(A2)计算出滞后距h对应的变差值.然后再计算出滞后距2h、3h、4h…对应的变差值.图A2(ac)中的小圆圈分别给出了x,y,z三个方向的变差的计算结果.在统计模拟中,变差函数采用的是一个解析函数,因此,计算出不同滞后距h的变差值V(h)后,还需要拟合出一个解析函数.变差函数的拟合有很多种方法,比如在指定出函数类型后,可以使用加权最小二乘法、线性规划法、遗传算法等.其中遗传算法能够拟合任何指定的函数类型,因此,本文采用遗传算法进行拟合.

图 A2 高斯模型拟合变差函数 Fig. A2 Fitting variogram function using Gaussian model

常用的解析函数类型有三种,分别是球状模型、高斯模型和指数模型.选用哪种函数拟合,一方面取决于样本数据,另一方面,也取决于对实际规律的掌握程度.这里我们采用高斯模型,它的形式是

变程a高斯=.图A2(ac)黑色实线给出了采用高斯模型拟合结果.可以得到:x方向上变程ax=330 m,块金值C0=500(kg·m-3)2,基台值C+C0=10760(kg·m-3)2y方向上变程ay=440 m,块金值C0=500(kg·m-3)2,基台值C+C0=11000(kg·m-3)2z方向上变程az=520 m,块金C0=320(kg·m-3)2,基台值C+C0=11000(kg·m-3)2.三个方向上的块金值求平均得到综合块金值C0=440(kg·m-3)2,基台值求平均得到综合基台值C+C0=10920(kg·m-3)2.变差函数的各项参数包含了模型的各向异性特征、几何属性,以及物性属性信息,也就是说先验地质信息以具体的参数形式包含在了变差函数中.由变差函数计算协方差矩阵的方法可以参考孙洪泉的《地质统计学及其应用》,这里不再赘述.附录B V-V绘制方法

V-V绘制方法实现步骤:

(1)选择一个变差模型,包括解析函数类型和三个参数:基台值、变程和块金值.由此,计算出密度协方差矩阵Cρρ.

(2)根据gαβ·gαβT计算出实验协方差矩阵,根据式(7)计算出理论协方差矩阵.

(3)对理论协方差从小到大进行排序并分组,用同样的排列顺序对实验协方差进行排序并分组.

(4)计算分完组后的理论协方差和实验协方差的平均绝对误差.

(5)修改变差模型参数使理论协方差和实验协方差的平均绝对误差最小,从而使二者之间具有较高的相关性.拟合过程中,结合已有的地质信息修改变差模型参数,有利于得到好的拟合结果.

参考文献
[1] Abdelrahman E M, Bayoumi A I, Elaraby H M.Short note: A least-squares minimization approach to invert gravity data. Geophysics, 1991, 56: 115-118.
[2] Asli M, Marcotte D, Chouteau M. 2000. Direct inversion of gravitydata by cokriging.//Kleingeld W,Krige D eds.Geostats 2000, Proceedings of the 6th International Geostatistics Congress, 64-73.
[3] ChasseriauP, Chouteau M. 2003. 3D gravity inversion using a modelof parameter covariance.Journal of Applied Geophysics, 52(1): 59-74.
[4] Chilès J P, Delfiner P. Geostatistics: Modeling Spatial Uncertainty, Second Edition, Hoboken, NJ, Wiley, 2012.
[5] Ennen C, Hall S. 2011. Structural mapping of the Vinton salt dome, Louisiana, using gravity gradiometry data.// 2011SEGAnnual Meeting Technical Program. Society of Exploration Geophysicists.Expanded Abstracts, 830-835.
[6] Essa K S.Gravity data interpretation using the s-curves method.J.Geophys.Eng., 2007, 4(2): 204-213.
[7] Geng M X, Huang D N, Yang Q J, et al. 2014. 3D inversion of airborne gravity-gradiometry data using cokriging.Geophysics, 79(4): G37-G47.
[8] Gloaguen E, Marcotte D, Chouteau M, et al. 2005. Borehole radar velocity inversion using cokriging and cosimulation.Journal of Applied Geophysics, 57(4): 242-259.
[9] Gloaguen E, Marcotte D, Giroux B, et al. 2007. Stochastic borehole radar velocity and attenuation tomographies using cokriging and cosimulation.Journal of Applied Geophysics, 62(2): 141-157.
[10] Guo L H, Meng X H, Shi L, et al. 2009. 3-D correlation imaging for gravity and gravity gradiometry data.Chinese J.Geophys. (in Chinese), 52(4): 1098-1106, doi:10.3969/j.issn.0001-5733.2009.04.027 .
[11] Guo Z H, Guan Z N, Xiong S Q. 2004. Cuboid ΔT and its gradient forward theoretical expressions without analytic odd points. Chinese J.Geophys.(in Chinese), 47(6): 1131-1138, doi:10.3321/j.issn:0001-5733.2004.06.029.
[12] Li Y G, Oldenburg D W. 1998. 3-D inversion of gravity data.Geophysics, 63(1): 109-119.
[13] Li Y G. 2001. 3-D inversion of gravity gradiometer data.// 2001 SEG Annual Meeting Technical Program.Society of Exploration Geophysicists.Expanded Abstracts,1470-1473.
[14] Li Y G. 2012. Recent advances in 3D generalized inversion of potential-field data. // 2001 SEG Annual MeetingTechnical Program.Society of Exploration Geophysicists.Expanded Abstracts, 1-7.
[15] Liang J, Gong J M, Cheng H Y. 2010. Control of salt rock distribution on oil and gas pooling in the Gulf of Mexico.Marine Geology Letters (in Chinese), 26(1): 25-30.
[16] Liu YK, Chang X, Jin D G, et al. 2011. Reverse time migration of multiples for subsalt imaging. Geophysics, 76(5): WB209-WB216.
[17] Ma G Q, Du X J, Li L L. 2012. Comparison of the tensor local wavenumber method with the conventional local wavenumber method for interpretation of total tensor data of potential fields.Chinese J.Geophys. (in Chinese), 55(7): 2450-2461, doi: 10.6038/j.issn.0001-5733.2012.07.029.
[18] Martinez C, Li Y G, Krahenbuhl R,et al. 2013.3D inversion of airborne gravity gradiometry data in mineral exploration: A case study in the Quadrilátero Ferrífero, Brazil.Geophysics, 78(1): B1-B11.
[19] Meng X H, Liu G F, Chen Z X, et al. 2012. 3-D gravity and magnetic inversion for physical properties based on residual anomaly correlation. ChineseJ.Geophys. (in Chinese), 55(1): 304-309, doi: 10.6038/j.issn.0001-5733.2012.01.030.
[20] Mickus K L, Hinojosa J H. 2001. The complete gravity gradient tensor derived from the vertical component of gravity: A Fourier transform technique. Journal of Applied Geophysics, 46(3): 159-174.
[21] Mohan N L, Anandababu L, Roa S.Gravity interpretation using Mellin transform.Geophysics, 1986, 52(1): 114-122.
[22] Myers D E. 1982.Matrix formulation of co-kriging.Mathematical Geology, 14(3): 249-257.
[23] Nettleton L L. 1976. Gravity and Magnetic in Oil Prospecting. New York: McGraw-Hill Book Co.
[24] Qruc B.Depth estimation of simple causative sources from gravity gradient tensor invariants and vertical component. Pure and Applied Geophysics, 2010, 167(10):1259-1272.
[25] Portniaguine O, Zhdanov M S. 1999.Focusing geophysical inversion images.Geophysics, 64(3): 874-887.
[26] Said R. 1962.The Geology of Egypt.Amsterdam: Elsevier.
[27] Salem A, Ravatz D, Mushayandebvu M F, et al.Linearized least-squares method for interpretation of potential-field data from sources of simple geometry.Geophysics, 2004, 69(3): 783-788.
[28] Shamsipour P, Marcotte D, Chouteau M, et al. 2010. 3D stochastic inversion of gravity data using cokriging and cosimulation.Geophysics, 75(1):I1-I10.
[29] Shaw R K, Agarwal B N P.A generalized concept of resultant gradient to interpret potential field maps.Geophysical Prospecting, 1997, 45(3): 513-520.
[30] Sun H Q. 1990. Geological Statistics and Application. Beijing: China University of Mining and Technology Press.
[31] Sun J J, Li Y G. 2010.Inversion of surface and borehole gravity with thresholding and density constraints.// 2010 SEG Annual Meeting Technical Program.Society of Exploration Geophysicists.Expanded Abstracts, 1798-1803.
[32] Yao C L, Hao T Y. Guan Z N. 2002.Restrictions in gravity and magnetic inversions and technical strategy of 3D properties inversion. Geophysical& Geochemical Exploration (in Chinese), 26(4): 253-257.
[33] Yao C L, Hao T Y, Guan Z N. 2003.High-speed computation and efficient storage in 3-D gravity and magnetic inversion based on genetic algorithms.Chinese J.Geophys. (in Chinese), 46(2): 252-258, doi:10.3321/j.issn:0001-5733.2003.02.020.
[34] Yao C L, Zheng Y M, Zhang Y W. 2007. 3-D gravity and magnetic inversion for physical properties using stochastic subspaces.Chinese J.Geophys. (in Chinese), 50(5): 1576-1583, doi:10.3321/j.issn:0001-5733.2007.05.035.
[35] Zeng H L. 1999. Present state and revival of gravity gradiometry. Geophysical& Geochemical Exploration (in Chinese), 23(1): 1-6.
[36] Zhdanov M S,Ellis R, Mukherjee S.2004. Three-dimensional regularized focusing inversion of gravity gradient tensor component data. Geophysics, 69(4): 925-937.
[37] Zhou H W. 2006. First-break vertical seismic profiling tomography for Vinton Salt Dome.Geophysics, 71(3): U29-U36.
[38] 郭良辉, 孟晓红, 石磊等. 2009. 重力和重力梯度数据三维相关成像. 地球物理学报, 52(4): 1098-1106, doi:10.3969/j.issn.0001-5733.2009.04.027.
[39] 郭志宏, 管志宁, 熊盛青. 2004. 长方体ΔT场及其梯度场无解析奇点理论表达式. 地球物理学报, 47(6): 1131-1138, doi:10.3321/j.issn:0001-5733.2004.06.029.
[40] 梁杰, 龚建明, 成海燕. 2010. 墨西哥湾盐岩分布对油气成藏的控制作用. 海洋地质动态, 26(1): 25-30.
[41] 马国庆, 杜晓娟, 李丽丽. 2012. 解释位场全张量数据的张量局部波数法及其与常规局部波数法的比较. 地球物理学报, 55(7): 2450-2461, doi: 10.6038/j.issn.0001-5733.2012.07.029.
[42] 孟小红, 刘国峰, 陈召曦等. 2012. 基于剩余异常相关成像的重磁物性反演方法. 地球物理学报, 55(1): 304-309, doi: 10.6038/j.issn.0001-5733.2012.01.030.
[43] 孙洪泉. 1990. 地质统计学及其应用.北京: 中国矿业大学出版社.
[44] 姚长利, 郝天珧, 管志宁. 2002. 重磁反演约束条件及三维物性反演技术策略. 物探与化探, 26(4): 253-257.
[45] 姚长利, 郝天珧, 管志宁等. 2003. 重磁遗传算法三维反演中高速计算及有效存储方法技术. 地球物理学报, 46(2): 252-258, doi:10.3321/j.issn:0001-5733.2003.02.020.
[46] 姚长利, 郑满元, 张聿文. 2007. 重磁异常三维物性反演随机子域法方法技术. 地球物理学报, 50(5): 1576-1583, doi:10.3321/j.issn:0001-5733.2007.05.035.
[47] 曾华霖. 1999. 重力梯度测量的现状及复兴. 物探与化探, 23(1): 1-6.