均匀任意形体正演方法

正演问题的模拟方法~

已知地电模型和场源分布，求解场的分布规律，称为电法勘探的正演问题。在学习电法勘探时，我们经常先对一些典型的地质模型进行正演模拟，从而建立地质模型与场分布特征之间的关系。因此，正演问题是学习电法勘探的重要基础。
解电阻率法正演问题有两个途径：一是通过物理模拟，即通过模型实验直接测量得到某种介质和场源条件下稳定电流场的分布情况；二是通过数学计算途径，即寻求满足表1.1.5边界条件下的拉普拉斯方程解。物理模拟方法主要有土槽、水槽、导电纸等手段。数值计算可分为解析法和数值计算方法两种。解析解有明确的表达式，但只有少数理想的地电模型才有解析解(这里指的理想模型是指数学意义上的，如均匀大地以及均匀大地中存在球体、两种介质垂直接触、水平均匀地层等，有些简单或规则模型不一定有解析解，如正方体等)，不能满足实际需要。目前，随着计算机和计算技术的发展，求解复杂条件下电磁场分布规律的数值计算方法已成为主要手段。主要数值计算方法有有限元法、有限差分法、边界元法和积分方程法，等等。
根据求解稳定电流场的唯一性定理，在求解稳定电流场分布的物理模拟和数值计算中，判断两个模型是否等价(即它们的电流场空间分布严格相似)，就要看这两个模型的稳定电流场是否可归结为相同的定解问题。显然，只要两个模型的尺寸按比例相似，各部分介质的电阻率也按比例变化，两个模型的稳定电流场就是相似的。换句话说，将野外模型中各长度因素按一定比例缩小，并按实际的电阻率比值，在实验室中建造物理模型并进行观测，这样所获得的电场和野外相同，这便是稳定电流场物理模拟的基本原理。

已知地电模型和场源分布，求解位场的分布规律，称为电法勘探的正演问题，它是电法反演问题的计算和观测资料解释的重要基础。在正演模拟问题中，地下电阻率分配被限定，其任务是计算地质体上方测线观测获得的视电阻率。正演实际上是反演的一部分，因为通过反演套路计算出的模型理论视电阻率值与实际观测值比较是否一致非常有必要，求解一个特定模型的视电阻率值方法主要有5种:(1)解析法；(2)有限差分法，(3)有限元法；(4)边界元法；(5)面积分法。解析法可能是最精确的方法，取得的结果具有典型意义，但是，解析法受多种条件限制，仅局限于相对简单的几何体电场分布问题才可以求解，如球体或圆柱体；边界元法比较灵活，但是不同电阻率值的区域数受到限制，一般少于10；近地表勘探时，地下可能存在任意分布的电阻率，因此，对于求解复杂条件下的电场分布规律，有限差分法和有限元法通常是可靠的选择方法，这些方法能将地下细分成数以万计的不同电阻率单元。当然，解析法和边界元法是非常有用的，可以单独用来检查有限差分法和有限元法的精度。
2.1.1 稳定电流场的边值问题
电阻率法的正演计算问题，均归结为求解稳定电流场的下述边值问题。
2.1.1.1 点源场中三维地电体的边值问题
在这种情况下，地下介质的电导率和电场均为三维空间的函数，即σ＝σ（x，y，z）和U＝U（x，y，z）。
（1）位函数U（x，y，z）所满足的微分方程
设在电导率为σ（x，y，z）的无限岩石中，位于A（xA，yA，zA）点处有一电流为Ⅰ的点电流源。由场论可知，在有源域内，可定义为

高密度电法勘探方法与技术

利用狄拉克δ函数的性质，则得

高密度电法勘探方法与技术

由j＝σE和E＝-▽U得稳定电流场位函数U所满足的微分方程为

高密度电法勘探方法与技术


高密度电法勘探方法与技术

当点电源A（xA，yA，zA）位于地表时，则上式应写成

高密度电法勘探方法与技术

对于地下岩石分区均匀情况，式（2.5）可写成

高密度电法勘探方法与技术

式（2.5）和式（2.6）分别为待求解的三维非均匀和分区均匀岩石中位函数U（x，y，z）所满足的微分方程。
（2）位函数U（x，y，z）所满足的边值和衔接条件
1）边值条件。由于电法所研究的稳定电流场分布于整个地下半空间，为了减少计算工作量，在求解边值问题时，通常把计算范围定在一个有限的区域内。这样，便需要在求解区域的边界Γ上，对电位函数U（x，y，z）赋予已知值。可见，为了求解位场问题，还必须知道所研究区域的边界上的位场分布状况，即边值条件。一般有三种类型的边值条件:
a. 第一类边值条件

高密度电法勘探方法与技术

式中:Г表示所研究区域Ω的边界；g（x，y，z）是定义在Γ上的已知函数。
b.第二类边界条件

高密度电法勘探方法与技术

式中:n是Γ的外法线。
c.第三类边值条件

高密度电法勘探方法与技术

式中:A为已知正数。
2）衔接条件。在所研究的区域Ω内，两种具有电导率为σ1和σ2的岩石交界面处，电位和电流密度法向分量满足以下衔接条件:
a. 由于电位的连续性，在交界面处有

高密度电法勘探方法与技术

b. 由于电流密度法向量的连续性，在交界面处有

高密度电法勘探方法与技术

式中:n为交界面的法线方向。
2.1.1.2 线场源中地电断面的边值问题
在此情况下，岩石的电导率和电场沿岩体走向（y方向）均无变化，只是坐标（x，z）的函数，σ＝σ（x，z）和U＝U（x，z）。此时，电位U（x，z）满足二维偏微分方程

高密度电法勘探方法与技术

式中:f＝-21σ（x-xA，z-zA）；Ⅰ为位于A（xA，yA）的线电流单位长度上的电流强度。
该问题中，边界上的位函数U（x，z）和已知函数g（x，z）同样满足式（2.7）～式（2.11）所示的边值和衔接条件。
2.1.1.3 点源场中二维地电断面的边值问题
在这种情况下，介质的电导率沿走向（y方向）无变化，仅是坐标（x，z）的函数，即σ＝σ（x，z），电场为三维空间的函数U（x，y，z）应满足方程

高密度电法勘探方法与技术

式中: 1= -2Iδ(Y-YA' Y - YA' Z-ZA) 。
对式（2.13）两端作傅里叶变换，得

高密度电法勘探方法与技术

式中: f1, = -1σ (x -xA， z -ZA) ; σ=σ( X， z);V=( X， λ ， z) λ 方空A波A。
这样，将三维的偏微分方程（2.13）变成了二维偏微分方程（2.14），即将（x，y，z）三维空间中的电位U（x，y，z）变换为（x，z）二维空间中的变换电位V（x，λ，z）。当分别对若干个给定的λ值求解方程（2.14），获得V（x，λ，z）后便可按下式作傅里叶逆变换，得到待求电位:

高密度电法勘探方法与技术

2.1.1.4 特殊条件下的边值问题
（1）无源域的边值问题
对于无源域，即在所求解的地电断面内无电源（点源或线源），方程（2.6），（2.12）和（2.14）中的位函数U（x，y，z），U（x，z）和V（x，λ，z）分别满足拉普拉斯方程和亥姆霍兹方程:

高密度电法勘探方法与技术

（2）分布均匀介质的边值问题
1）有源情况。由于所求解域内介质的电导率σ为常量，故有 ，所以式（2.6），式（2.12），式（2.14）可以分别写成

高密度电法勘探方法与技术

式中: 1= -2Ipδ(M-A); λ= -Ipδ(M -A) 。
2）无源情况。由于所求解域内介质的电导率σ为常量，且无电源，故式（2.16），式（2.17），式（2.18）可简写成

高密度电法勘探方法与技术

综上所述，计算传导类电法勘探的正演问题，或求解稳定电流场的边值问题，就是要根据情况确定未知的电位函数U（或其变换函数），使其在已知电导率σ分布的求解区域内满足相应的微分方程和边值条件。
2.1.2 解析法
在简单地电条件下求解位场分布时，常用解析法中的分离变量法和类比法。
（1）分离变量法
分离变量法是求解电磁场边值问题常用的解析方法，该法的具体内容是假定待求的位函数由2个或3个仅含1个坐标变量的函数乘积所组成，然后将该乘积关系代入所求解的偏微分方程，经过整理，通过分离常数的联系，原来的偏微分方程可转化为几个常微分方程，其数目等于独立坐标变量数。解这些常微分方程并以给定的边值（包括边界上的衔接条件）决定其中的待定常数和函数后，即可求得待定的位函数。
（2）类比法
在边值问题的分析计算中，根据位场解答的唯一性定理，各种物理场，不论其对应物理量的意义是否相同，只要它们具有相似的微分方程和边值，则所得的解在形式上必完全相似，因而就把某一位场的分析计算结果，推广到一切相似的位场中去。
另外，在特定条件下，求解电法勘探的边值问题常用的解析方法还有镜像法和保角变换法等。
2.1.3 数值计算方法
解析法只能求解少数规则形体的边值问题，对于复杂条件下电场分布，常需借助数值计算方法。目前，用于电法勘探的数值计算方法有:有限差分法、有限单元法、面积积分法和边界单元法等，下面将对这些数值计算方法的原理进行简要的介绍。
2.1.3.1 有限差分法
有限差分法是以差分原理为基础的一种数值计算方法，它用各离散点上函数的差商来近似替代该点的偏导数，把要解的边值问题转化为一组相应的差分方程，然后，解出差分方程组（线性代数方程组）在各离散点上的函数值，便得到边值问题的数值解。现以二维等步长差分格式为例，说明有限差分法的原理和方法步骤。
（1）区域离散和网格剖分
如图2.1所示，用平行于坐标轴的两组直线族将地下划分成正方形网格，相邻两坐标线的距离为h，则任一节点的（x，z）坐标为

高密度电法勘探方法与技术


图2.1 二维等步长正方形网格

每个正方形为一单元，其边长为h，称为步长，网格的交点称为节点。任一节点的坐标（x，z）可表示为（ih，kh），或简化为（i，k），用阶梯状折线代替原来的曲线段。在边界线以内的节点称为内节点，边界上的节点称为边界节点。
（2）微分方程离散化和构组差分方程某一内节点（i，k）处的电位为U（i，k），由于h很小，可将节点（i，k）四周的电位在节点处展成泰勒级数:

高密度电法勘探方法与技术

式中:Ux，Uxx…和Uz，Uzz…分别表示U对x和z的一阶导数、二阶导数等。
将前两个式子相加，并且忽略h的四次项和更高次项，经整理可得

高密度电法勘探方法与技术

同理可得

高密度电法勘探方法与技术

将上述Uxx和Uzz代入含源分区均匀岩石中位函数U所满足的微分方程（2.19）的第二式，即得二维函数U（x，z）的差分方程

高密度电法勘探方法与技术

对于无源分区的均匀介质，位函数U（x，z）所满足的微分方程（2.20）的差分方程为

高密度电法勘探方法与技术

（3）线性方程组的形成与求解
对于边界节点，其相应的差分方程可根据边界条件给出。全部节点所建立差分方程（2.26）和（2.27）的总和可分别写成以下矩阵形式:

高密度电法勘探方法与技术

和

高密度电法勘探方法与技术

式中:［A］是方程组的系数矩阵，它是与电阻率分布有关的函数；｛U｝ 是电位U的列向量，其分量为所有节点上的电位；｛F｝ 是常向量。
当给定电阻率分布及边界条件后，解线性方程（2.28）和（2.29）便可求得电位的空间分布。
电位 ｛U｝ 值的计算精度与步长h的大小有很大关系。一般来说，网格划分越细，即h值越小，｛U｝ 值与理论值就越接近。但是，此时节点数目也急剧增加，因而，所需的计算机内存和计算时间也就会增大。解决计算机速度与精度这一矛盾的较好方法是采用变步长，即在近区将网格划分得密些，远区影响较小可划分得稀些。
图2.2为一个采用不同装置的2-D高密度电法勘探视电阻率拟断面有限差分正演模拟例子，其中图2.2a为地质模型，图2.2b为温纳（Wenner）装置获得的视电阻率拟断面，图2.2c为单极-单极（pole-pole）装置获得的视电阻率拟断面，图2.2d为偶极-偶极（dipole-dipole）装置获得的视电阻率拟断面。
图2.3为一3-D高密度电法勘探视电阻率拟断面有限差分正演模拟例子，其中图2.3a为5×15布设的电极测网（即225根电极）和3-D地质模型，在电阻率为50Ω·m的模型介质中嵌入了4个四方形棱柱，图为穿过地质体的水平切片。图2.3b为单极-单极（pole-pole）装置的视电阻率值水平方向的切片，注意，在网格中心，从1.0～3.2 m深处有一电阻率为10Ω·m的低阻体，该低阻异常体是由于极距小于4 m造成的；但是，对于电极距大于6m，低阻棱柱成了高阻异常，这是一个由于C1和P1电极间，在近地表区域为负灵敏度出现了 “异常翻转” 的现象。
2.1.3.2 有限单元法
有限单元法是以变分原理和剖分插值为基础的数值计算方法。用这种方法解稳定电流场问题时，首先利用变分原理把所要求解的边值问题转化为相应的变分问题，也就是所谓泛函的极值问题。然后，与有限差分法相似，使连续的求解区域离散化，即按一定的规则将求解区域剖分为一些在节点处相互连接的网格单元，进而在各单元上近似地将变分方程离散化，导出以各节点电位值为未知量的高阶线性方程组，最后求解此方程组，算出各节点的电位值，以表征稳定电流场的空间分布。
（1）稳定电流场的变分问题
由前面讨论可知，稳定电流场的计算可归结为求解电位所满足的偏微分方程边值问题式（2.5）至式（2.20）。为了用有限单元法求上述边值问题的数值解，首先将其转化为相应的变分问题，即构造相应的泛函表达式。
1）点源三维情况。根据前述稳定电流场的边值问题，点源三维稳定电流场的计算，可归结为求解下列的边值问题:

图2.2 一个采用不同装置的2-D高密度电法勘探视电阻率拟断面有限差分正演模拟例子


图2.3 一个3-D高密度电法勘探正演模拟例子（测网为15×15）


高密度电法勘探方法与技术

式中:Γ2为求解区Ω的地面边界；Γ1为求解区Ω的其余边界。
与式（2.30）问题等价的变分问题为

高密度电法勘探方法与技术

顺便指出，边值问题式（2.30）中的边值条件 在变分问题式（2.31）中没有出现，因为根据变分原理，这一边值条件将自然得到满足。
2）线源二维情况。与点源三维情况相似，二维地电条件下，稳定电流场的边值问题为

高密度电法勘探方法与技术

与其等价的变分问题为

高密度电法勘探方法与技术

式中:S为二维（平面上）的求解区；Г＝Г1 ＋Г2构成求解区S的全部边界线。
与三维情况相似，应用变分，可以证明式（2.32）与式（2.33）的等价性；并且，同样可以看到，式（2.32）中的边值条件 在变分问题式（2.33）中也是自然被满足。
3）点源二维情况。对于点源二维地电条件，所求解的稳定电流场的边值问题为

高密度电法勘探方法与技术

式中
与二维偏微分方程边值问题式（2.14）等价的变分问题为

高密度电法勘探方法与技术

求出变换电位V（x，y，z）之后，便可按下式作傅里叶逆变换计算电位:

高密度电法勘探方法与技术

（2）区域离散和网格剖分
有限单元法与有限差分法类似，需将连续求解离散化和作网格剖分。在二维情况下，最常用的方法是将求解区剖分为一系列互不重叠的三角形单元，每个单元的顶点称为节点。如图2.4所示，剖分求解区域的规则是:单元的节点只能与相邻单元的节点相重，而不能成为相邻单元的边内点；每个单元内介质的电导率为常数，在不同电导率介质的分界面（内边界）上或求解区的（外）边界上，使三角形单元的一个边互相衔接，以这样形成的折线逼近边界线；三角形单元的三条边长（或三个顶角）之间差异不要过大；在电场变化剧烈、电位参数的二阶或高阶导数值大的地方，应使单元剖分得细一些，而在电场变化平缓的地方，可剖分得稀些。
（3）线性插值
为了计算二维变分问题式（2.33）中积分形式的泛函J（U），需要知道求解域内的U函数值，通常是利用各节点函数值在各单元内作线性内插，按单元逐个计算U值。
如图2.5所示，设第e个单元的3个节点按逆时针方向编号分别为i，j，m，其坐标分别记为（xi，zi），（xj，zj），（xm，zm），对应的节点函数值为Ui，Uj，Um。假设各单元内，函数U（x，z）是线性变化的，即

图2.4 有限元网格剖分示例


图2.5 三角形单元


高密度电法勘探方法与技术

式中:系数a，b，c可由单元上3个节点的函数值和坐标算出。对于单元e，3个节点上的U都应满足式（2.37），故有

高密度电法勘探方法与技术

按克莱姆法则求解上述方程组，得

高密度电法勘探方法与技术


高密度电法勘探方法与技术

将式（2.39）代入式（2.37），便得到单元e内函数U的线性插值表示式:

高密度电法勘探方法与技术

式中:

高密度电法勘探方法与技术

如果在每个单元上都按此给出U（x，z）的近似表达式，便可得到整个求解区内U（x，z）的总体近似函数，这个函数在各个单元内是线性的，即它是一个分片线性函数。对于任意两个相邻单元来说，近似函数在公共边上的值被两个端点的节点函数值唯一确定，所以，总体近似函数在整个求解区自然是连续的。
（4）变分问题的离散化
二维变分问题式（2.33）中的泛函——整个求解区上的积分J（U），可以分解为各个单元上的积分Je（U）之和:

高密度电法勘探方法与技术

式中:Je（U）不难根据各单元e内函数U的线性插值近似式（2.40）算出。
1）单元分析。根据式（2.33），并考虑到单元内σ为常量，则可将单元e上的泛函写成:

高密度电法勘探方法与技术

式中:σe为单元e内的电导率值；ⅠA为电流源强度。
由式（2.40）可得

高密度电法勘探方法与技术

因而式（2.42）可用矩阵表示为

高密度电法勘探方法与技术

式中:

高密度电法勘探方法与技术


高密度电法勘探方法与技术
;n 方以i 方A庶的A元小敏。
由式（2.44）可见，［Ke］为单元系数矩阵，它是一个对称矩阵，其元素的一般表达式为

高密度电法勘探方法与技术

2）总体合成。按式（2.41）对所有单元的Je（U）相加，并将式（2.44）代入，便可得到整个求解区的泛函J（U）关于节点电位的离散表达式:

高密度电法勘探方法与技术

这样，泛函J（U）就完全离散化为多元二次函数:

高密度电法勘探方法与技术

所以，将变分问题式（2.33）化为二次函数的极值问题，即

高密度电法勘探方法与技术

根据函数极值理论，应有

高密度电法勘探方法与技术

则由式（2.48），即得

高密度电法勘探方法与技术

或表示为矩阵形式

高密度电法勘探方法与技术

式中: ，其元素为所有单元矩阵［Ke］的相应元素之和；列向量 ，其元素也为所有单元场源列矢量 ｛Ie｝ 的相应元素之和；｛U｝＝｛U1，U2，…，UN｝T。
式（2.52）经边界条件处理，应修改为

高密度电法勘探方法与技术

方程组的系数（即刚度矩阵［K'］的元素）决定于节点坐标及电导率的分布，是已知的；右端项也是给定的，于是求方程组（2.53），便可获得待求电位值。
图2.6为一2-D高密度电法勘探视电阻率拟断面限单元法正演模拟例子，采用温纳（Wenner）β装置，图2.6a为视电阻率拟断面，图2.6b为地质模型。
2.1.3.3 边界单元法
边界单元法是以边值问题控制微分方程的基本解为基础，首先建立边界积分方程，然后，在区域的边界上划分单元，进行数值求解。
（1）以加权剩余格式建立边界积分方程
加权剩余法已成为求解数学物理方程的普遍方法，为了一般化，现以加权剩余格式，由控制方程及边值条件出发，建立边界积分方程。使用以下记号

高密度电法勘探方法与技术

上两式分别表示函数u和v的乘积在域Ω和边界г的积分，称（u，v）Ω和（u，v）г为u和v的内积。
对于自伴随线性微分算子，通过分部积分，总可以得到如下形式的内积公式:

高密度电法勘探方法与技术


图2.6 一有限单元法高密度电法勘探正演模拟例子

式中:S为基本边界条件算子，G为自然边界条件算子。该算子相应的微分方程边值问题可列为

高密度电法勘探方法与技术

式中:g，s是对应于相应边界算子的给定函数；Г1，Г2则代表给定s或g的边界。
当用加权余量法求解时，对微分方程及对应边界条件，选取近似解u代替精确解u0，从而得到余量（或误差）函数:

高密度电法勘探方法与技术

适当选取权函数w，利用加权，使余量接近于零，从而总体意义上满足域内及边界要求，即设

高密度电法勘探方法与技术

如果在域Ω及其各边界段分别设权函数w为

高密度电法勘探方法与技术

利用式（2.56）和式（2.64），可将式（2.63）写成

高密度电法勘探方法与技术

该式更紧凑的形式为

高密度电法勘探方法与技术

在边界元法中，采用算子的基本解作为权函数w*，即采用满足如下方程的奇异函数:

高密度电法勘探方法与技术

式中:δi是狄克函数，它在所考虑的“i” 点邻域上的积分等于1，而在其他地方的积分为零，且有

高密度电法勘探方法与技术

将式（2.67），（2.68）代入方程（2.66）中，即得

高密度电法勘探方法与技术

式中:ci是与Ω的边界所确定的张角成正比的系数。
若记Bi＝〈f，w*〉，则式（2.69）可写成

高密度电法勘探方法与技术

式（2.70）即为待求的边界积分方程。
1）点源场三维地电体位场问题。假设电阻率分别为ρ1和ρ2的围岩介质和三维地电体分布在域Ω1和Ω2内，其中位函数分别表示为U0（x，y，z）和V0（x，y，z）。具有电流为Ⅰ的点电源位于光滑地表A处。位函数精确解U0（x，y，z）和V0（x，y，z）满足边值问题:
控制微分方程:
自然边界条件:
本质边界条件:
当用位函数的近似解u和v分别替代式（2.71）至式（2.74）和相应的衔接条件式（2.10），式（2.11）中的精确解U0和V0，并取权函数w*＝u*，S（w*）＝u*，G（w*）＝ 时，代入方程（2.70）后，可得该边值问题的边界积分方程:

高密度电法勘探方法与技术

式中:u*＝1/4πr为拉普拉斯方程的基本解； 。
2）点源场二维地电体的位场问题。假设具有电流为I的点电流源A位于光滑地表，分布在域Ω1和Ω2的围岩介质和二维地电体的电阻率分别为ρ1和ρ2，其内的电位函数分布为U0（x，y，z）和V0（x，y，z），如图2.7所示。三维位函数U0和V0满足如下的边值问题:
控制微分方程:
自然边界条件:
本质边界条件:
式中: 分别为Г1和Г2边界上的已知函数。

图2.7 点源场二维地电断面

对上述各式进行余弦傅里叶变换，可将三维的位函数U0（x，y，z）和V0（x，y，z）转变为二维的位函数u0（x，y，z）和v0（x，y，z），并满足边值问题:
控制微分方程:
自然边界条件:
本质边界条件:
式中:λ为傅里叶变换量（或波数）。
当用位函数的近似解u和v分别代替式（2.81）～式（2.84）和衔接条件式（2.10），（2.11）中的精确解u0和v0，并取与点源场三维地电体位场问题中相同形式的权数w*，S（w*），G（w*）时，代入方程（2.70）后，可得该边值问题的边界积分方程:
式中 是亥姆霍茨方程和基本解。
3）均匀场二维地电体位场问题。给定的地电断面如图2.8所示，地下围岩介质的电阻率ρ1，矿体的电阻率ρ2，并分别属于域Ω1和Ω2，其电位精确解分别用u0和v0表示。根据式（2.20），位函数u0（x，z）和v0（x，z）应满足如下边值问题:
控制微分方程:
自然边界条件:
本质边界条件:
式中: 分别为Г1和Г2上的已知边值函数。
当用位函数的近似解u和v分别替代式（2.87）～式（2.90）中的精确解u0和v0，并取相应的权函数，代入方程（2.70）后，可得该边值问题的积分表达式:

高密度电法勘探方法与技术

式中: 是二维拉普拉斯方程的基本解。

图2.8 均匀、稳定电流场二维地电断面

（2）确定离散方程与基本方程
为了利用边界积分方程（2.70）计算待求电位ui值，需要将其离散化为线性代数方程组。为此，将边界Г划分成许多小单元，称之为边界元。假设在各单元上，取u和q值为常数，且用中结点值表示，则称常数元（图2.9a）；若u和q值在单元上呈线性变化，则称线性元（图2.9b）；如用二次曲线来逼近边界，则称为二次元（图2.9c）。下面以线性元解法为例来说明离散方程与基本方程的确定方法。

图2.9 边界元类型

因为线性元的u和q在单元上呈线性变化，故在两个直线元的交点（端结点）上赋已知值（图2.10）。假定在Г1的端结点上，q是已知的，在Γ2的端结点上u是已知的，根据积分的可加性，对于结点“i”，可将式（2.70）中的积分化为每个单元Гj的积分和，即得离散方程:

图2.10 线性元


高密度电法勘探方法与技术

对于平面域的边界线元Гj，因为u和q在单元上是线性变化的，所以u和q在单元上任一点的值，可借助与端结点的值和两个线性插值基函数来确定:

高密度电法勘探方法与技术

式中:ξ是无量纲坐标， 。
将式（2.94）和（2.95）式代入（2.93）后，其中的单元积分可写成:

高密度电法勘探方法与技术

式中:

高密度电法勘探方法与技术

要写出对应于节点“i”的离散方程形式，需要把两个相邻单元Гj和Гj-1的贡献加起来，合成一项，以决定j节点上uj和qj的系数:

高密度电法勘探方法与技术

于是，式（2.93）对应于“i” 的装配方程，可写成:

高密度电法勘探方法与技术

对于空间域Ω的边界面元Гj，假定三角形单元Гj的3个顶点坐标为（x1，y1，z1），（x2，y2，z2）和（x3，y3，z3），则三角形内任一点的坐标（x，y，z）为

高密度电法勘探方法与技术

由于假定u和q值在每个单元呈线性变化，所以，可借助与三角形端结点的值和线性插值基函数，确定该单元内任一点的u和q:

高密度电法勘探方法与技术

式中:φk＝ξk。于是，方程（2.93）中的单元积分，可写成

高密度电法勘探方法与技术

式中:
利用式（2.103），式（2.104），对应于节点i，式（2.93）可写成:

高密度电法勘探方法与技术

对于方程（2.99）和方程（2.105），若记 ，则上述方程可分别写成如下更紧凑的形式:

高密度电法勘探方法与技术

同常数元解法一样，若把全部“i” 点合起来，方程（2.106）和方程（2.107）分别构成N阶和M阶线性代数方程组，可表示为矩阵方程

高密度电法勘探方法与技术

对于线积分hij和gij运用四点或七点高斯求积，用高斯消元法求解线性方程组（2.108），便得待求的二维电位。对于点源二维问题，尚须使二维电位U（x，λ，z）通过傅里叶逆变换，以获得要求的三维电位值U（x，y，z）。最后，根据所采用的电阻率法装置类型和大小，按公式ρs＝K△UP1P2/I计算视电阻率异常值。
2.1.3.4 面积分法

图2.11 用面积分法求解位场示意图

面积分法是从积累电荷的概念出发，通过求解积分方程，以确定电场空间分布的一种数值计算方法。如图2.11所示，在电阻率为ρ1的半无限介质中，埋有电阻率为ρ2的三维地质体，由场论中静电场与稳定电流场的类比关系可知，求解稳定电流场，可利用电阻率分界面上积累电荷的面分布。为此，设地下某一点的电位U（M）为

高密度电法勘探方法与技术

式中:U0（M）为供电电极A在M点产生的正常电位。即

高密度电法勘探方法与技术

Ua（M）为地质体表面的积累电荷对M点产生的异常电位:

高密度电法勘探方法与技术

式中:σ（Q）为地下矿体表面S上任一点Q的积累电荷面密度。
为镜像矿体表面 积累电荷对M点产生的异常电位，有

高密度电法勘探方法与技术

因而，地下任一点M的电位表达式（2.109）可写为

高密度电法勘探方法与技术

可见，根据边界条件求出积累电荷的面分布，便可有上式计算出电位值。

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

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

（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）即得整个测区内所有测点的正演值。

---

毕节市18745828353： 什么是正演模型?定义和使用方法 - ？
汲该已烯： 在地球物理磁法勘探的理论研究中,根据磁性体的形状、产状和磁性数据,通过理论计算、模拟计算或模型实验等方法,得到磁异常的理论数值或理论曲线,统称为正演问题.在异常推断和钻探验证的基础上,常常根据已掌握的磁性体的资料进行正演计算,用来修正推断的结果,使其更符合客观实际. 希望对你有用!

毕节市18745828353： 谁能区分一下“正演”和“反演”啊? - ？
汲该已烯： 在不同学科正反演定义有微小差别,但大体上一致的. 举个几个的例子: 在地球物理中,已经知道地球介质的性质,如(地震波传播速度等),求地震波的走时(即地震波在地球中的传播时间等)这是正演;又如,已知地下介质结构,物质特...

毕节市18745828353： 如何用MATLAB做出地震正演模型 - ？
汲该已烯： 1.在matlab的命令窗口里输入Simulink,回车,打开Simulink的功能模块函数库窗口. 2.在功能模块函数库窗口,选择菜单File-New-Model,就可以打开一个空白设计区域. 3.在功能模块函数库窗口中双击“Sources”打开一个子函数库,用鼠标把“Sine wave...

毕节市18745828353： 求AX=B的正演matlab程序 - ？
汲该已烯： >>A=zeros(10) %产生一个10*10的全零矩阵>>for i=1:10; for j=1:10; A(i,j)=1/(i+j); end;end;>>A>>X=zeros(1,10); for i=1:10; X(i)=(-1)^i; end;>>X; B=A*(X') 这个是正演的程序.反演???有点麻烦,还没想好.

毕节市18745828353： 京剧别誉为国粹,它的基本表演手段有哪些?行当有哪些? - ？
汲该已烯： 表演手段: 四功,唱念做打 五法,手眼身发步 京剧行当的划分,除依据人物的自然属性(性别、年龄)和社会属性(身分、职业)外,主要是按人物的性格特征来分类.京剧班社旧有“七行七科”之说 : 七科:为音乐、剧装、容妆、盔箱、剧...

毕节市18745828353： 元杂剧的艺术表现形式是什么? - ？
汲该已烯： “四折一楔子”是元杂剧最常见的剧本结构形式,合为一本,每个剧本一般由四折戏组成,有时再加一个楔子,演述一个完整的故事.少数作品也有一本分为五折或六折的,还有用两个楔子的.通常一本就是一部戏,个别情节过长的戏,可写成...

毕节市18745828353： jason里怎么做正演 - ？
汲该已烯： 跟反演的流程是一样的,先要建立模型,在输出的时候选择合成记录选项!

毕节市18745828353： 地球物理正反演的含义是什么?作用如何? - ？
汲该已烯： 通过这些受局限的观察信息推演相关过程发生的原因或机制,就称为“反演”.由于事件或过程发生在先,而结果或信息接收在后,对自然事件或过程发生的描述和预测被称为“正演”...

毕节市18745828353： 国粹京剧是大众( )的曲艺形式,人物分为()四种.锣鼓敲响,他们一个个(),在台上(),唱的(),演的 - ？
汲该已烯： 国粹京剧是大众(喜闻乐见)的曲艺形式,人物分为(生旦净丑)四种锣鼓敲响,他们一个(粉墨登场),在台上(唱念做打),唱的(字正腔圆),演的(惟妙惟肖).

毕节市18745828353： 京剧花脸分为几类 - ？
汲该已烯： 一、主要分为五类,有关的解释如下: 花脸里又分大花脸、铜锤花脸、架子花脸、武花脸、摔花花脸……京剧花脸脸谱大花脸是正净的俗称,过去讲大花脸多以铜锤花脸指代,唱功较多,重头戏较多,过去的架子花脸的地位偏低,由于几代花脸...