数值模拟方法和反演

第三篇数值模拟方法和反演引言已知场源分布和地下介质的物性分布，要求计算地表或地下场的分布，这就是地球物理勘探的正演问题，它是地球物理勘探资料解释的基础。

反过来根据所测到的地表或地下场的分布以及场的分布，来推断地下介质的特性分布，就是地球物理勘探中的反演问题，它就是对地球物理勘探资料的解释。

地球物理的正演问题，除了少数简单物性分布可用场论中的经典方法解析求解外，极大部分都需要用数值模拟方法求解。

这些方法包括有限单元法、有限差分法和边界元法。

有限差分法（Finite Difference Method—FDM）是从地球物理场所满足的偏微分方程和边界条件出发，将微分方程转变为差分方程。

其研究步骤是：首先将研究区域按一定方式离散化，例如对二维断面可以将研究范围划分为一系列小的长方形单元，然后在每个单元内假设物性均匀，地球物理位、场呈线性变化，因而微分方程中的微分就可用差分来代替，于是就可以建立一组线性差分方程，最后求解此线性方程组即得到相应的位、场分布。

有限单元法（Finite Element Method—FEM）也是从地球物理场所满足的偏微分方程和边界条件出发，根据微分方程的解与泛函极小问题的等价性，将微分方程和其边界条件转化为相应泛函的变分问题。

其研究步骤仍然是将研究的区域按一定方式划分离散化，例如对二维情况，可以将断面分割成一系列小的长方形单元或三角形单元，同样设单元内位、场呈线性变化，物性参数均匀，这是泛函是各节点位、场的二次函数，利用求极小的必要条件，即泛函对各节点位、场的变分为零，二次函数的变分为一次函数，由此可得到一个线性方程组，解此线性方程组便可求得各节点的位、场值。

有限元法由于可采用三角形网络划分研究区域，所以能比较容易弯曲的物性界面和起伏地形，这是它比有限差分法优越之处。

当然有限差分法也可通过较密的网格剖分来近似地模拟弯曲的物性界面和起伏地形，但这时要增加计算工作量。

边界单元法（Boundary Element Method-BEM）是在有限单元法以后发展起来的一种数值方法。

在边界单元法以前，就有所谓的边界积分方程法，它通过格林公式，将由偏微分方程表示的边值问题，转变成区域边界上的积分方程，然后，将场值集中在边界的若干点上（即用边界单元法中的所谓零次插值），求解积分方程。

有限单元法发展起来以后，其中的单元划分和插值函数的概念，被用来解边界积分方程法，发展成为边界单元法。

可以说，边界单元法是解边界积分方程的有限单元法。

地球物理反演问题求解方法可分为线性反演方法与非线性反演方法两类。

地球物理线性反演方法有最小二乘（最小方差）法、线性规划法、L 范数解以及广义线性反演法等。

非线性反演方法有:梯度法、尝试法、蒙特卡洛法、人工神经网络法、模拟退火法、遗传算法等。

第一章 地球物理模拟的有限差分法有限差分法是近似求解偏微分方程边值问题最常用的方法。

有限差分法解椭圆型偏微分方程的基本思想，是用离散的、只含有限个未知数的差分方程去近似代替连续变量的微分方程及边界条件，并把相应的差分方程的解作为边值问题的数值近似值。

有限差分法将连续问题离散化的步骤是首先对求解区域作网格剖分，用有限个网格节点代替连续区域，用数值微商代替微分，或以数值积分代替积分进行离散，从而把微分方程的定解问题化为线性代数方程组的求解问题。

本章将通过介绍二维均匀直流电场的有限差分法模拟来介绍有限差分法的具体作法。

一、差分方程的导出如图（3-1-1），均匀二维直流电场的边值问题可由下式表示BCnuAD n uAB E n uCD E n uΩz u z x u x ∈=∂∂∈=∂∂∈=∂∂∈-=∂∂∈=⎪⎭⎫⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂ 0 0000σσ （ 3-1-1）图4-1-1 均匀二维电场有限差分模拟区域及网格上式中假设y 轴平行地质体走向，E 0为已知的x 方向上均匀电场强度，故方程中没有y 变量。

u 为待求的电位。

有限差分法求解(3-1-1)式的作法是，首先将给定的二维地电断面离散化，通常采用平行于ｘ轴和平行于ｚ轴的直线划分断面成许多小长方形，每个小长方形称为单元，两平行线组的交点称为节点。

可将x 方向节点编号为i=1，2，…，N ；z 方向节点编号为j=1，2，…，M 。

将研究区域离散以后，由中心差商公式1121---∆-=⎪⎭⎫⎝⎛∂∂i i i i x u u x u ，i i i i x u u x u ∆-=⎪⎭⎫ ⎝⎛∂∂++121（3-1-2）则()⎥⎦⎤⎢⎣⎡∆-+∆-∆+∆=∆+∆⎪⎭⎫⎝⎛∂∂-⎪⎭⎫ ⎝⎛∂∂=⎪⎪⎭⎫⎝⎛∂∂+-----+i i i i i i i i i i i i ix u u x u u x x x x x u x u x u 11111212122221（3-1-3）式中各量的意义及其相应位置见图3-1-2。

图4-1-2 计算差分时的网格点位示意图相应于（3-1-1）式左端各项的差分形式为][2)(11111ii i i i i i i i i i x u u x u u x x x u x ∆-+∆-∆+∆=∂∂∂∂+----σσσ ][2)(11111jj j j j j j j j j j z u u z u u z z z u z ∆-+∆-∆+∆=∂∂∂∂+----σσσ 由此可得与（3-1-1）式对应的差分方程为0,1,1,,1,1=+++++-+-j i ijP j i ij B j i ij T j i ij R j i ij L u u u u u c c c c c（3-1-4）式中⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫+++-=∆∆+∆=∆∆+∆=∆∆+∆=∆∆+∆=--------)()(2,)(2)(2,)(21,111,1,11,1c c c c c c c c c ijB ij T ij R ij L ij P j j j j i ij B j j j j i ij T i i i ji ijR i i i ji ijL z z z z z z x x x x x x σσσσ （3-1-5）这里c ij L，c ij R，c ijT，cijB为相邻节点的连接系数，cij P为自身连接系数。

从上式可见，节点（i,j ）的电位u i,j 值只与(i-1,j )，（i+1,j ），(i,j-1)和(i,j+1)周围四个节点的u 值有关，其连接系数仅与网格形状大小和位置以及相应的介质导电率σ值有关。

二、边界条件的处理对于边界节点，考虑到各种边界条件，式（3-1-3）可进一步简化。

（1）对于地面节点即j =1行，满足诺依曼条件，即可设空中虚行j =0，j =0行的各节点为j =2行各节点对地面的镜象，设对应的u 值和σ值均相等，故有10z z ∆=∆ 2,0,i i σσ=2,0,i i u u =（3-1-6）（2）对地下边界节点j=M 行，同样满足诺依曼条件，此时可直接令1,,-=M i M i u u（3-1-7）（3）对左面和右面边界，由一阶中心差商（3-1-2）公式，可得1,1,1,2,1 --∆-=-∆=-N j N j N j j x E u u x E u u （3-1-8） 利用上面各式，可以建立每个节点的差分方程，于是可得到线性方程组，其矩阵形式为S CU = （3-1-9）式中C 为连接系数所组成的矩阵，U 为待求的各节点电位值组成的向量，S 为与源有关的右端项向量。

求解线性方程组（3-1-9），可得网格各节点电位值。

第二章 有限元法有限单元法是50年代从弹性力学中发展起来的方法。

主要优点是适用于物性参数复杂分布区域的数值模拟。

它的解题步骤如下：（1）给出地球物理边值问题；（2）将边值问题转化为变分问题；（3）用有限元方法解泛函极值问题；（4）解线性方程组。

下面将通过对于直流线源二维地电问题的有限元数值模拟来介绍有限元方法。

一、直流线源二维地电断面的边值问题设线源方向和地质体走向相同，并与y 轴并行，这样电位仅是x,z 的函数。

若线源位于地面P 点，则可得到电位满足下列微分方程Ω∈-=∇⋅∇ )(2)(p I u δσ（3-2-1）为简单计，假设计算区域Ω很大，从而边界处电位满足如下边界条件：Γ∈=∂∂ 0nu（3-2-2） （3-2-1）（3-2-2）就是直流线源的边值问题。

二、直流线源二维地电断面的变分问题方程（3-2-1）两边同乘以变分u δ后在区域Ω内积分：)(2)(Ω-=Ω∇⋅∇⎰⎰⎰ΩΩud p I ud u δδδσ（3-2-3）由于：)()()(u u u u u u δσδσδσ∇⋅∇+∇⋅∇=∇⋅∇用积分公式，并注意到边界条件（3-2-2）可得Ω∇⋅∇-=Ω∇⋅∇-Γ∂∂=Ω∇⋅∇-Ω∇⋅∇=Ω∇⋅∇⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰ΩΩΓΩΩΩd u u d u u d nuu d u u d u u ud u )( )( )( )()(δσδσσδδσδσδσ （3-2-4） 将（3-2-4）代入（3-2-3），并将变分号提出，可得：0])(2[2=Ω-∇⎰⎰Ωd u p I u δσδ也即边值问题（3-2-1）（3-2-2）与下面的变分问题等价⎪⎩⎪⎨⎧=Ω-∇=⎰⎰Ω0)(])(21[)(2u I d u p I u u I δδσ （3-2-5）三、有限元法第1步：区域剖分将积分区域剖分成许多三角单元，总节点数为N （图3-2-1），则（3-2-5）式中对区域Ω的积分可分界为对各三角单元△积分之和。

⎰∑∆ΩΩ-∂∂+∂∂=d u z x I zu x u u I p p })()(2])()[({21)(22δδσ (3-2-6)(a) (b)图3-2-1 区域Ω有限元网格剖分（a ）及三角单元（b ）示意图第2步：线性插值假设三角单元内电导率为σ，电位u 线性变化，即在每个三角单元△ ∆=++=U N T 332211u N u N u N u（3-2-7）这里，下标1、2和3是三角单元的三角节点号（参见图3-2-1b ）； T U ),,(321u u u =∆ 是三角单元节点上的电位；)(21,),(3,21l l l l c z b x a N N N N ++∆==TN 是x 和z 的线性函数),3,2,1(=l 其中⎩⎨⎧=个节点上时）位于三角单元另外二当（上时）位于三角单元节点当（x,z l x,z z x N l 01),( )(1221b a b a -=∆是三角单元的面积。

可以得出：122133113223321123312231213132321,,,,,,z x z x c z x z x c z x z x c x x b x x b x x b z z a z z a z z a -=-=-=-=-=-=-=-=-=式中))(,(2,211z x z x 和)(3,3z x 是三角单元三角节点的x 和z 坐标。