首先确定观测数据与模型参数之间的函数关系，构建目标函数；为使改进后的算法有效应用于实测数据的反演，用L曲线的拐点值代替拉格朗日经验参数；为对模型空间结构加以约束，引入粗糙度矩阵，从而改善反演中存在的病态性问题；引入深度加权函数，使核函数可以反映出各个深度真实的权值；最终利用改进后的预条件共轭梯度法作为反演算法，高效地完成反演计算。

反演中，目标函数φ由数据拟合差函数φ_d和模型目标函数φ_m构成，解决反演问题即在数据拟合差函数满足条件的前提下，寻找模型向量m使求得的目标函数φ值最小^[1]：

式中，Δm=m-m₀代表模型参数向量m与初始模型m₀之间的改正量；G代表核函数，即表示从模型单元到观测点的线性投影算子；Δd代表相对应的观测值的改正量；W_d为一个对角矩阵; σ_i为第i个数据的标准偏差。

预条件共轭梯度反演方法中，使用拉格朗日乘子作为正则化参数。当求解大规模矩阵时，对拉格朗日乘子的求解需要大量运算，采用经验值作为正则化参数的相对最佳取值。由于正则化参数起到平衡数据拟合泛函与模型拟合泛函的作用，取值过大，会导致模型响应与观测值相差较大; 取值过小，会导致模型与实际情况不符, 使反演结果不可靠。

本文使用L曲线的曲率拐点值作为正则化参数。L曲线是基于实际数据拟合与模型目标函数之间进行权衡的一种准则，适合大规模问题的求解。当把模型目标函数与数据拟合函数画在对数-对数图上，所呈现的曲线有一个明显的拐点; 当正则化参数值大于拐点时, 模型目标函数变化不大, 数据拟合差函数变化很大；当正则化参数值小于拐点时，模型目标函数快速增长，而数据拟合差函数呈缓慢递减趋势^[8]。所以将拐点作为最佳的正则化参数，L曲线的曲率k定义为:

式中，$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \rho } = \ln {\varphi _d};\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \eta } = \ln {\varphi _m}; $符号“′”和“″”分别表示函数的一阶和二阶导数。因此反演中会根据实际数据所构建的模型寻求L曲线的最大曲率。

为了使反演图像能够在3个方向上连续变化，光滑过渡，将粗糙度矩阵引入到重力梯度张量的三维成像反演中^[9]。定义三维模型的粗糙度矩阵R为模型向量m在x、y、z 3个方向上一阶偏微分的平方和:

式中，v为模型向量所在的体积空间。

由于重力梯度值与场源到观测点距离的立方成反比，重力梯度数据没有固有的深度分辨率，所构建的核函数矩阵为线性的；随着深度的增加，核函数快速衰减，使反演结果集中于地表附近，难以反映出深部异常体的真实位置。通过引入深度加权因子，可以抵消核函数的快速衰减，采用文献[10]方法引入深度加权函数如下：

式中，Z₀取决于观测面的高度以及单元体的尺寸大小；Z为每层网格单元中心点的埋深；β为常数，对于重力梯度数据β=3，函数图像与重力梯度核函数图像衰减一致。

在网格单元划分时，浅层采用高分辨率格网，深层为较低分辨率格网，可以有效增加深层块体在地表处的重力梯度效应，减弱核函数的快速衰减，同时也可以减少反演模型的参数个数，改善反演中存在的病态性问题。

将粗糙度矩阵以及深度加权函数代入到模型目标函数中^[11]：

式中，$ {\alpha _s}、{\alpha _x}、{\alpha _y}、{\alpha _z}$为模型目标函数中各项的加权因子。

通过利用有限差分法代替微分，将式(7)改写成矩阵形式：

式中，α_i为目标函数中各项的权重系数；R_i为每一个分量的差分算子；D代表离散化的深度加权函数矩阵；s、x、y、z分别表示式(7)中加权因子参与计算的4个分量。将构建的模型函数代入到目标函数中：

改写成方程组形式：

用$\mathit{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{G}}\\ {\sqrt \mu {\mathit{\boldsymbol{W}}_i}} \end{array}} \right]$代表雅克比矩阵，$\mathit{\boldsymbol{b}} = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{d}}\\ {\bf{0}} \end{array}} \right]$，则式(11)写成：

位场反演中存在病态性与非唯一性，且经常需要对大规模问题求解，而共轭梯度法(conjugate gradient algorithm, CGA)只需利用函数的一阶导数信息，存储量小，算法简便，收敛速度比最速下降法快，适合求解大规模问题。CGA主要用于求解二次方程$\varphi = \frac{1}{2}{\mathit{\boldsymbol{m}}^{\rm{T}}}\mathit{\boldsymbol{Am}} - {\mathit{\boldsymbol{m}}^{\rm{T}}}\mathit{\boldsymbol{b}}$，解为$\mathit{\boldsymbol{m = }}{\mathit{\boldsymbol{A}}^{ - 1}}\mathit{\boldsymbol{b}}$，与方程(12)等价。该算法利用初始点的梯度作为最初的共轭方向，沿着这组方向而非负梯度方向去搜索目标函数的极小点，改善解的稳定性^[12]。反演过程如下：

1) 已知观测值d^obs和初始模型m₀^k, 其中k=1, 2 …N_max, N_max为外循环最大迭代次数。

2) 计算观测数据和由初始模型计算得到的理论数据之差Δd^k=d^obs-Gm₀^k。

3) 设定初始变量: $\mathit{\boldsymbol{A}} = {\left[ {\mathit{\boldsymbol{G}}\sqrt {{\mu ^{ - 1}}} {\mathit{\boldsymbol{W}}_i}} \right]^{\rm{T}}}, {\mathit{\boldsymbol{b}}^k} = \left[ {\Delta {\mathit{\boldsymbol{d}}^k}\;\;0} \right], \Delta \mathit{\boldsymbol{m}}_0^k = {\bf{0}}, {\mathit{\boldsymbol{r}}_0} = {\mathit{\boldsymbol{A}}^{\rm{T}}}\left( {\mathit{\boldsymbol{b}} - \mathit{\boldsymbol{A}}\Delta \mathit{\boldsymbol{m}}_0^k} \right), {\mathit{\boldsymbol{z}}_0} = \mathit{\boldsymbol{S}}{\mathit{\boldsymbol{r}}_0}, {\mathit{\boldsymbol{p}}_0} = {\mathit{\boldsymbol{z}}_0}, \mathit{\boldsymbol{S}}$代表预条件矩阵。

4) 设定i=1, 2 …N_CGmax, 其中N_CGmax为内循环最大迭代次数，进行预条件共轭梯度算法计算，其中q_i=Ap_i, α_i=(r_i-1^Tz_i-1)/q_i^Tq_i, Δm_i^k=Δm_i-1^k+α_ip_i, r_i=r_i-1-α_iA^Tq_i, z_i=Sr_i。z_k为改善的搜索方向; β_i=(r_i^Tz_i)/(r_i-1^Tz_i-1), p_i=z_i-1+β_ip_i-1。

5) 共轭梯度反演内循环计算结束, 设定m₁^k=m₀^k+Δm^k。

6) 设定m₀^k=m₁^k, 重建初始模型并返回步骤2)，继续迭代计算，直到满足收敛条件。

预条件矩阵S=(A^TA)^-1可以改善条件数，有效提高迭代速度，但求解大规模矩阵时对矩阵S的求解较为困难，在实际中用深度加权函数W(z)作为预条件因子，同样可以达到改善条件数的作用^[11]。

为验证改进后的预条件共轭梯度算法应用于密度结构反演的有效性，利用加入随机噪声的重力梯度响应作为观测值，并将地下模型改进为不等间隔模型，经过单分量反演后，将结果和真实模型进行比较分析。为充分利用重力梯度张量数据，改善反演中解的非唯一性，将5个独立分量进行联合反演，与单分量的反演结果对比，突出联合反演的优势。为表明将反演算法应用于航空重力梯度数据的可靠性，设计的真实模型中包含多个形状大小、密度值和质心埋深各异的异常体，对其进行联合反演并对结果分析。

为了测试模型目标函数中依次加入正则化参数、粗糙度矩阵以及深度加权函数对反演结果的影响，根据构建的模型进行反演试算。图 1(a)代表所构建的地下模型是1 000×1 000 m²区域范围内埋藏的单个异常体，记为异常体A。将地下模型划分为10层，使浅层模型单元的中心点与观测点所在的位置对应，深层模型单元的宽度逐渐变大；1~3层模型单元的行数及列数均为25，4~5层为20，6~8层为15，9~10层为10；每层模型单元大小从40 m递增到100 m，共由3 550个模型单元构成。模型尺寸统计信息见表 1。地表观测点的间隔为50 m，总共有400个观测点，观测面高度为80 m，初始模型矩阵m的值为0，将加入5%随机噪声的T_zz作为观测值，对反演结果进行比较分析。

图  1  反演算法测试

Figure 1.  Test of Inversion Algorithm

图 1(b)仅加入正则化参数，结果显示异常体集中于地表附近。图 1(c)再加入粗糙度矩阵，可以明显改善模型单元体的平滑度，避免相邻单元密度值的突跳现象。图 1(d)在此基础上加入深度加权函数，克服重力场固有的“趋肤”效应，使异常体归置于相应深度。将三者结合后，反演效果最佳，异常体集中，更接近于真实模型。

重力勘探中，重力梯度数据在浅层比重力数据具有更高的分辨能力。虽然从理论上同一方法不同分量无法改变由位场性质决定的多解性，但在同一高度不同分量的观测数据可以使数据更加丰富完整，改善反演中的多解性^[13]。模型统计信息见表 1红色异常体B和蓝色异常体C。初始模型矩阵m的值为0，在观测值中加入5%的随机噪声。图 2(a)代表所构建的地下模型是1 000×1 000 m²区域范围内埋藏的两个并排的正负异常体，图 2(b)~2(g)列出了重力梯度张量单分量的反演结果。

图  2  单分量反演与多分量联合反演结果及其分析

Figure 2.  Inversion of Independent Component and Joint-Inversion of Five Independent Components

由图 2(b)至图 2(g)可知，与真实模型相比，所得异常体的剩余密度值、形状、位置、边界均存在不同程度的偏差，T_xx、T_yy、T_xz、T_yz、T_zz为5个单独分量。利用T_xx和T_yy得到的反演结果在水平方向上互补，但反演得到的剩余密度值偏差较大，振幅可达±0.5 g/cm³；利用T_xz和T_yz得到的反演结果集中于异常周围, 可以反映地下密度异常体的边界信息，得到的剩余密度值偏低; T_zz反演结果与真实模型对应较好, 但异常体中心位置偏下，得到的剩余密度值偏低。为充分利用各单分量反演结果的优势并改善反演中解的不唯一性，利用重力梯度的5个分量联合反演的结果如图 2(h)所示，得到的反演结果更精确，异常体在深度方向上更集中于中心，剩余密度值与真实模型密度值更接近，边界明显更容易确定出异常体的形状。为进一步分析反演结果精度，图 2(i)代表包含噪声的观测值T_zz，等值线代表由联合反演得到的剩余密度值计算出的T_zz′，图形吻合较好，差值较小; 计算T_zz和T_zz′的残差值，图 2(j)代表两者之间的残差直方图，残差值的标准差为±0.977 7 E，表明反演结果可靠、有效。

为了增强航测重力梯度数据反演结果的可靠性，增加真实模型的复杂性，设计的模型由多个形状大小各异的异常体组成。图 3(a)代表所构建的地下模型是1 000×1 000 m²区域范围内埋藏的多个异常体。初始矩阵m的值为0，在观测值中加入5%随机噪声。利用重力梯度的5个分量进行联合反演，图 3(b)~3(c)表示真实模型在不同深度的横截面；图 3(d)~3(f)表示联合反演结果及其在不同深度的横截面；图 3(g)中的等值线为联合反演结果计算得到的T_zz′；图 3(h)为联合反演结果计算得到的T_zz′与观测值T_zz的残差直方图。

图  3  原始模型与反演结果对比及其分析

Figure 3.  Comparison Between Original Model and Inversion Result with Analysis of Inversion Results

由图 3(e)可以看出, 反演结果在浅部反映出异常体的形状和位置更接近真实模型。随着深度的增加，由于梯形异常体由几个异常块体叠加而成，形状较为复杂，由图 3(f)可以看出反演结果的形状、位置和真实模型稍有偏差，小的异常体反演结果的形状、位置和真实模型较为相近。整体反演结果可以准确反映出异常体的数量、形状、大小及其所处的中心位置。对反演结果精度的分析表明，图 3(g)中观测值T_zz和剩余密度值m计算出的T_zz′之间的差值较小，图 3(h)中残差值的标准差为±1.003 1 E，在有效范围内，表明算法可用于复杂模型的反演计算。

实测数据来源于澳大利亚西澳洲地区2009年建立的航空重力/重力梯度测量Kauring试验场，西澳洲矿产资源丰富，以富产铁矿、金矿、氧化铝等而全球闻名，且具有储藏量大、分布广、埋藏浅和易开发等特点。该试验场位于珀斯简达科特机场东北方向115 km处(图 4)，该地区地势平坦，有利于试验的开展。

图  4  Kauring重力梯度试验场位置及该地区地形图

Figure 4.  Location of Kauring Gravity Gradient Test Site and the Topographic Map of Area

Fugro Airborne Surveys于2011年7月至2012年2月间，利用FALCON航空重力梯度仪(airborne gravity gradient, AGG)在该地区进行了航空重力梯度测量，在Kauring试验场飞行的总测线长度为1 265.3 km，飞行垂直高度为80 m，使用的坐标系统为参心大地坐标系。实测重力梯度数据已经过地形改正，为突出浅部异常，经过参数试验并结合之前的研究成果^[14]，图 5(a)~5(f)使用相匹配的滤波法消除由深部异常产生的影响。北向范围为6 469~6 473 km，东向范围为504~508 km，所选测区范围大小为4 000 m×4 000 m。

图  5  巴特沃斯四阶低通滤波后的重力梯度全张量

Figure 5.  Gravity Gradient Tensor by Using Butterworth Four Order Low Pass Filter

实测航空重力梯度数据的分辨率为5.1 m，精度为5 E；在采集数据中共选取80×80个点作为观测数据，观测点间隔为50 m，利用5个独立分量进行联合反演，故由观测值组成的矩阵大小为32 000×1；将模型划分为10层，每层厚度为80 m。参照数值模拟中模型设计的方法，地下模型单元大小从80 m递增到150 m，地下模型共由11 546个模型单元组成，核函数G为32 000×11 546的矩阵，初始模型矩阵m的值设为0。图 6(a)给出反演结果的三维模型图，图 6(b)至图 6(d)给出异常体中心位置的横纵截面。为突出剩余密度异常值较大的异常体，在三维模型图上仅给出密度值大于1.0×10³ kg/m³的异常块体。相比之前的结果，在反演得到的三维模型图 6(a)中发现中心异常块体周围还存在着多个大小不等的密度异常体单元，与之前已发现的中心块体沿线呈一定角度分布，这些异常体的密度值在1.25×10³ kg/m³左右。在y=1 700 m处给出纵截面图，图 6(b)表明反演结果的中心异常块体位于地下深度240~400 m范围，反演出异常体的中心位置与刘金钊等^[15]之前对该地区密度异常深度探测的结果相吻合。反演中心块体位置位于坐标系正东方向505.11~505.45 km范围内，正北方向6 470.45~6 470.89 km范围内，反演出中心块体的密度值在2.3×10³ kg/m³左右，与之前用航空重力反演该地区异常体的密度值相近，验证了反演方法应用于实测数据的可靠性。

图  6  反演结果

Figure 6.  Inversion Results

由于从反演出的三维模型图上看出异常块体主要位于第五层320~400 m范围，所以给出360 m的横截面切片图。从图 6(c)发现，几块明显的异常块体沿坐标点(2 000，0)呈arctan0.5=26.5°分布，异常块体分布的沿线距离长度在2.3 km左右。如图 6(d)所示，可将横截面显示出的异常块体划分为4块，反演出异常块体的中心分别位于(505 280, 6 470 670)、(505 610, 6 470 170)、(505 950, 6 469 500)和(506 170, 6 469 110)点处，异常体大小从300 m×300 m至300 m×500 m不等。由于异常块体中心深度位于360 m左右，距离地表较近，且该地区的地势平坦，在之前的地质研究中并未发现断层的存在，因反演的初值m均为0，并未加任何约束，反演出异常体的密度值最大可达到2.41×10³ kg/m³。而该地区矿藏具有埋藏浅的特点，故推断是由于矿体的存在导致该地区存在明显的异常体。对比之前重力反演的结果^[14]，发现中心块体周围还存在着明显的异常块体。

利用改进算法，可反演出异常体的位置、大小、密度异常值，从而判定密度异常体的分布特征。

本文利用改进后的预条件共轭梯度算法对地下密度异常体进行三维密度结构反演，通过构建模型重力梯度场验证了方法的有效性和可靠性。并利用该算法对澳大利亚西澳洲Kauring地区实测重力梯度场进行了反演，给出了测区内密度异常体的分布信息。国内航空重力梯度测量还未大范围开展, 希望通过本文对国外实测航空重力梯度数据的处理与解算, 为今后国内重力梯度试验场数据的应用提供经验和基础；该方法在重力梯度矿产资源调查中具有广泛的应用前景，可以得到更加完整的密度异常体分布信息，从而进行更加深入的地质结构分析和解释工作。