0 引言

随着高精度磁异常资料解释工作的迅速发展，磁异常剖面反演已经难以适应精细地质构造的刻画，所以三维反演和解释成为一种重要的地球物理方法，并且得到了广泛应用 (Li and Oldenburg, 1996; Pilkington, 1997; Portniaguine and Zhdanov, 2002; Fregoso and Gallardo, 2009; Shamsipour et al., 2011; 陈召曦等，2012；Luo et al., 2015)，这些文献均采用长方体正演公式进行磁异常计算.对于真实的地质体，其空间结构往往较为复杂，例如断层、褶皱、岩脉等，如果采用长方体进行模拟，为达到逼近精度，需要缩小长方体的尺寸，这样就增加了计算量.为此，有必要对均匀磁化多面体的磁异常进行研究.

关于均匀磁化多面体的磁异常前人已经取得了许多成果.Barnett (1976)提出将磁异常体积分转换为多边形面积分，再将多边形面积分转换为若干个三角形面积分的叠加，以得到解析解.Coggon (1976)的方法是将多面体的每一个面等价为一个无限延伸的棱柱，其轴线垂直于相应的多边形，且磁化强度方向平行于轴线，最后将所有多边形的磁效应叠加求和.Singh和Guptasarma (2001)在不需要坐标变换的条件下通过定义一种新的矢量将多面体的面积分转换为线积分，然后将线积分叠加求和.Okabe (1979)提出了一种将体积分转变为面积分的任意多面体重磁位场正演公式，为计算方便，在计算每一个面积分时首先进行坐标变换使积分面与Z坐标轴垂直，待积分完成后再进行坐标反变换，最后由位场的叠加原理，将所有积分面产生的磁异常进行求和.王邦华等 (1980)、贾春生 (1999)、骆遥和姚长利 (2007)等基于Okabe的方法分别提出了改进形式.上述方法均没有解决观测点位于多面体外延面时存在的奇点问题，而这一问题在磁异常正反演中非常重要.本文重点讨论了多面体磁异常的奇点问题，并给出了无奇点的解析表达式.最后通过模型试验说明方法的正确性和应用前景.

1 磁荷面积分基本方程

根据磁荷理论，可以将均匀磁化多面体的体积分转换为面积分，然后对所有面积分求和，也即：

(1)

式中，F_k是多面体的磁场的某个分量，F_kj是第j个磁荷面所产生的该分量的磁场值.

首先讨论与Z轴垂直的水平梯形的磁荷面的积分解.根据磁荷理论，水平磁荷面磁场表达式为

(2a)

(2b)

(2c)

其中σ为磁荷面密度；x，y，z为观测点坐标值；ξ，η，ζ₁为磁荷面上任意点的的坐标值.坐标系选为Z轴垂直向下的右手直角坐标系.上述磁场的计算可以归结为以下两类积分的解，公式为

(3a)

(3b)

对于图 1所示的水平梯形，(3a) 式可写为

(4)

其中，a₁=a′₁，b₁=b′₁，a₂=a′₂-a′₁cotθ₁，b₂=b′₂-a′₁cotθ₂，a₃=cotθ₁，b₃=cotθ₂.其积分结果见附录A中的公式A1.

公式 (3b) 在图 1情形下可改写为

(5)

a₁，b₁，a₂，b₂，a₃，b₃的意义与f(μ, v, ω) 式同.其积分结果见附录A中的公式A2.

2 新坐标系的建立

多面体的各个面与笛卡尔坐标系存在各种空间关系，为了使用上述推导的水平梯形磁荷面磁场的解析解，必须将多边形转化为与Z轴垂直，为此需要建立新的坐标系.设多面体第j个多边形平面ΔS_j外表面顶点A₁、A₂、A₃、…、A_N按逆时针方向顺序排列，以矢量A₁A_N为新坐标系的X′轴方向，多边形平面ΔS′_j内法线-n_j(A₁A_N × A₁A₂) 为Z′轴方向，矢量 (A₁A_N × A₁A₂)×A₁A_N为Y′轴方向 (图 2)，得到计算所需要的新坐标系.

如果新坐标系X′、Y′、Z′轴与原坐标系X、Y、Z轴之间的夹角依次为α₁、β₁、γ₁；α₂、β₂、γ₃；α₃、β₃、γ₃，则将原坐标系中矢量或空间中任意一点变换到新坐标系下，有以下坐标变换，公式为

(6)

通过该变换，可以将原坐标系中的任意观测点P(x, y, z) 和任意多边形顶点A(ξ_i, η_i, ζ_i) 转换为新坐标系中的观测点P′(x′, y′, z′) 和多边形顶点A′(ξ′_i, η′_i, ζ′_i).

3 多边形水平磁荷面的磁场

为方便起见，先对任意四边形磁荷面的磁场进行计算，再进行推广.

假设任意四边形A₁A₂A₃A₄(见图 3a)，各角点坐标值为A₁(ξ′₁, η′₁, ζ′₁)，A₂(ξ′₂, η′₂, ζ′₂)，A₃(ξ′₃, η′₃, ζ′₃) 和A₄(ξ′₄, η′₄, ζ′₄)，且θ₁₂，θ₂₃，θ₃₄和θ₄₁分别表示A₁A₂，A₂A₃，A₃A₄和A₄A₁与x′轴正向的夹角 (0≤θ_{i, i+1}≤π).

经推导 (见附录B)，N变形水平磁荷面磁场H′_ax表达式为 (王邦华等，1980；骆遥和姚长利，2007)：

(7)

其中：

J为磁化强度矢量的模，I和A分别表示磁倾角和磁方位角.

以类似的方式，可以得到水平N边形磁荷面的磁场H′_ay和Z′_a分量.在计算H′_ay时，需要将多边形划分为许多上下底边与y′轴相平行的梯形 (见图 3b)，其表达式为 (王邦华等，1980；骆遥和姚长利，2007)：

(8)

(9)

由 (7) 至 (9) 式，即可求出水平磁荷面在三个方向上的磁场分量.

4 关于奇点问题的讨论

水平磁荷面在三个方向上的磁场分量表达式前人已经进行了深入研究，例如Okabe (1979)、王邦华等 (1980)、贾春生 (1999)、骆遥和姚长利 (2007)等.但是这些研究都没有指出磁场表达式中可能存在的让表达式无意义的奇点问题.本文将对此进行分析，并得到无解析奇点的理论表达式.

首先考察公式 (7)，若对数项内部的真数为零，例如R_i+r_i为零，则：

将上式第二项及第三项移到等式右边，然后两边平方，可得：

合并同类项，整理上式，可得：

最终，可以得到R_i+r_i为零的条件 (记为Cond1) 为

(10)

由上式中的条件1和3可知，观测点 (x′, y′, z′) 其实是在线段A_iA_i+1的延长线上 (不包括A_iA_i+1线段)，设观测点 (x′, y′, z′) 在线段A₁A₂的延长线上，即：

(11)

若R_i+1+r_i+1=0，则有：

(12)

将上式中的条件1和条件3带入附录B1中第一项，则其积分H′_axk为

上式中的第一项记为H′_axk1，其积分表达式为

H′_axk1是由A₁A′₂段积分得到的结果，另外再加上A′₂A₄段的积分结果，就可以得到A₁A₄边的积分表达式.而将 (11) 式中条件1带入H′_axk中第二项，其结果记为H′_axk2，可得：

因为观测点 (x′, y′, z′) 在线段A₁A₂的延长线上 (不包括A₁A₂线段)，于是有 (η′₁-y′)(η′₂-y′)≥0恒成立.

当 (η′₁-y′)>0且 (η′₂-y′)>0时:

当 (η′₁-y′) < 0且 (η′₂-y′) < 0时:

当 (η′₁-y′)(η′₂-y′)=0时，由于测点 (x′, y′, z′) 位于线段A₁A₂的延长线上，所以必有 (y′-η′₁)=0且 (y′-η′₂)=0，从而得到η′₁=η′₂，也即直线A₁A₂平行于x′坐标轴，从而三角形A₁A₂A′₂所围面积为零，所以必有H′_axk=0.

由以上分析，可以得到包含公式 (10) 情况 (记为Cond1) 的无解析奇点水平多边形的磁场分量H′_ax的表达式为

(13)

其中：

下面讨论水平多边形的磁场分量H′_ay的表达式，使得H′_ay无意义的奇点条件依然是R_i+r_i为零，也即满足公式 (10).此时，将多边形划分为许多上下底边与y′轴相平行的梯形 (见图 3b)，并设测点 (x′, y′, z′) 位于线段A₁A₂的延长线上，于是有：

将上式中的第一项记为H′_ayk1，类似于H′_ax1的推导方法，可得其积分表达式为

H′_ayk1是由A₂A′₁段积分得到的结果，另外再加上A′₁A₃段的积分结果，就可以得到A₂A₃段完整的积分值.

而将 (12) 式中条件1带入上式第二项，其结果记为H′_ayk2，可得：

当 (ξ′₂-x′)>0且 (ξ′₁-x′)>0时：

当 (ξ′₂-x′) < 0且 (ξ′₁-x′) < 0时：

当 (ξ′₁-x′)(ξ′₂-x′)=0时，由于测点 (x′, y′, z′) 位于线段A₁A₂的延长线上，所以必有 (ξ′₁-x′)=0且 (ξ′₂-x′)=0，从而得到ξ′₁=ξ′₂，也即直线A₁A₂平行于y′坐标轴，从而三角形A₁A₂A′₁所围面积为零，所以必有H′_ayk=0.

由以上分析，可以得到包含公式 (10) 情况 (记为Cond1) 的无解析奇点水平多边形的磁场分量H′_ay的表达式为

(14)

其中：

最后讨论水平多边形的磁场分量Z′_a的表达式，使得Z′_a无意义的奇点条件依然是 (ζ′_i-z′) 为零或sinθ_{i, i+1}为零，由积分公式 (2c) 可知，当 (ζ′_i-z′) 为零时，被积函数为零，所以其积分值为零.从宏观场论也可知，当观测点位于磁荷面的延伸面上时，磁场在垂直于磁荷面的方向上分量为零.

当sinθ_{i, i+1}为零时，例如sinθ_{1, 2}为零，直线A₁A₂平行于x′轴，过点A₂所作平行于x′轴的直线与直线A₁A₂相交于A′₂，则三角形A₁A₂A′₂的面积为零，因此面积分公式 (2c) 的积分值为零.综上所述，有：

(15)

其中：

由 (13)、(14)、(15) 三式组成新的无解析奇点的水平磁荷面的三个方向的磁场分量.

设沿原坐标系三个轴向磁场分量为H_ax、H_ay和Z_a，根据坐标变换关系，可以将新坐标系下的H′_ax、H′_ay和Z′_a转换到原坐标系，公式为

(16)

利用如下公式，可以得到任意产状的磁荷面的总磁场ΔT为

(17)

最后将各磁荷面的磁场线性叠加，即可得到整个均匀磁化多面体的磁场.

5 模型试验

为了检验方法的正确性，将本文推导的 (13) 至 (15) 式组成的无解析奇点的多面体正演公式与王邦华等 (1980)、骆遥和姚长利 (2007)的公式进行对比.具体方法是，将本文方法与其他两种方法同时编写程序，并应用于相同的模型体及观测网格.本文模型试验观测的平面范围均为20000 m×20000 m，测点网格大小均为100 m×100 m，观测平面的高度为Z=0.模型体均匀磁化，磁化倾角60°，磁化偏角10°.

图 4a为标准四棱台模型，上顶面为A′B′D′C′，下底面为ABDC，均为垂直于Z轴的矩形，上下两个矩形对应边平行，且中心点在XY平面上的投影重合，模型试验中根据图 4a设

计了三个模型，各模型所对应的顶点坐标分别见表 1、表 2和表 3，分别被命名为模型一、模型二和模型三.

模型一对应的角点坐标见表 1，共采用四种方法对模型一的磁场进行正演，分别为王邦华等方法 (1980)(Wang法)，骆遥和姚长利方法 (2007)(Luo法)，水平剖分近似薄长方体方法 (Horizontal thin cuboid方法，简称HTC法)，以及本文的无解析奇点方法 (No analytical singularity方法，简称NAS法).图 5a、5b、5c和5d分别为Wang法、Luo法、HTC法和NAS法正演结果的等值线图，单位为nT.Wang法和Luo法均存在奇点，图中用黑色点标出，奇点坐标见表 1.图 5a、5b除了奇点处没有正演值外，其余各观测点处的正演值与图 5d完全一致.而HTC法和NAS法不存在奇点，且这两种方法在Wang法和Luo法的奇点位置处的正演值一致，可达0.0001 nT的精度.但由于HTC法需要将模型剖分成薄长方体，计算量很大，模型试验表明，若使HTC法与NAS法的均方误差在0.001 nT以内，则HTC法计算机耗时约为NAS法的十倍以上.这表明本文推导的NAS法解决了奇点问题，并且相对于HTC法计算速度快.

模型二对应的角点坐标见表 2，相对于模型一，模型二的顶面和底面分别上移1000m和下移1000m.图 6a、6b、6c和6d分别对应Wang法、Luo法、HTC法和NAS法的正演等值线图.Wang法和Luo法均存在奇点，图中用黑色点标出，奇点坐标见表 2.表 1和表 2的奇点坐标不同，这说明随着四棱台顶底面的上下移动，奇点位置也随之改变.其余结论与模型一的结论相同.

模型三的角点坐标见表 3，相对于模型一，模型三的顶面四条边和底面四条边分别外扩1000m，图 7a、7b、7c和7d分别对应Wang法、Luo法、HTC法和NAS法的正演等值线图.Wang法和Luo法均存在奇点，图中用黑色点标出，奇点坐标见表 3.表 1和表 3的奇点坐标不同，这说明四棱台顶底面的大小不同，奇点位置也相应改变.其余结论与模型一的结论相同.

图 8a为垂直四棱柱模型，四条侧边平行于Z轴，顶面与底面分别由两个三角形组成，且三角形A′B′C′与三角形B′C′D′不在同一平面上，三角形ABC与三角形BCD不在同一平面上，此模型记为模型四，其角点坐标见表 4.分别采用Wang法和Luo法 (见图 9a)、本文的NAS法 (见图 9b)、垂直剖分近似直立长方体方法 (Vertical slim cuboid方法，简称VSC法，见图 9c、9e) 四种方法进行正演模拟.其中垂直剖分和近似方式见图 8，作四条平行直线与顶面交于s₁、s₂、s₃、s₄，与底面交于s₅、s₆、s₇、s₈，其中s₁、s₂、s₃、s₄及s₅、s₆、s₇、s₈在XY平面上投影的顺序连线均为平行于坐标轴的正方形.取s₁、s₂、s₃、s₄四个点纵坐标的平均值为图 8b中的h₁，取s₅、s₆、s₇、s₈四个点纵坐标的平均值为图 8b中的h₂.分别对图 8a进行水平方向上2×2剖分和4×4剖分，即分别剖分成4个和16个形如图 8b的小立方体，其对应的VSC法正演的等值线图分别为图 9c和图 9e.Wang法和Luo法均存在奇点 (见图 9a)，图中用黑色点标出，奇点坐标见表 4.图 9d、9f分别为2×2剖分和4×4剖分的VSC法与NAS法正演差值的等值线图.由图 9及表 5可知，VSC法不存在奇点，但相对于NAS法，其计算时间较长，剖分越细，VSC法结果越接近NAS法的结果，但计算机耗时越长 (见表 5，计算机参数为Intel^® Core (TM) 2CPU 2.4GHz 2.4GHz，4.00GB内存，64位操作系统).

上述模型试验说明，模型几何参数的变化会导致奇点位置的改变，如果将包含奇点的多面体磁异常理论用于三维磁界面反演，则由于奇点位置的不确定性而导致反演过程无法顺利进行.而如果采用横向的或纵向的立方体剖分近似方法，则具有耗时较多的缺点，给使用造成不便.

6 结论与展望

本文针对已有的多面体磁异常表达式存在的解析奇点问题，讨论了奇点存在的条件，并根据磁荷面积分理论推导了无解析奇点的多面体磁异常表达式.模型试验表明，随着模型几何参数和位置的变化，已有的多面体磁异常的奇点位置也随之改变.本文推导的多面体磁异常表达式没有解析奇点，正演结果表明了其结果的正确性.横向剖分和纵向剖分的长方体近似方法正演的结果与多面体无解析奇点表达式的正演结果近似，且近似程度与模型剖分精细度程度正相关，但随着剖分程度的提高，计算机耗时增加较多.在实际的三维界面反演中，如果采用含奇点的多面体磁异常表达式，由于奇点位置在每次迭代过程中不断改变而让反演无法进行.如果采用立方体近似复杂曲面，则模型剖分程度要求较高，从而增加计算机耗时.因此，本文的多面体无解析奇点表达式在三维磁界面反演中具有重要作用.

附录A  f(u, v, w) 与φ(u, v, w) 的积分结果

文中公式 (4) 的积分结果为

(A1)

文中公式 (5) 的积分结果为

(A2)

其中：

附录B   H′_ax分量推导过程

如图 3a所示，计算H′_ax分量时，过A₂、A₄作x′的两条平行线，分别与对边相交于A′₂(ξ″₂, η′₂, ζ′₂)，A′₄(ξ″₄, η′₄, ζ′₄) 两点，于是四边形被划分为顶边与x′轴平行的三个梯形 (三角形是梯形的特例).这样计算四边形磁荷面的磁场H′_ax的积分可表示为

(B1)

上式中三项积分均属于f(μ, ν, ω) 类型，令u=x′，v=y′，w=z′和τ₁=ζ′₁，并且将各积分限具体数值替换附录A1式中相应的a₁，a₂，a₃，b₁，b₂，b₃值 (例如在第一个积分式中，a₁=η′₁，b₁=η′₂，a₂=ξ′₁-η′₁cotθ_{1, 2}，b₂=ξ′₁-η′₁cotθ_{4, 1}，a₃=cotθ_{1, 2}，b₃=cotθ_{4, 1})，经过简化后可得：

(B2)

其中ζ′₁=ζ′₂=ξ′₃=ξ′₄.

致谢 感谢审稿专家提出的修改意见和编辑部的大力支持！