均匀任意形体正演方法

Answer 1

目前国内外用电子计算机作正演计算的方法很多，本节讨论计算密度或磁性均匀的地质体的重磁异常的方法。物性均匀对一般的地质地球物理问题都可以近似适用，而对于“非均匀”体也可以近似地分解成若干个“均匀”体的组合，且可根据精度要求选择划分的精细程度，因此均匀地质体异常的计算具有重要的意义。

从前面一套正演公式可见，正演方法归纳为计算一系列三重积分或面积分。对于任意形体要靠解析方法求出这些积分是困难的，所以采用数值解法求其近似解，根据近似方法的不同，大致可以分成以下几类。

（1）“点元”法：将一个任意形体按适当的方法划分为若干个规则几体形体（长方体、正方体），每一个均视作“点元”，先用解析方法求出各个点元的三重积分值，再累加求和即得整个形体的三重积分的近似值，近似程度取决于全部“点元”与该形体的吻合程度。

（2）“线元”法：用两组相互垂直的平行面把任意形体分割成很多棱柱体，每一棱柱体的作用，以位于其柱中心线的“线元”来代替，用解析法求出各“线元”的作用值，然后在垂直于“线元”的截面上作二重数值积分，即得到整个形体的三重积分近似值，其近似程度除了取决于全部棱柱体与该形体的吻合程度以外，还取决于所采用的数值积分方法。

（3）“面元”法：用一组相互平行的平面去分割任意形体，每个截面内用一个多边形去代替该形体在截面内的形状，用解析方法求出多边形域的二重积分值，然后在垂直于截面的方向上，用数值积分求出第三重积分，即得三重积分近似值。其近似程度取决于各多边形吻合该形体的各截面的程度及所采用的数值积分方法。

（4）表面积分法：根据（1.1.11）式，积分是在包围形体的全表面进行的。采用一系列多边形水平面的组合来近似全表面，用解析方法分别计算出每一个多边形水平面的积分值，然后累加求和。其近似程度取决于多边形水平面对该形体外表面的吻合程度。

综上所述，用电子计算机作重、磁正演计算的特点是：先分解，再合成，即分解成比较容易计算的若干个部分，然后再合成整个三重积分（或表面积分）的值。

1.2.1 “点元”法

“点元”法所取的各个点元的体积可以相同，也可不同。各个点元的物性可以相同，也可不同。通常是将勘探剖面之间的地质体用适当的长方体或立方体来近似，确定出各个点元的角点坐标，即可计算出该点元的三重积分值。

对于一个点元而言，其计算公式如下，设观测点 P 位于坐标原点，那么由式（1.1.3）和式（1.1.7）、（1.1.8）、（1.1.10）可得：

地球物理数据处理教程

式中：V₀、V₁、V₂、V₃、V₄、V₅、V₆分别为

地球物理数据处理教程

地球物理数据处理教程

图1.5 计算点元示意图

式中x₁、x₂、y₁、y₂、z₁、z₂分别为点元角点在x、y、z坐标轴方向上的投影坐标，即积分的上下限，如图1.5所示。

需要指出的是，当所选的“点元”各边不平行于坐标轴时，通常用坐标旋转变换来实现。

当给定了每个点元的位置，即给出了x₁、x₂、y₁、y₂、z₁、z₂及物性参数后，代入式（1.2.1）即可求出每个点元的三重积分值，然后累加求和，得到式（1.1.8）中的各值，用式（1.1.7）得到磁场的各个分量值，用式（1.1.3）得到重力值。

1.2.2 “线元”法

图1.6 “线元”法计算示意图

本方法是用一组平行于x轴的垂直面与另一组平行于y轴的垂直面去切割任意形体Q，于是得到一系列直立矩形棱柱体（图1.6），当分割异常体的这两组垂直截面很密时，所得到的这一系列直立棱柱体横截面很小，每个棱柱体可近似看作为直立棒状体而用解析方法计算出它在观测点P的作用值，然后对这一系列的作用值在垂直于棱柱体轴线的平面内进行二重数值积分，就可得到P点正演值。

在式（1.1.3）和（1.1.8）的三重积分中，把沿 z 方向的积分积出来（设棱柱体轴线是沿z轴方向），并令观测点P为坐标原点，于是就得到下列公式：

Δg=GσV₀

地球物理数据处理教程

ΔT=ΔXcosIcosD+ΔYcosIsinD+ΔZsinI

式中：

上式中R₁=（x²+y²+

）^1/2，R₂=（x²+y²+

）^1/2，z₁和z₂分别为棱柱体上顶面和下底面的z坐标值。

对上述的二重积分采用二维梯形公式或辛普生公式作数值积分，二维梯形公式或辛普生公式形式如下

地球物理数据处理教程

式中：F_k（x_i，y_j）为V₀，V₁…V₆中的被积函数；W_i和C_j分别为x和y方向的数值积分系数。

当观测点P（即坐标原点）与某一线元恰好位于同一铅垂线上时，该线元的x_i=y_j=0，此线元的F₁、F₂、F₄分别应为

地球物理数据处理教程

用“线元”法作实际正演计算时，一般用一水平矩形网格分割任意形体为一系列直立矩形棱柱体，矩形网格的大小视所要求的精度而定，给出每个棱柱体中心坐标x_i、y_j和上下顶底面的z₁和z₂坐标，并给出物性参数，计算出各线元的F_k值，用式（1.2.3）作数值积分得到正演值，数值积分是先沿y方向，再沿x方向进行。

1.2.3 “面元”法

本方法通常是用平等于yoz的一组平面去分割该任意形体为若干片，每一个片视为一个截面，每一个截面内用一个多边形去近似吻合形体在该截面内的形状，当适当地选取多边形的边数和形状时，就可以较为精确地表示出形体的截面形状。然后用解析的方法分别计算出每个截面内的多边形域的二重积分值，最后对这一组平行截面的积分值进行一维数值积分以求得正演值。

在实际计算时，我们将相互平行的勘探剖面作为这组截面，勘探所得到的地质体截面形状用多边形去吻合即可。所以该方法使用方便。其近似程度取决于多边形对地质体截面的吻合程度和所采用的数值积分方法。

图1.7 垂直于x轴的平面对地质体所截的截面

假设截面为垂直于x轴的yoz平面，如图1.7所示。并设对应于式（1.1.3）和式（1.1.8）中的三重积分V₀、V₁、V₂、V₃、V₄、V₅、V₆沿y和z方向的二重积分值为S₀、S₁、S₂、S₃、S₄、S₅、S₆，即是

地球物理数据处理教程

地球物理数据处理教程

我们下面就来具体求上述面积分的解析式。首先设某一截面与x轴交点P′的x值为x_α，即该截面内所有点的x坐标为x_α，用多边形KLMN……K来近似该截面内任意形体的截面形状，如图1.8所示。其中KL为多边形的第i条边，P_i为P′（x_α，0，0）到KL的垂直距离，g_i为KL在y轴上的截距，C_i是KL在Z轴上的截距，此外，距离R_i、R_i+1和角度α_i、β_i的意义见图1.8所示。

图1.8 x=x_α截面内用多边形近似的示意图

图1.9 计算多边形截面作用值的示意图

为了计算多边形域的二重积分，我们可以将多边形积分域分为各条边与y轴所围成的梯形的代数和，任意第i边与y轴所围成的梯形的积分值设为S_i，那么KL边与y轴围成的梯形为KK′L′L，如图1.9，它的作用值可以计算出来，首先计算S_i，从式（1.2.4）的第一式可得：

地球物理数据处理教程

先计算第一项（S_oi）₁，线段KL，可由两式的直线方程得

（y-y_i）/（z-z_i）=（y_i+1-y_i）/（z_i+1-z_i）

令

，

，则z=m_iy+c_i

代入（S_oi）₁中，有：

地球物理数据处理教程

式中：

，

而第二项（S_oi）₂有：

地球物理数据处理教程

同样可以求出多边形其他各边与y轴围成的梯形面积的积分值，然后求它们的代数和，我们规定多边形角点的编号方向是沿顺时针方向增加的。有

地球物理数据处理教程

可以看出

（S_oi）₂对于封闭的多边形KLMN……K，正好各项抵消为0，所以

地球物理数据处理教程

下面讨论S_i的计算，先将S₁化为极坐标形式如下：

地球物理数据处理教程

将多边形积分域分解成各条边与 P′（x_α，0，0）点所组成的三角形的代数和，设第i边KL所组成的三角形P′KL的积分值为S_1i，如图1.10，从上式可得

地球物理数据处理教程

图1.10 计算KL边所构成的三角形作用值的示意图

由图1.10中可以看出

地球物理数据处理教程

代入上式则有：

地球物理数据处理教程

式中各量计算如下：

地球物理数据处理教程

同样可以计算出其他各边的三角形域二重积分值，将它们取代数和即为多边形截面KLMN……K的二重积分值：

地球物理数据处理教程

类似的做法可以求出（1.2.4）式中的其他各二重积分S₂、S₃、S₄、S₅、S₆分别如下：

地球物理数据处理教程

地球物理数据处理教程

上列式中各量表示如下：

地球物理数据处理教程

注意：当多边形截面的任意一边平行于y轴或z轴时，以上公式中的分母就出现0的情况，从而使计算溢出，为此，我们将以上公式作一变形，将n_i、g_i、c_i、m_i直接代入，消去分母中只包含（y_i+1-y_i）、（z_i+1-z_i）的项，于是得到式（1.2.5）的另一等价形式：

地球物理数据处理教程

根据拉普拉斯方程，有

V₁+V₄+V₆=0

可得：

S₁+S₄+S₆=0

于是

S₁=-（S₄+S₆）

上面各式中：

地球物理数据处理教程

Y=y_i+1-y_i

Z=z_i+1-z_i

当二重积分求出以后，得到一组相互平行的截面上的S值，然后沿x方向作数值积分，如果我们用抛物线公式（即辛普生积分），那么有

地球物理数据处理教程

式中：j=0，1，2……6，积分限x₁～x_n为该地质体沿x轴方向上的两个端点的坐标值；S_j1、S_j2、S_j3分别为地质体在x=x₁、x₂、x₃截面上S_j的二重积分值；k为数值积分求和次数，即

地球物理数据处理教程

当n为偶数时，x₁～ x_n-1区间用辛普生积分公式，x_n-1～x_n区间用梯形积分公式。

计算得V₀、V₁、V₂……V₆以后，用式（1.1.3）和式（1.1.7）即可得到各种正演值。

1.2.4 “表面积分法”（多面体法）

该法是用一系列多边形小平面来逼近任意形体的全部表面形态，即用多面体代替任意形体，根据式（1.1.11）磁位是对全部表面进行积分而求得的，因此首先计算出每个多边形小平面的积分值，然后叠加求和，即为全表面的积分值。

首先计算边数为N_j的第j个多边形的积分值，为此建立新坐标系，使该多边形平面在新坐标系下能视为“水平面”，求得新坐标系中的三个轴向磁场分量，再通过坐标变换成原坐标系各磁场分量，再将所有的多边形平面的磁场值累加求和，即为所需正演值。

图1.11 新坐标系示意图

先讨论新坐标系的建立，在第j个多边形平面内，任选三个角点，分别为A_i-1、A_i、A_i+1，选定

为新坐标系的X′轴的方向，利用矢量积关系而得到另外两个坐标轴Y′和Z′的方向

地球物理数据处理教程

互相垂直，如图1.11所示，将它们分别作为新坐标系的X′、Y′、Z′轴的三个方向，过原点O分别以上述三个矢量为轴方向建立新坐标系为OX′Y′Z′，其中Z′轴垂直于第j个多边形平面，X′OY′平面平行于第j个多边形平面，因而这第j个多边形就在新坐标系中可视为“水平面”，新坐标系的三个轴在原坐标系中的方向余弦应为：

地球物理数据处理教程

式中：cosα₁、cosβ₁、cosγ₁为X′轴在原坐标系中的方向余弦；cosα₂、cosβ₂、cosγ₂为Y′轴在原坐标系中的方向余弦；cosα₃、cosβ₃、cosγ₃为Z′轴在原坐标系中的方向余弦；ξ_i-1、η_i-1、ζ_i-1及ξ_i、η_i、ζ_i分别为A_i-1和A_i角点在原坐标系中的坐标值；而B₁₁、B₁₂、B₁₃和C₁₁、C₁₂、C₁₃分别是下列B、C行列式中的i、j、k的代数余子式

地球物理数据处理教程

α₁、β₁、γ₁和α₂、β₂、γ₂及α₃、β₃、γ₃分别表示X′Y′Z′轴与X、Y、Z轴之间的夹角。

利用新旧坐标间的9个方向余弦，可将第j个多边形平面上的各角点和观测点的坐标转换成新坐标如下

地球物理数据处理教程

于是在新坐标系中，我们来求水平多边形磁荷面的磁场，首先将它分解成一系列与x′轴（或y′轴）相平行的梯形（即梯形的顶、底与轴平行），利用面磁荷积分公式，可计算磁场的各分量如下。

例如计算第j个多边形磁荷平面的x′轴方向水平分量

地球物理数据处理教程

式中：ζ′₁表示第j个平面的z′坐标值（常数）；S表示第j个多边形平面域，S₁、S₂……S_k表示所分解的共k个梯形域，将各梯形的积分限代入，不难得出

地球物理数据处理教程

式中：

同样y′轴和z′轴方向的分量为

地球物理数据处理教程

在计算出新的坐标系下的各分量以后，利用如下关系求得原坐标下的各分量

地球物理数据处理教程

求得第j个多边形磁荷平面的磁场各分量以后，按同样的计算方法计算出其他各个面的场值，然后累加求和得整个形体的场值。用一规格化公式表示如下

地球物理数据处理教程

式中：F代表H_αx、H_αy、Z_α中的任一个；其他符号除上述已解释的以外，还有

Ay_i=η′_i-y′

Ay_i+1=η′_i+1-y′

Az_i=ζ′_i-z′

Az_i+1=ζ′_i+1-z′

地球物理数据处理教程

n_j为第j个多边形平面的角点个数；m表示包围该形体的全表面的多边形平面总数（即多面体的面数）。

F为H_αx时：K₁=cosα₁，K₂=cosα₂，K₃=cosα₃

F为H_αy时：K₁=cosβ₁，K₂=cosβ₂，K₃=cosβ₃

F为Z_α时：K₁=cosγ₁，K₂=cosγ₂，K₃=cosγ₃

具体计算步骤归纳如下：

（1）建立一个z轴向下的原始坐标系，原点O可以任选，给出在该坐标系下，用来近似任意形体的多面体各面内所有角点的坐标值（ξ、η、ζ）。

（2）计算第j个面内的多边形平面：在该面内任选三个角点，以其中任意两个角点的连线作新坐标系中x′轴的方向，利用式（1.2.8）计算出新坐标轴的9个方向余弦。

（3）利用式（1.2.9），将第j个平面内的所有角点坐标转换成新坐标系的新坐标值（ξ′_i、η′_i、ζ′_i），并求出各边与x′轴的夹角θ_i，i+1和该面的σ值。

（4）用式（1.2.10）和（1.2.11）计算出H′_αx、H′_αy、Z′_α，即第j个磁荷面在观测点上沿新坐标轴的场值。

（5）用式（1.2.12）计算出第j个磁荷面沿原始坐标轴的各分量H_αx、H_αy、Z_α。

（6）改换平面为第j+1个，重复上述计算步骤（2）～（5），可逐一求得多面体的所有各面的场值，然后累加求和得到一个观测点的场值。

（7）改换观测点，重复上述步骤（2）～（6）即得整个测区内所有测点的正演值。