0 引言

在地球物理勘探中，模型体的正演在位场异常的解释中有着重要的意义，对其产生异常特征的认识是掌握场与场源对应关系的切入点，是进行位场反演、地质解释的基础，历来深受地球物理学家的关注。重力矢量异常不仅反映重力在垂直方向上的变化，而且提供了更加丰富更加细致乃至于能反映地球时变的中高频重力场信息^{[1, 2]}。

针对三度体重力异常的正演计算，目前主要是通过已有的解析公式直接计算、逼近目标体的方法来实现。这类方法以牛顿万有引力定律为出发点，对具有简单形态和均匀密度分布的地质体直接利用解析公式计算获得；对具有复杂形态和非均匀密度分布的地质体采用不同的方式对复杂形体进行分割，使之转化为一系列简单形体的组合，计算简单形体的重力异常，通过求和得到整个复杂形体的异常，如基于直立矩形柱体、直立多边形棱柱体、多面体模型的三维正演技术^{[3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]}。一般而言，利用直立矩形柱体模拟复杂的三维非均匀密度模型是最简单的，如果矩形体足够小，则每个矩形体均可看成具有单一的密度分布。该方法虽然理论上简单，但是在实际应用中却稍显笨重。一方面是因为实际的地质体通常难以用矩形柱体进行建模；另一方面，如果与某个矩形柱体相邻的其他柱体与该矩形柱体具有相同的密度分布，则完全没必要将其他柱体切割成矩形柱体进行重复计算。为改进直立矩形柱体的上述缺点，Talwani^[4]提出了直立多边形棱柱体模型，通过对一系列无限薄的水平薄板的叠加模拟实际地质体，每个薄板的形状用多边形近似。但该方法对具有任意形状的地质体模拟却稍显不足。为此，Okabe^[7]提出了更加实用的针对具有任意形状地质体的多面体法，用一系列不同的多边形围成的多面体来逼近任意形状的地质体，其逼近的程度取决于围成多面体多边形的数量与顶点的选取。由此可见，基于均匀直立矩形柱体、直立多边形棱柱体、多面体模型模拟复杂地质形体的三维正演技术包括了两个方面的近似：一是分割方式与实际形体的拟合程度；二是数值积分代替解析积分的近似程度。这种基于简单规则几何体剖分的重力异常计算方法，无论是计算精度还是计算效率均不高，易导致对后续的反演处理与解释产生严重影响。

有限元法作为地球物理数据正演的方法之一，最大优势在于可计算任意形状、任意密度的复杂地质体。随着计算机技术的发展，CAD(computer aided design)/CAE(computer aided engineering)的集成及ANSYS等大型有限元软件提供了强大的有限元建模及网格剖分功能。对于三维变密度的复杂地质体模型，只要将连续的求解区域离散为一组有限个、且按一定方式互相联结在一起的单元组合体，再将不同单元赋以不同密度值，就可以模型化形状复杂的求解域，通过叠加原理，很快就能计算出不同模型在整个空间的重力矢量。目前，基于有限单元法对三维地质体重力矢量的正演研究还未见相关文献报道。但是仅仅针对重力异常及二维常密度地质体梯度张量的正演计算却取得了可喜的进展^{[14, 15, 16, 17, 18]}，如Cai和Wang^[15]以引力位在场源外满足拉普拉斯方程、在场源内满足泊松方程以及相关边界条件为出发点，利用有限元法能真实地再现复杂情况下地球物理场的分布，有效地处理不均匀介质和求解边界形状复杂问题的独特性能，实现了利用有限元法计算非均匀复杂形体重力异常的快速算法。其模型试验表明，有限元法受密度分布的非均匀性影响远远小于传统的数值积分法，而且在网格单元相同时，极大地提高了重力异常正演计算的精度；更为可贵的是相对传统的数值积分法，计算效率得到了极大的提高。

本文依据重力位的泊松方程，得到与计算三度体重力矢量边值问题对应的变分问题，进而利用有限元法实现了对三度体重力矢量的正演计算。

1 三度体重力矢量的变分问题

1.1 边值问题

在位场勘探中，引力源一般位于地下，在地表上部的无源空间内满足拉普拉斯方程，在场源内部满足泊松方程，即：

式中：U为重力位；G为引力常数；ρ是场源体与围岩的密度差。对式(2)分别在x，y和z方向求导，则有 式中：式(3)即为重力矢量满足的偏微分方程。

如果给定区域Ω上的ρ分布以及边界S上重力矢量的某类边界条件，例如第三类边界条件，则求解区域内的U_i与以下的边值问题相对应：

式中： n为边界S的外法向方向；α、β 是常量参数。

1.2 变分问题

利用有限单元法求解边值问题式(3)中的U_i，与之相适应的变分问题为对泛函F(U_i)求极值(推导过程见附录A)，即

需要指出的是,在对式(4)求解的过程中，参数α,β是未知的。在本文的研究中做以下近似处理：将边界区域S取得足够远，这时场源在边界产生的重力位可看作是质量集中于场源体质心的质点在边界上的重力位，可表示为 式中：c为常数；r为质心到边界单元面中心的距离。则对式(6)在X,Y,Z方向上求偏导数可得 进一步地，U_i对边界面法线方向n的偏导数为 由此可知

在实际计算中，式(9)中α第二项的分母i代表x，y，z，可表示为x=rcosθcosD,y=rcosθsinD,z=rsinθ，如图 1所示。

将式(9)代入式(5)中即可得到与边值问题式(4)相对应的变分问题：

进而可利用有限单元法来求解式(10)的变分问题。

2 数值实现

2.1 单元插值函数及积分实现

取如图 2a所示的六面体母单元，取8个角点为节点，其编号及坐标如图中所示。则可构造如下形函数：

其中：(ξ_j,η_j,ζ_j)是j点(j=1,2,…,8)的坐标；ξ,η,ζ∈^{[-1, 1]}。单元中U_i的插值函数是 式中，U_i,j是六面体单元顶点j的待定U_i值。

为避免在等参单元积分中高斯数值积分点与计算节点不在同一点的情况，在本文的研究中，根据六面体单元的特点，令，则

式中：a，b，c为六面体子单元的棱边长；x₀，y₀，z₀是子单元的中点(图 2b)。据此，对式(10)第一个积分中的第一项计算可表示为 其中：为六面体子单元中8个节点处的待定U_i 值；；j,k=1,2,…,8。将式(11)代入 K_{1 e}中进行积分，即可得到K_j,k。且 K _1,e是对称矩阵，K_j,k=K_k,j。

式(10)中第一个积分中的第二项可表示为

式中

对式(10)中的第二个积分进行计算，则可得到以下结果：

若单元e 的底面位于边界面上，则 K _2e 具有如下的形式：

由此，将各个单元的节点组成的矩阵相加，得

式中：K和P是由所有节点构成的扩展矩阵。故而，由式(17)对式(10)取极值，且考虑总体合成矩阵 K 的对称性，有

式中，δ表示变分。由于δU_i 的任意性，得线性方程组： 解线性方程组(18)可得各个节点的重力矢量 U _i 。

2.2 数据存储方式

对任何地球物理问题的有限元分析最终都归结为求解大型线性方程组问题，即对具有式(18)形式的方程组的求解。对于三维问题，一般情况下系数矩阵 K 的维数非常大。例如，将求解区域分为60×60×60的网格单元，则有226 981个网格节点，形成的矩阵 K 大小为226981×226981。如此巨大的矩阵，常规存储方式很难实现其数值计算。因而系数矩阵的存储结构及方式成为有限元求解效率的关键之一。考虑到在地球物理问题中系数矩阵的稀疏性及对称性，很多学者提出了有效的节约存储空间的存储方法，如半带宽、变带宽、coordinate storage、CRS(compressed row storage)等^{[19, 20]}。为进一步节约存储空间，笔者根据求解重力矢量三维有限元的系数矩阵特点，提出一种适合迭代算法的改进一维变带宽压缩存储方法。因系数矩阵的对称性，在存储过程中，仅存储矩阵的下三角部分。假设矩阵 A 有如下的形式：

则笔者提出的改进存储方法可按下列步骤实现：

1)用一维矩阵 B 按行依次存储各非零元素，非零元素的总个数为N_d。在 B 的N_d+1 位置，存放矩阵 A 的行数。

2)生成另一矩阵 I_A 用来存放两个相邻的非零元素间的间隔数，I_A 的长度为N_d 。

由此，A 的存储可表示为表 1。

需要指出的是，笔者提出的存储方法一方面避免了CRS等存储方法中需要单独存入对角线元素位置及列号的缺点，特别当 A 数组维数较大时，存储的列号数据要远远大于本文中 I_A 的数据。此外，就存储的数据量而言，CRS为2N_d+N，本文为2N_d+1，较之CSR存储方法，笔者提出的存储方法 I_A 的数据小，进一步地节约了存储空间。另一方面，在计算 A 的转置矩阵时是非常方便的，只要将 B 中的数据按列，以 I_A 所示的间隔数依次存储即可得到。由此可见，利用笔者提出的存储方法更易于得到上三角阵，这为利用预条件处理技术求解线性方程组提供了方便。

2.3 预条件处理技术

对式(18)所示的大型稀疏线性方程组，利用不完全Cholesky共轭梯度(ICCG)迭代方法能够取得较好的效果，将其应用于地球物理正演具有较大潜力。对于ICCG法，文献^{[21, 22]}中均有详细的描述，文中不再赘言，仅给出其迭代过程。

其中：j=0,1,2,…为迭代次数，例如 U _ij即表示第j次迭代的 U _i值； C 为下三角矩阵，是不完全Cholesky分解因子，即 K ≈ CC ^T。需要指出的是，在ICCG方法中，C矩阵与K 矩阵具有完全一样的稀疏性，仅对角线元素值不同，从而能够简单快速地计算出 C 矩阵；此外，在存储过程中，不需要另外的存储空间来存储 C 矩阵。

3 模型试验

前文从理论上说明了利用有限元法实现对重力矢量计算的可行性。下面以顶面埋深为20 m、边长为40 m的理想立方体模型为例，研究利用有限元法对重力矢量进行正演计算的可靠性，探讨网格剖分间距大小及边界距离立方体中心远近对计算结果的影响。计算过程中立方体的剩余密度取为1.0 g/cm³。网格的节点编号如图 3所示。

图 4给出了利用有限单元法计算的重力矢量。由图 4可以看出异常具有明显的对称性，这与立方体模型的重力矢量异常分布相符，证明了计算结果的正确性。图 4中：U_x异常等值线呈椭圆形分布，南高北低，有两个符号相反的极值中心，极值大小分别为0.986 4和-0.986 4 g.u.；而U_y异常则沿东西向呈对称状分布，其余特征与U_x类似，不再赘述；U_z异常等值线呈近圆形状分布，有一个极大值，位置在立方体中心的地表投影处，大小为2.465 9 g.u.，向四周等值线逾稀疏。

误差限为0.001时，不同立方体单元边长及单元总数、总体系数矩阵非零元素数、占用空间、计算用时等参数见表 2。表 2中的计算用时均为程序在主频为3.10 GHz、内存为3.25 GB的微机上运行时间，边界长度指的是2倍边界单元的外表面到场源体中心的垂直距离。从表 2中可以看出：随着单元数及节点数的增加，总体系数矩阵中非零元素数倍增，相应的存储空间亦呈逐渐增大的趋势；因采用前文中所述及的存储方法，整体系数矩阵及其联系数组所占用的空间并不大，当节点总数为226 981时，所占用的存储空间仅为23.49 Mb，一般的微机均能满足存储要求。此外，从ICCG迭代次数看，计算U_x的迭代次数要大于U_z，但二者差别不大，在3次左右；从计算效率看，随计算节点总数的增加，计算用时愈来愈多，且计算U_x的用时要略大于U_z。

较之传统算法中的棱柱体法而言^[15]，若将本文模型中的立方体按照文中提出的有限元单元法计算过程中所划分的单元数量，计算同等数量的节点重力矢量，本文的计算用时要远小于棱柱体法，详见表 2。究其原因是因为在本文的方法中，假如求解区域为M个单元，场源体为N个单元，二者的比值M/N=Q为10~100，则可以同时计算出所有节点处的重力矢量，所需的积分次数则为M。而棱柱体算法中，计算任意一点的重力矢量均需要进行N次积分运算，若求解区域的节点数为P，则所要进行的积分总数为N·P。这就表明比值N·P/M=P/Q可以用来作为衡量棱柱体算法与本文提出的方法之间计算效率的标准。由此可见，文中提出的方法就计算效率而言，要远远高于传统的棱柱体算法。

4 计算精度分析

考虑到前文中立方体模型U_x与U_y的呈转置关系的特点，文中仅分析U_x和U_y的计算精度。为方便精度分析，选取过U_x和U_y极值中心的主剖面AB(图 4)作为分析对象。图 5给出了边界长度为200 m，立方体单元边长分别为20、10、5和4 m的有限元法计算结果及其与理论值的相对误差和相对百分比误差，其中模型的重力矢量理论值按文献^[23]推导的理论计算公式计算。由图 5a可以看出：有限元法的计算的U_z值小于理论值，而U_x总体上呈正值端小于理论值、负值端略大于理论值、0值处与理论值相等的特征；此外，随着单元立方体边长的减小，其与理论值的拟合程度逾来逾高。由图 5b可以看出：U_z的相对误差绝对值在边界处较小，而在极值处则与理论值有较大偏差，总体呈两端平直，中部略微下凹的特征；而U_x则在总体上呈线性特征。从百分比误差看，U_x和U_z的百分比误差绝对值均在端部较大，而在观测平面的0值附近，二者的百分比绝对值误差均小于5%(图 5c)。

图 6给出了单元立方体边长为5 m，边界长度分别为100、150、200、250和300 m时的有限元计算结果及其与理论值的相对误差和相对百分比误差。由图 6a可以看出，随着边界长度的增加，U_x和U_z与理论值的拟合程度逐渐增高。这在相对误差曲线(图 6b)上表现得尤为明显： U_z误差曲线在边界长度较小时呈“V”型特征，随着边界长度的增加，这种特征逐渐消失，曲线呈逐渐被拉直的特征；U_x误差曲线则呈过0点的斜率大于0的近直线特征，边界长度越大，直线的斜率越小。从百分比误差曲线(图 6c)看：U_x和U_z均呈现出倒“U”型的特征，其中U_x的“U”型曲线陡度大于U_z，即其在靠近端部的百分比误差绝对值要大于U_z；随着边界长度的增加，二者的倒“U”型曲线底部越来越平直，百分比误差绝对值越来越小。需要特别指出的是，当边界长度增加到200 m时，再增加边界长度时，U_x和U_z端部的百分比误差减小得并不多，表明当边界长度达到5倍场源体边长时，边界长度对U_x和U_z端部值的影响越来越小。为进一步说明边界长度对有限元法计算结果的影响，文中统计了在0~50 m范围内的百分比误差数值，如表 3所示。

由以上可以看出，利用有限单元法计算重力矢量的计算精度主要受剖分网格单元的边长大小及边界长度的影响。一般情况下，边界处的百分比误差较大。究其原因，其一是因为在计算边界单元的贡献时，假定了场源在边界产生的重力位是质量集中于场源体质心的质点在边界上的重力位；其二是因为边界取的仍不够远；其三则因为在本文的研究中单元立方体节点间插值函数取的是线性插值，相对精度不是太高。若要获得理想的计算结果，一方面要考虑网格单元的大小、网格单元总数，另一方面需要权衡计算时间，这在实际应用中尤为重要。结合本文的计算结果，笔者建议对场源体的剖分尽量细化，单元的边长至少应小于场源体边长的1/10，而边界长度则至少应为场源体长度的7.5倍。

5 结论

1)文中提出的总体系数矩阵存储方式相对于半带宽、变带宽、coordinate storage、CRS等传统存储方式能够进一步地节约存储空间；且在利用ICCG法解方程组时，更易于得到系数矩阵的转置矩阵，这为更快速地得到线性方程组的解提供了保障。

2)重力矢量正演在采用同一插值函数的前提下，其精度主要取决于单元总数和单元网格的尺寸大小。其计算效率则与所要计算的节点总数息息相关。此外，对大型稀疏线性方程组的求解算法的优劣成为提高计算效率与否的关键。

3)利用有限单元法进行重力矢量的正演计算，当单元的边长小于场源体边长的1/10、边界长度大于场源体长度的7.5倍时，能够获得较为理想的结果，除边界附近百分比误差较大外，其余位置的百分误差一般小于1%。就计算效率而言，有限单元法要远远高于传统的棱柱体算法。

附录 A

针对式(4)的边值问题，构造泛函

则其变分为

根据格林公式

令U=δU_i，V=U_i，则(A-2)式的第一项可写为

又(A-2)式的第二项可写为

将(A-3)和(A-4)式代入(A-2)式中，可得

考虑到(A-5)式中的第一项因重力矢量满足的泊松方程，其值为0，第三项中由于在边界处无密度体分布，故而其值亦为0，因此，有

δI(U_i)-_SδU_iU_indS=0 。(A-6) 将(A-1)式中的I(U_i)以及第三类边界条件代入(A-5)式中，即可得到与(4)式表述的边值问题相对应的变分问题，即