0 引言

任意复杂条件下的重力异常的高效、高精度数值模拟，在重力异常反演成像及人机交互解释、建模中至关重要。目前重力异常的正演方法按照求解域的不同一般分为空间域和波数域两类。

在空间域正演中，Hubbert^[1]给出了二度体重力异常线积分公式；Talwani等^[2]推导了截面为多边形的常密度二度体重力异常解析式，并指出对于任意二度体都可以用该模型逼近；Rao^[3]推导了截面为矩形、密度随深度呈二次变化的重力异常解析表达式。后来关于截面为多边形、物性连续变化的二度体重力异常解析式的计算方法取得了较大进展^[4-8]。李明等^[9]采用级数展开计算了重力垂直一次导数；徐世浙^[10]在二维重力场垂直分量及梯度张量的正演计算中引入了有限单元法；姜效典等^[11]利用样条插值函数计算了变密度地质体重力异常；金刚燮等^[12]利用等参数有限元的高斯数值积分实现了复杂形体重力异常正演；Reeder等^[13]提出一种采用卷积算子提高二度体重力异常计算效率的有限元方法；Ren等^[14]提出了一种基于非结构化网格的重磁正演快速多极子算法，很好地解决了复杂边界和几何形状情况下重磁异常的正演问题；肖宝泽等^[15]采用均匀多边形截面二度体重力异常公式实现了复杂模型的重力异常数值模拟。

在波数域正演中，Sharma等^[16]给出了任意断层面倾角的断层重力场的波数域表达式；Rao等^[17]推导一个截面为等腰三角形的二度体模型的波数域表达式；Bhattacharyya等^[18]推导了任意二度体的重力异常场波数域表达式；在此基础上，Pedersen^[19]首次推导了以三角形组合作为其截面形状的二度体重力异常波数域表达式，并分析了在波数域如何以简单方式建模，从而解决波数域反演中的一些问题；吴宣志^[20]将均质物性的任意二度体模型推广至密度随深度呈线性或指数变化的模型；柴玉璞^[21]运用偏移抽样方法提高了反傅里叶变换数值精度；Tontini等^[22]实现了重力异常波数域正演，并研究了快速傅里叶变换(FFT)扩边法与误差的关系；Wu等^[23]利用Gauss-FFT提高了反傅里叶变换的数值精度，降低了FFT引起的强制周期化、边界震荡效应等影响；商宇航等^[24]推导了双曲线密度模型的波数域重力异常表达式。

空间域方法数值模拟精度高，网格剖分灵活，但计算大量位置点的异常场数据时，速度较慢。相对空间域方法，波数域方法速度较快，但精度相对较低，且垂向网格等间隔采样不利于模拟复杂地形或复杂模型。为了兼顾数值模拟的精度与效率，充分利用空间域方法和波数域方法的优势，本文针对任意形状、任意密度分布的二度体重力异常计算效率低的问题，提出一种空间—波数混合域二度体重力异常正演方法。该方法对空间域引力位满足的二维偏微分方程沿水平方向进行一维傅里叶变换，把空间域引力位满足的二维偏微分方程转化为不同波数相互独立的一维常微分方程，将一个复杂问题分解为多个小问题。该方法实现了不同波数之间常微分方程相互独立，具有高度并行性。该方法在垂向上为空间域，便于浅层网格剖分适当加密，深层网格剖分适当稀疏，有利于适应复杂地形的模拟，兼顾了计算精度与计算效率；采用有限单元法求解不同波数满足的常微分方程，并充分利用追赶法求解定带宽线性方程组的高效性，进一步提高了计算效率。

1 理论方法 1.1 控制方程

引力位满足二维泊松方程^[25]

$ {\nabla ^2}U\left( {x,z} \right) = - 4{\rm{ \mathsf{ π} }}G\Delta \tilde \rho \left( {x,z} \right) $

(1)

式中：U为引力位；G为万有引力常数；Δρ为异常体剩余密度；(x, z)表示空间任意一点坐标。

对式(1)沿x方向进行一维傅里叶变换，得到

$ \frac{{{\partial ^2}\tilde U\left( {k,z} \right)}}{{\partial {z^2}}} - {k^2}\tilde U\left( {k,z} \right) = - 4{\rm{ \mathsf{ π} }}G\Delta \tilde \rho \left( {k,z} \right) $

(2)

式中：k为波数；$\widetilde U$为空间—波数混合域引力位；Δ$\widetilde \rho $为空间—波数混合域异常体剩余密度。式(2)即为空间—波数混合域引力位常微分方程。

特别地，当k=0时，空间—波数混合域引力位满足常微分方程

$ \frac{{{\partial ^2}\tilde U\left( {0,z} \right)}}{{\partial {z^2}}} = - 4{\rm{ \mathsf{ π} }}G\Delta \tilde \rho \left( {0,z} \right) $

(3)

1.2 边界条件

在无源区域，式(2)的通解为

$ \tilde U = A{{\rm{e}}^{kz}} + B{{\rm{e}}^{ - kz}} $

(4)

式中A和B为任意常数。在笛卡尔坐标系中，令坐标轴z向下为正，模型的上边界z=z_min，下边界z=z_max，如图 1所示。则上、下边界条件为

$ \left\{ {\begin{array}{*{20}{l}} {{{\left. {\frac{{\partial \tilde U}}{{\partial z}}} \right|}_{z = {z_{\min }}}} = k\tilde U}\\ {{{\left. {\frac{{\partial \tilde U}}{{\partial z}}} \right|}_{z = {z_{\max }}}} = - k\tilde U} \end{array}} \right. $

(5)

1.3 重力场和重力梯度张量的计算

在无源区域，可以通过频率域求导得到重力场和重力梯度张量。空间域重力场g、重力梯度张量T与引力位分别为

$ \mathit{\boldsymbol{g}} = \left[ {\begin{array}{*{20}{l}} {{g_x}}\\ {{g_z}} \end{array}} \right] = - \left[ {\begin{array}{*{20}{l}} {\frac{{\partial U}}{{\partial x}}}\\ {\frac{{\partial U}}{{\partial z}}} \end{array}} \right] $

(6)

$ \mathit{\boldsymbol{T}} = \left[ {\begin{array}{*{20}{l}} {{g_{xx}}}&{{g_{xz}}}\\ {{g_{zx}}}&{{g_{zz}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{l}} {\frac{{\partial {g_x}}}{{\partial x}}}&{\frac{{\partial {g_x}}}{{\partial z}}}\\ {\frac{{\partial {g_z}}}{{\partial x}}}&{\frac{{\partial {g_z}}}{{\partial z}}} \end{array}} \right] $

(7)

由式(6)可得空间—波数混合域引力位与重力场的关系式

$ \mathit{\boldsymbol{\tilde g}} = \left[ {\begin{array}{*{20}{c}} {{{\tilde g}_x}}\\ {{{\tilde g}_z}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} { - {\rm{i}}k\tilde U}\\ {k\tilde U} \end{array}} \right] $

(8)

式中$\mathit{\boldsymbol{\widetilde g}}$为空间波数混合域重力场。

由式(7)可得空间—波数混合域引力位与重力梯度张量的关系式

$ \mathit{\boldsymbol{\tilde T}} = \left[ {\begin{array}{*{20}{l}} {{{\tilde g}_{xx}}}&{{{\tilde g}_{xz}}}\\ {{{\tilde g}_{zx}}}&{{{\tilde g}_{zz}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{k^2}\tilde U}&{k\left( { - {\rm{i}}k\tilde U} \right)}\\ { - {\rm{i}}k\left( {k\tilde U} \right)}&{ - {k^2}\tilde U} \end{array}} \right] $

(9)

式中$\mathit{\boldsymbol{\widetilde T}}$为空间—波数混合域重力梯度张量。

在有源区域，重力场和重力梯度张量，即异常体内部的$\frac{{\partial \mathit{\widetilde U}}}{{\partial z}}$可以通过形函数求导^[25]得到。

2 常微分方程求解

式(2)为引力位在空间—波数混合域满足的常微分方程，本文采用基于二次插值的一维有限单元法对其进行数值模拟。该方法一方面能实现复杂模型和复杂地形的高精度模拟，另一方面在一维条件下利用追赶法可实现对角线性方程组(五对角)的高效、高精度求解。

联立式(2)与式(5)可得空间—波数混合域引力位的边值问题

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{{\partial ^2}\tilde U}}{{\partial {z^2}}} - {k^2}\tilde U = - 4{\rm{ \mathsf{ π} }}G\Delta \tilde \rho \left( {k,z} \right)}\\ {{{\left. {\frac{{\partial \tilde U}}{{\partial z}}} \right|}_{z = {z_{\min }}}} = k\tilde U}\\ {{{\left. {\frac{{\partial \tilde U}}{{\partial z}}} \right|}_{z = {z_{\max }}}} = - k\tilde U} \end{array}} \right. $

(10)

基于变分原理^[25]构造泛函，可得到与边值问题(式(10))等价的变分问题

$ \left\{ \begin{array}{l} F\left( {\tilde U} \right) = \int_{{z_{\min }}}^{{z_{\max }}} {\left[ {{{\left( {\frac{{\partial \tilde U}}{{\partial z}}} \right)}^2} + {{\left( {k\tilde U} \right)}^2} + 8{\rm{ \mathsf{ π} }}G\Delta \tilde \rho k\tilde U} \right]{\rm{d}}z} + \\ \;\;\;\;\;\;\;\;\;\;\;\;k\left( {k\tilde U_{{z_{\max }}}^2 + k\tilde U_{{z_{\min }}}^2} \right)\\ {\rm{ { δ} }}F\left( {k\tilde U} \right) = 0 \end{array} \right. $

(11)

在如图 1所示的笛卡尔坐标系中，沿z方向进行单元剖分。在每个单元内，采用二次插值函数，即空间—波数混合域引力位${\mathit{\widetilde U}}$和空间—波数混合域剩余密度Δ$\widetilde \rho $在每个单元内是二次变化的。对变分问题(式(11))逐项进行单元分析、总体合成(见附录A)，得到扩展矩阵或列阵

$ F\left( {\tilde U} \right) = {\mathit{\boldsymbol{u}}^{\rm{T}}}\mathit{\boldsymbol{Ku}} - 2{\mathit{\boldsymbol{u}}^{\rm{T}}}\mathit{\boldsymbol{P}} $

(12)

式中：u为由某个波数在z方向各个节点的引力位${\mathit{\widetilde U}}$组成的列向量；K为最后组合成的五对角矩阵；P是源项；$\boldsymbol{K}=\sum \boldsymbol{K}_{1 q}+\boldsymbol{B}_{1}+\boldsymbol{B}_{2}, \boldsymbol{P}=\sum-\boldsymbol{P}_{1 q} \boldsymbol{s z}_{q}$，其中单元积分矩阵K_1q、B₁、B₂、P_1q和列阵sz_q的求解见附录A。

对F(${\mathit{\widetilde U}}$)求变分，得

$ {\rm{ { δ} }}F\left( {\tilde U} \right) = {\rm{ { δ} }}{\mathit{\boldsymbol{u}}^{\rm{T}}}\left( {2\mathit{\boldsymbol{Ku}} - 2\mathit{\boldsymbol{P}}} \right) = 0 $

(13)

由于δu≠0，故

$ \mathit{\boldsymbol{Ku}} = \mathit{\boldsymbol{P}} $

(14)

对该五对角线性方程组(式(14))的求解参见文献[26]。通过重力场、重力梯度张量与引力位之间的关系式(式(6)和式(7))，得到空间—波数混合域的重力场和重力梯度张量，采用Gauss-FFT对空间波数混合域引力位、场及梯度张量进行反变换，可得到空间域的引力位、重力场和重力梯度张量。

3 数值试验

为验证本文算法的正确性及可靠性，分别设计了截面为矩形的常密度和变密度二度体模型及任意密度分布的Teplow模型^[13]。最后，对本文算法的计算效率随网格剖分大小的变化规律进行了统计分析。以下算法均利用Fortran95语言编程串行计算实现，笔记本电脑配置为：CPU-Inter Core i4-4790，主频为4.0GHz，内存为32GB。

3.1 常密度二度体模型

设计如图 2所示的截面为矩形的常密度二度体模型。研究区域为：x方向[-500m，500m]，z方向[0，500m]。网格数为200×100，水平和垂向采样间隔均为5m。异常体范围沿x方向为[-100m，100m]，z方向为[200m，300m]，异常体剩余密度为100kg/m³。重力异常及梯度张量的解析解^[27]和数值解及与地面位置各观测点的相对误差分别如图 3、图 4所示。

从图 3可以看出，重力场x分量的数值解与解析解吻合度较高，地面位置的最大相对误差小于1.0×10^-4。从图 4可以看出，重力梯度张量的数值解与解析解吻合较好，地面各观测点的重力梯度z分量在零值点附近相对误差最大，约为2.0×10^-3，这是零值点附近数值计算误差大引起的，在其他位置相对误差均小于1.0×10^-3。综合图 3和图 4可以看出，重力场和重力梯度张量的相对误差较小，验证了本文算法的正确性和高精度。

3.2 变密度二度体模型

设计一个截面为矩形的变密度二度体模型，研究区域为：x方向[-500m，500m]，z方向[0, 500m]。网格数为200×100，水平和垂向采样间隔均为5m。异常体沿x方向范围为[-100m，100m]，z方向为[150m，300m]。异常体的密度随深度的关系为

$ \rho \left( z \right) = 1.54 + 0.24z - 0.035{z^2}\;\;\;\;\;150 \le z \le 300 $

当z=0时(即地面)，重力异常和梯度张量解析解^[27]及相对误差如图 5所示。可以看出，数值解与解析解形态十分吻合，重力异常和重力梯度张量的相对误差均小于2.0×10^-4，与常密度模型计算结果相比，变密度模型的计算结果精度更高。说明本文算法更适用于变密度二度体模型的重力异常数值模拟。

3.3 Teplow密度分布模型

为了进一步检验本文算法对任意形状截面、任意密度分布情况下二度体重力异常计算的适应性，引入Reeder等^[13]的Teplow密度分布模型(图 6)。研究区域范围为：x方向[0，4268.3m]，z方向为[0，1829.3m]，剖分网格数为994×743。

采用本文算法计算地面位置的重力异常，结果如图 7所示。由图可知，传统空间域累加算法^[27]与本文算法的计算结果吻合很好，表明本文算法能够计算复杂密度分布情况下的重力异常，且精度较高。

3.4 计算效率统计

计算效率是评价数值模拟方法好坏的一个重要指标。图 8给出了本文方法的计算时间随网格剖分数目变化的曲线。可以看出，计算时间随着网格剖分数目的大小呈近似线性增长。当网格剖分数目为500×500时，即251001个节点，采用本文算法计算整个剖面耗时约0.24s，采用传统空间域累加算法计算地面501个节点需耗时38.73s^[27]，表明本文算法计算效率更高。

4 结论

本文提出了一种高效、高精度的空间—波数混合域二度体重力异常正演方法。设计了二度体模型开展数值试验，通过对比数值解与解析解验证了本文方法的正确性和可靠性。通过计算Teplow密度分布模型重力异常，表明了该算法对任意密度分布的二度体模型具有很好的适应性。数值算例结果表明，本文算法具有较高的计算精度和计算效率，为计算任意复杂形体的重力异常提供了一种新途径，对提高二度体重力异常反演成像的效率具有现实意义。

附录A

求解文中变分问题(式(11))时，需将整个区域的单元积分分解为各个单元的积分之和，则式(11)变为

$ \begin{array}{l} F\left( {\tilde U} \right) = \sum {\int_q {{{\left( {\frac{{\partial \left( {\tilde U} \right)}}{{\partial z}}} \right)}^2}} {\rm{d}}z} + \sum {\int_q {{{\left( {k\tilde U} \right)}^2}} {\rm{d}}z} + \\ \;\;\;\;\;2\sum {\int_q {\tilde J\tilde U{\rm{d}}z} } + k\tilde U_{{z_{\max }}}^2 + k\tilde U_{{z_{\min }}}^2 \end{array} $

(A-1)

其中

$ \tilde J\left( {k,z} \right) = 4{\rm{ \mathsf{ π} }}G\Delta \tilde \rho $

(A-2)

令

$ \left. {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{N}} = {{\left( {{N_j},{N_p},{N_m}} \right)}^{\rm{T}}}}\\ {{\mathit{\boldsymbol{u}}_q} = {{\left( {{u_j},{u_p},{u_m}} \right)}^{\rm{T}}}} \end{array}} \right\} $

(A-3)

式中：j、m分别为单元两端节点坐标；p为单元中点坐标。有

$ \mathit{\boldsymbol{u}} = {\mathit{\boldsymbol{N}}^{\rm{T}}}{\mathit{\boldsymbol{u}}_q} = \mathit{\boldsymbol{u}}_q^{\rm{T}}\mathit{\boldsymbol{N}} $

(A-4)

式(A-1)中第一项单元积分为

$ \begin{array}{*{20}{c}} {\int_q {\frac{{\partial \mathit{\boldsymbol{u}}}}{{\partial z}}\frac{{\partial \mathit{\boldsymbol{u}}}}{{\partial z}}{\rm{d}}z} = \mathit{\boldsymbol{u}}_q^{\rm{T}}\left( {\int_q {\frac{{\partial \mathit{\boldsymbol{N}}}}{{\partial z}}\frac{{\partial {\mathit{\boldsymbol{N}}^{\rm{T}}}}}{{\partial z}}{\rm{d}}z} } \right){\mathit{\boldsymbol{u}}_q}}\\ { = \mathit{\boldsymbol{u}}_q^{\rm{T}}{\mathit{\boldsymbol{K}}_{1q}}{\mathit{\boldsymbol{u}}_q}} \end{array} $

(A-5)

其中

$ {\mathit{\boldsymbol{K}}_{1q}} = \int_q {\frac{{\partial \mathit{\boldsymbol{N}}}}{{\partial z}}\frac{{\partial {\mathit{\boldsymbol{N}}^{\rm{T}}}}}{{\partial z}}{\rm{d}}z} = \frac{1}{{3l}}\left[ {\begin{array}{*{20}{c}} 7&{ - 8}&1\\ { - 8}&{16}&{ - 8}\\ 1&{ - 8}&7 \end{array}} \right] $

(A-6)

式中l为单元长度。

式(A-1)中第二项单元积分为

$ \int_q {{k^2}{\mathit{\boldsymbol{u}}^2}{\rm{d}}z} = \mathit{\boldsymbol{u}}_q^{\rm{T}}\left( {\int_q {{k^2}\mathit{\boldsymbol{N}}{\mathit{\boldsymbol{N}}^{\rm{T}}}{\rm{d}}z} } \right){\mathit{\boldsymbol{u}}_q} = \mathit{\boldsymbol{u}}_q^{\rm{T}}{\mathit{\boldsymbol{K}}_{2q}}{\mathit{\boldsymbol{u}}_q} $

(A-7)

其中

$ {\mathit{\boldsymbol{K}}_{2q}} = \int_q {{k^2}\mathit{\boldsymbol{N}}{\mathit{\boldsymbol{N}}^{\rm{T}}}{\rm{d}}z} = \frac{{l{k^2}}}{{30}}\left[ {\begin{array}{*{20}{c}} 4&2&{ - 1}\\ 2&{16}&2\\ { - 1}&2&4 \end{array}} \right] $

(A-8)

式(A-1)中第三项单元积分为

$ 2\int_q {\tilde J\mathit{\boldsymbol{u}}{\rm{d}}z} = 2\int_q {{j_{az}}u{\rm{d}}z} $

(A-9)

利用形函数将j_az表示为

$ {j_{az}} = {j_{azj}}{N_j} + {j_{azp}}{N_p} + {j_{azm}}{N_m} = {\mathit{\boldsymbol{N}}^{\rm{T}}}\mathit{\boldsymbol{s}}{\mathit{\boldsymbol{z}}_q} = \mathit{\boldsymbol{sz}}_q^{\rm{T}}\mathit{\boldsymbol{N}} $

(A-10)

其中

$ \mathit{\boldsymbol{s}}{\mathit{\boldsymbol{z}}_q} = {\left( {{j_{azj}},{j_{azp}},{j_{azm}}} \right)^{\rm{T}}} $

(A-11)

则式(A-9)中的单元积分为

$ \int_q {{j_{az}}\mathit{\boldsymbol{u}}{\rm{d}}z} = \mathit{\boldsymbol{u}}_q^{\rm{T}}\left( {\int_q {\mathit{\boldsymbol{N}}{\mathit{\boldsymbol{N}}^{\rm{T}}}{\rm{d}}z} } \right)\mathit{\boldsymbol{s}}{\mathit{\boldsymbol{z}}_q} = \mathit{\boldsymbol{u}}_q^{\rm{T}}{\mathit{\boldsymbol{P}}_{1q}}\mathit{\boldsymbol{s}}{\mathit{\boldsymbol{z}}_q} $

(A-12)

其中

$ {\mathit{\boldsymbol{P}}_{1q}} = \int_q {\mathit{\boldsymbol{N}}{\mathit{\boldsymbol{N}}^{\rm{T}}}{\rm{d}}z} = \frac{l}{{30}}\left( {\begin{array}{*{20}{c}} 4&2&{ - 1}\\ 2&{16}&2\\ { - 1}&2&4 \end{array}} \right) $

(A-13)

式(A-1)中z=z_min、z=z_max分别为第一个和最后一个节点，可将其分别扩展成

$ \left\{ {\begin{array}{*{20}{l}} {ku_{{z_{\min }}}^2 = {\mathit{\boldsymbol{u}}^{\rm{T}}}{\mathit{\boldsymbol{B}}_1}\mathit{\boldsymbol{u}}}\\ {ku_{{z_{\max }}}^2 = {\mathit{\boldsymbol{u}}^{\rm{T}}}{\mathit{\boldsymbol{B}}_2}\mathit{\boldsymbol{u}}} \end{array}} \right. $

(A-14)

其中

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{B}}_1} = \left[ {\begin{array}{*{20}{c}} k&{}&{}&{}\\ {}&0&{}&{}\\ {}&{}& \ddots &{}\\ {}&{}&{}&0 \end{array}} \right]\\ {\mathit{\boldsymbol{B}}_2} = \left[ {\begin{array}{*{20}{c}} 0&{}&{}&{}\\ {}& \ddots &{}&{}\\ {}&{}&0&{}\\ {}&{}&{}&k \end{array}} \right] \end{array} \right. $

(A-15)

综上，可将式(A-1)表示为

$ F\left( \mathit{\boldsymbol{u}} \right) = {\mathit{\boldsymbol{u}}^{\rm{T}}}\mathit{\boldsymbol{Ku}} - 2{\mathit{\boldsymbol{u}}^{\rm{T}}}\mathit{\boldsymbol{P}} $

(A-16)