地表以及航空重力测量目前广泛应用于地质构造调查、能源勘探、工程测量等诸多领域，重力反演是通过地表数据获知地下密度分布信息的重要手段，特别是在上述应用的精细分析阶段，重力反演更是不可或缺.

在过去的几十年中，重力反演发展了一系列的方法，根据反演结果的不同，可以将其归纳为两类：第一类是直接应用线性或者非线性算法反演地下介质的密度异常分布，线性反演是通过最小化目标函数、求解方程，来直接反演地质体密度异常的分布(Li and Oldenburg， 1996，1998; 姚长利等， 2003，2007;冯娟等，2014;李淑玲等，2014)，非线性反演则包括应用蒙特卡洛算法反演重力异常(Bosch et al., 2006)，在重力反演中应用协方差作为先验信(Chasseriau and Chouteau， 2003)，以及地质统计学方法等(Shamsipour et al., 2010a，2010b，2011).上述方法在重力数据的反演计算中都有广泛的应用.

另外一类重力反演方法称为概率层析方法，该方法不是直接反演密度异常，而是通过+1和-1值，称为等效密度，来表征表示地下密度异常的分布.该方法最早应用于自然电位数据的反演中，目前在电法、地磁、重力等数据反演中都得到了推广应用 (Patella， 1997a，1997b; Mauriello et al., 1998; Mauriello and Patella， 1999a，1999b，2001，2005，2008; Bosh and McGaughey， 2001; Boulanger and Chouteau， 2001; Guo et al., 2011a，2011b; 郭良辉等，2009; 孟小红等，2012).

在重力数据的反演中，反演的结果存在非唯一性，即不同的地下密度异常分布可产生相同的地表测量数据(Bosh and McGaughey， 2001; Boulanger and Chouteau， 2001；祁光等，2012；邓震等，2012；王万银等，2013；陈晓红等，2013)，因此，如果使重力数据反演的地下信息有真正的地质意义，则需要在反演过程中引入必要的约束.对比上述的第一和第二种方法，第二种方法虽然方法简单，稳定，但是反演的结果是等效密度值，很难将已知构造或者岩性信息与之建立联系.另外，在实际的反演过程中，我们更希望得到的是真正的密度分布.基于此，本文提出了一种由概率层析成像结果到密度异常分布的迭代计算方法，同时为了提高成像结果的聚焦效果，在反演过程中加入了密度体上下界限制的约束.通过两个简单模型对本文方法进行了测试计算，证明了方法的有效性.

1 概率层析反演基础

在反演计算中，通常是把反演空间剖分成若干长方体网格，并假定每个网格内的重力值恒定.在正反演计算中，这些长方体成为基本计算单元.在一个坐标系中，水平方向为坐标系的(x－y)面，z轴正方向向下，重力数据在的z=0表面上进行测量.对于地下的某一长方体q，密度异常值为Δρ_q，其地表测量的布格重力异常值Δg_q(x，y，z)为(如图 1a)：

假如地下剖分的长方体数量为Q，则地表测量的布格重力异常值是成像空间剖分后所有长方体所产生重力值的和(如图 1b):

图 1

图 1(a)剖分长方体密度为ρ1背景密度为ρ2； (b)成像区间被剖分成若干长方体.Fig. 1(a)rectangular density is ρ1 and background density is ρ2； (b)imaging space is divided to rectangular.

如果长方体为单位质量，则此时地表测量的布格重力异常值为

2 由概率层析结果到密度异常的转化

正如上文所述，概率层析方法得到的是等效密度，它用介于+1到-1之间的值来表征密度异常相对于背景的盈余或者亏损，该方法简单，稳定，计算速度快，但是，在实际的重力数据反演计算中，更多时候希望能得到真实的密度分布，基于此，本文提出了一种迭代计算流程，用以实现等效密度到真实密度物性的转化.

假设某模型地质背景密度为ρ，对成像区剖分后的某个长方体，其密度为ρ+Δρ，对该模型进行重力正演计算得重力数据data1，对改数据进行概率成像，这个长方体的概率密度值为C，如果其为正，说明Δρ为正，反之为负.以小的密度值乘以C，此时长方体的密度值ρ+Δρ2，对此模型进行正演计算，得到重力数据data2.此时，对data1－data2的差值进行概率成像，结果为C1，而实际上，其为Δρ2－Δρ的概率成像结果.乘以一个小的密度值Δρ3，可以将C1转化为密度，它便可以理解为反演过程中的密度扰动，而Δρ2+C1*Δρ3便为反演计算更新后的密度模型.

概括上述方法，概率成像结果到密度之间的转化可以表达为如下流程：

(1)采集重力数据data1，并对其进行概率层析计算，得成像结果为C.

(2)将C乘以小密度值得到密度模型m1，并将其作为密度反演的初始值.

(3)对模型m1进行正演计算，得到重力数据data2.

(4)对data1－data2进行概率层析，其结果为C1.

(5)用m1+C1*Δρ更新密度模型.

(6)重复(3)-(5)，直到正演的数据与采集的数据的相对误差小于一个设定的门限值.

3 模型测试

首先采用一个2D模型对本文方法进行测试，该模型背景异常为0，一个密度为0.5 g/cm³镶嵌于模型中上部，在正反演计算中，将模型剖分成51×51=2601个正方形，正方形的大小为50×50 m，首先对该模型进行概率层析，模型及概率层析的结果如图 2.

图 2

图 2 2D简单模型(a)及其概率成像结果(b). Fig. 2 2D simple model(a) and probability imaging result.

通过结果可以看出，在异常体位置，概率层析结果出现了正的概率成像值，基本能够表征异常体的存在，但是，异常范围偏大，成像聚焦效果较差.应用上述转化方法，我们将概率值转化为密度值，其结果如图 3.

图 3

图 3 本文方法反演结果Fig. 3 inversion result by introduced method

对比两个结果我们可以看出，本文方法反演的结果首先将概率值转化为了真正的密度值，同时改善了异常体的聚焦效果，更好的实现了异常体的定位.

测试模型二为一个三维的模型，模型背景密度值为0，其中存在两个阶梯型的异常体，其密度为恒定的0.1 g/cm³，模型被剖分成21×21×21=9261个长方体，单个长方体大小为50 m×50 m×25 m.在地表采集重力数据，采集网格为21×21个，间距为50 m.模型和数据如图 4.

图 4

图 4 3D模型(a)及其正演数据(b).Fig. 4 3D model(a) and the forward gravity data.

应用本文方法对数据进行反演计算，图 5是概率成像结果，其中彩色部分为成像结果三维体切片，切片上方的黑线是模型的真正位置，从反演结果可以看出，基本能够确定有异常体的存在，但整体聚焦效果较差.图 6是本文的反演方法结果，对比图 5可以看出，聚焦效果有进一步的提高，并且获得的是实际密度物性的分布.

图 5

图 5 概率层析成像切片 (a)y=350 m;(b)y=550 m;(c)z=100 m;(d)=175 m.Fig. 5 Probability tomography imaging slices (a)y=350 m;(b)y=550 m;(c)z=100 m;(d)=175 m.

图 6

图 6 本文反演数据体切片 (a)y=350 m;(b)y=550 m;(c)z=100 m;(d)=175 m.Fig. 6 inversion result slices (a)y=350 m;(b)y=550 m;(c)z=100 m;(d)=175 m.

为了缩小反演解的范围，减少解的非唯一性，提高反演分辨率，重力反演计算经常加入已知的地质信息对反演进行约束.最简单的办法就是为每个反演的长方体设置一个解的区间[ρ_min，ρ_max]，在反演迭代计算中，设置如下的限制：

我们将这个方法应用于模型二的反演过程中，其中为每个长方体设置一个统一的区间为[0，0.15]，则图 7是加入此限制的反演结果.对比图 6，反演的聚焦性得到大幅度的提高，框定异常体的范围也进一步缩小，同时反演的重力值也更接近于模型的真实重力值.

图 7

图 7 加约束后本文反演切片 (a)y=350 m；(b)y=550 m；(c)z=100 m；(d)=175 m.Fig. 7 constrained inversion result (a)y=350 m；(b)y=550 m；(c)z=100 m；(d)=175 m.

4 结 论

本文提出了一种迭代算法来实现等效密度到真正密度的转化.其以采集数据概率成像结果乘以一个小的密度值为初始反演模型，并将正演和采集数据差的概率成像结果乘以一个密度值作为模型的更新值，并在反演过程中加入已知信息的约束，反演迭代直到正演结果和反演结果的相对误差满足某个设定的限制为止.通过两个模型的测试，本文方法实现了概率成像结果到密度转化的同时，也增加了异常体成像的聚焦效果，更好的框定了异常范围.该计算保留了概率层析计算稳定、简单的的特点，计算效率是概率层析的迭代次数倍，该方法可以成为通用方法的另一个选择，或可为通用方法提供一个初始的计算模型.

致 谢 感谢审稿专家和编辑部老师的指导和帮助.