三度体重力矢量的有限单元法正演计算

蒋甫玉，谢磊磊，常文凯，黄 岩，张作宏

1.河海大学地球科学与工程学院，南京 210098 2.江苏省地质勘查技术院，南京 210048



三度体重力矢量的有限单元法正演计算

蒋甫玉1，谢磊磊1，常文凯1，黄 岩2，张作宏2

1.河海大学地球科学与工程学院，南京 210098 2.江苏省地质勘查技术院，南京 210048

以重力位在场源内部满足泊松方程为依据，以重力矢量满足第三类边界条件为切入点，推导了与三度体重力矢量满足的边值问题相对应的变分问题，进而利用有限单元法实现了对变分问题的求解。立方体模型试验结果表明：文中提出的新的系数矩阵存储方式较之传统方式能够更有效地节约存储空间，且为利用预条件共轭梯度技术更加快速地求解线性方程组提供了保障；重力矢量的计算精度与边界长度及单元网格的边长息息相关，其计算效率则主要取决于所要计算的节点总数和大型稀疏线性方程组求解算法的优劣；一般情况下，当单元的边长小于场源体边长的1/10、边界长度大于场源体长度的7.5倍时，能够获得理想的结果。

变分问题；重力矢量；有限单元法；数据存储；计算精度；三度体

0 引言

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

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

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

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

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

1.1 边值问题

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

(1)

(2)

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

(3)

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

(4)

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

1.2 变分问题

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

(5)

需要指出的是,在对式(4)求解的过程中，参数α,β是未知的。在本文的研究中做以下近似处理：将边界区域S取得足够远，这时场源在边界产生的重力位可看作是质量集中于场源体质心的质点在边界上的重力位，可表示为

(6)

式中：c为常数；r为质心到边界单元面中心的距离。则对式(6)在x,y,z方向上求偏导数可得

(7)

进一步地，Ui对边界面法线方向n的偏导数为

(8)

由此可知

(9)

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

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

(10)

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

图1 计算重力矢量的区域Ω和边界SFig. 1 Region Ω and boundary S used in computing gravity vector

2 数值实现

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

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

(11)

其中：(ξj,ηj,ζj)是j点(j=1, 2, …, 8)的坐标；ξ,η,ζ∈[-1,1]。单元中Ui的插值函数是

(12)

式中，Ui,j是六面体单元顶点j的待定Ui值。

a. 母单元及其坐标；b. 子单元。图2 六面体母单元及子单元Fig.2 Mother units and subunits of hexahedron

(13)

式中：a，b，c为六面体子单元的棱边长；x0，y0，z0是子单元的中点(图2b)。据此，对式(10)第一个积分中的第一项计算可表示为

(14)

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

(15)

式中：

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

(16)

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

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

三是社交性。 不再像音乐电台那样，用户是被动且封闭的信息受传者，在音乐的网络传播过程中创作者和受众、创作者之间、受众之间可以及时进行交流。 Beats1作为苹果音乐下的分支，间接地享受到网络传播模式带来的良好社交性。

(17)

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

式中，δ表示变分。由于δUi的任意性，得线性方程组：

(18)

解线性方程组(18)可得各个节点的重力矢量Ui。

2.2 数据存储方式

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

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

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

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

表1 系数矩阵非零元素存储

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

2.3 预条件处理技术

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

令r0=-P-KUi0,q0=(CCT)-1r0, 则

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

3 模型试验

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

Nz、Ny分别为z方向和y方向一条直线内的节点总数。图3 单元立方体编号示意图Fig.3 Schematic of cell cube node number

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

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

AB为极值中心主剖面。图4 立方体模型重力矢量(单位：g.u.)Fig.4 Gravity vector of cube model (unit: g.u.)

单元总数节点总数总体系数矩阵非零元素数占用空间/Mb迭代次数UxUz计算用时/sUx本文算法棱柱体算法Uz本文算法棱柱体算法边界长度201031331155610.121191.22.21.12.1200m,不1020392611181210.90181536.072.832.672.2同立方体5403689219202417.0228252781.44528.52761.64518.6边长/m4503132651178780113.64322810942.225432.810821.625056.3立方体边10020392611181210.90191648.598.447.197.5长5m,不150303297913916812.992320441.81025.3423.21011.8同边界200403689219202417.0228252781.49875.62761.69859.3长度/m250503132651178780113.64322810847.630895.810799.230786.4300603226981307836123.49363242290.0109865.042100.2109400.0

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

4 计算精度分析

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

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

表3 0～50 m范围内的相对百分比误差统计

a.与理论值的对比；b.相对误差；c.相对百分比误差。图5 不同单元立方体边长的Uz和UxFig.5 Calculated Uz and Ux using different length of cell cube

a. 与理论值的对比；b. 相对误差；c. 相对百分比误差。图6 不同边界长度的Uz和UxFig.6 Calculated Uz and Ux using different boundary length

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

5 结论

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

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

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

附录 A

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

(A-1)

则其变分为

(A-2)

根据格林公式

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

(A-3)

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

(A-4)

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

(A-5)

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

(A-6)

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

(A-7)

[1] 曾华霖. 重力场与重力勘探[M]. 北京: 地质出版社, 2005. Zeng Hualin. Gravity Field and Gravity Exploration[M]. Beijing: Geological Publishing House, 2005.

[2] 李姗姗, 吴晓平, 吴星. 重力矢量及张量信息在地球重力场中的应用[J]. 测绘学院学报, 2005, 22(2): 94-96. Li Shanshan, Wu Xiaoping, Wu Xing. Applications of Gravity Vector and Tensor in the Gravity Field[J]. Journal of Institute of Surveying and Mapping, 2005, 22(2): 94-96.

[3] Nagy D. The Gravitational Attraction of a Right Rectangular Prism[J]. Geophysics, 1966, 31(2): 362-371.

[4] Talwani M. Computer Usage in the Computation of Gravity Anomalies[J]. Methods in Computational Physics, 1973, 13(1): 343-389.

[5] Kwok Y K. Singularities in Gravity Computation for Vertical Cylinders and Prisms[J]. Geophysical Journal International, 1991, 104(1): 1-10.

[6] Kwok Y K. Gravity Gradient Tensors Due to a Polyhedron with Polygonal Facets[J]. Geophysical Prospecting, 1991, 39(3): 435-443.

[7] Okabe M. Analytical Expressions for Gravity Anomalies Due to Homogeneous Polyhedral Bodies and Translations into Magnetic Anomalies[J].Geophysics, 1979, 44(4): 730-741.

[8] Holstein H, Ketteridge B. Gravimetric Analysis of Uniform Polyhedral[J]. Geophysics, 1996, 61(2): 357-364.

[9] Commer M. Three-Dimensional Gravity Modelling and Focusing Inversion Using Rectangular Meshes[J]. Geophysical Prospecting, 2011, 59(5): 966-979.

[10] Li X, Chouteau M. Three-Dimensional Gravity Mo-deling in All Space[J]. Surveys in Geophysics, 1998, 19(4): 339-368.

[11] 蒋甫玉, 高丽坤, 黄麟云. 油气模型的重力梯度张量研究[J]. 吉林大学学报: 地球科学版, 2011, 41(2): 545-551. Jiang Fuyu, Gao Likun, Huang Linyun. Study on Gravity Gradient Tensor of the Oil-Gas Model[J]. Journal of Jilin University: Earth Science Edition, 2011, 41(2): 545-551.

[12] 骆遥, 姚长利. 复杂形体重力场、梯度及磁场计算方法[J]. 地球科学: 中国地质大学学报, 2007, 32(4): 517-522. Luo Yao, Yao Changli. Forward Method for Gravity, Gravity Gradient and Magnetic Anomalies of Complex Body[J]. Earth Science: Journal of China University of Geosciences, 2007, 32(4): 517-522.

[13] 陈召曦, 孟小红, 郭良辉, 等. 基于GPU并行的重力、重力梯度三维正演快速计算及反演策略[J]. 地球物理学报, 2012, 55(12): 4069-4077. Chen Zhaoxi, Meng Xiaohong, Guo Lianghui, et al. Three Dimensional Fast Forward Modeling and the Inversion Strategy for Large Scale Gravity and Gravimetry Data Based on GPU[J]. Chinese Journal of Geophysics, 2012, 55(12): 4069-4077.

[14] Dave A M, Matthew G K. Optimal, Scalable Forward Models for Computing Gravity Anomalies[J]. Geophysical Journal International, 2011, 187(1): 161-177.

[15] Cai Yongen, Wang Chiyun. Fast Finite-Element Cal-culation of Gravity Anomaly in Complex Geolo-gical Regions[J]. Geophysical Journal International, 2005, 162(3): 696-708.

[16] Cruz F A, Knepley M G, Barba L A. PetFMM: A Dynamically Load-Balancing Parallel Fast Multipole Library[J]. International Journal for Numerical Methods in Engineering, 2010, 85(4): 403-428.

[17] Farquharson C G,Mosher C R W.Three-Dimensional Modelling of Gravity Data Using Finite Differences[J]. Journal of Applied Geophysics, 2009, 68(3): 417-422.

[18] 朱自强, 曾思红, 鲁光银, 等. 二度体的重力张量有限元正演模拟[J]. 物探与化探, 2010, 34(5): 668-672. Zhu Ziqiang, Zeng Sihong, Lu Guangyin, et al. Finite Element Forward and Simulation of the Two Dimensional Gravity Gradient Tensor[J]. Geophysical and Geochemical Exploration, 2010, 34(5): 668-672.

[19] 徐世浙. 地球物理中的有限单元法[M]. 北京: 科学出版社, 1994. Xu Shizhe. Finite Element Method in Geophysics[M]. Beijing: Science Press, 1994.

[20] 张力. 坑道直流聚焦超前探测有限元数值模拟研究[D]. 长沙：中南大学, 2011. Zhang Li. Research on Numerical Simulation for DC Focusing Advanced Detection in Tunnel with Finite Element Method[D]. Changsha: Central South University, 2011.

[21] 吴小平, 徐果明. 不完全Cholesky共轭梯度法及其在地电场计算中的应用[J]. 石油地球物理勘探, 1998, 33(1): 89-94. Wu Xiaoping, Xu Guoming. Incomplete Cholesky Conjugate Gradient Method and Its Application in Geoelectric Field Computation[J]. OGP, 1998, 33(1): 89-94.

[22] 吴小平, 徐果明, 李时灿. 利用不完全Cholesky共轭梯度法求解点源三维地电场[J]. 地球物理学报, 1998, 41(6): 848-855. Wu Xiaoping, Xu Guoming, Li Shican. The Calculation of Three-Dimensional Geoelectric Field of Point Source by Incomplete Cholesky Conjugate Gradient Method[J]. Chinese Joural of Geophysics, 1998, 41(6): 848-855.

[23] Nagy D, Papp G, Beneded J. The Gravitational Potential and Its Derivatives for the Prism[J]. Journal of Geodesy, 2000, 74: 552-560.

Forward Calculation of Three Dimensional Gravity Vector Using Finite Element Method

Jiang Fuyu1, Xie Leilei1, Chang Wenkai1, Huang Yan2, Zhang Zuohong2

1.SchoolofEarthSciencesandEngineering,HohaiUniversity,Nanjing210098,China2.GeologyExplorationTechnologyInstituteofJiangsuProvince,Nanjing210048,China

Variational problem of three dimensional gravity vector was deduced to meet the boundary value based on Poisson equation and the third boundary condition, and the solution of variational problem is further implemented by using the finite element method. The results of the cubic model test show that the proposed new coefficient matrix storage strategy is more effective to save storage space than a traditional approach; this, in turn, makes it possible to quickly solve liner equations by using the preconditioned conjugate gradient technology. The calculation precision of the gravity vector is closely related to the boundary length and unit grid; while the computational efficiency mainly depends on the total number of nodes and the algorithm used in solving a large sparse system of linear equation. In general, when the length of unit grid is less than 1/10 of the body length, and the boundary length is greater than 7.5 times of the length of the source, a desired result can be achieved.

variational problem; gravity vector; finite element method; data storage; calculation precision; three dimensional body

10.13278/j.cnki.jjuese.201504301.

2014-10-30

江苏省自然科学基金项目 (BK20140844)；江苏省地质矿产勘查局科研技改项目(2014-KY-15)

蒋甫玉(1981--)，男，讲师，博士，主要从事固体地球物理学研究，E-mail:jiangfy@hhu.edu.cn。

10.13278/j.cnki.jjuese.201504301

P631.1

A

蒋甫玉，谢磊磊，常文凯，等. 三度体重力矢量的有限单元法正演计算.吉林大学学报:地球科学版,2015,45(4):1217-1226.

Jiang Fuyu, Xie Leilei, Chang Wenkai, et al. Forward Calculation of Three Dimensional Gravity Vector Using Finite Element Method.Journal of Jilin University:Earth Science Edition,2015,45(4):1217-1226.doi:10.13278/j.cnki.jjuese.201504301.