直接解拉氏方程位场转换[5]
发布网友
发布时间：2022-10-01 04:11
共1个回答
热心网友
时间：2024-11-05 09:13
这是采用以下的边界条件
地球物理数据处理教程
来求解拉氏方程。并取地形曲面z=f（x，y）为上、下半空间的分界面，代替传统的取无穷大平面为上、下半空间的分界面。
于是，得到这组新边界条件下的解析解
地球物理数据处理教程
这个解比较简单，因而在实际资料计算时比较方便，也有较高的精度。
利用这一组解析解，可以进行各种位场转换。
3.2.1 以曲面为边界的拉氏方程解
重磁场的位和各分量在没有场源的空间，是满足拉普拉斯方程的，一切问题总是可以归纳为按具体的边界条件求拉氏方程的解。以往考虑该问题的范围是整个上半空间，但实际测量大多是在地面上进行的，而地面常是起伏的曲面，而且异常值仅存在于局部空间，因此考虑整个上半空间没有必要，而且往往会引起具体异常边部的各种计算减低精度，所以改用以下的边界条件更为合适。
地球物理数据处理教程
式中：u为位或场的分量值；Lx为异常在x方向的长度；Ly为异常在y方向的长度。
在上述4个边界范围内的空间，以实际地形面划分为上、下两半部分，设实际地形面的方程为z=fS（x，y），且
地球物理数据处理教程
由于实测场值为离散值，故有
u（xR，yR，zR）=φ（xR，yR）
式中：R=1，2，3，…，ω，ω 为总测点数。在分界面以上，仍取 z→∞为最后一个边界，有
u（x，y，z）z→∞=0 （3.2.4）
下面讨论拉普拉斯方程的解
地球物理数据处理教程
通过分离变量法，有
u（x，y，z）=R（x）·S（y）·T（z）
代入式（3.2.5），有
地球物理数据处理教程
从而得到三个常微分方程
地球物理数据处理教程
而σ2=λ2+δ2
可由所取的边界条件求这三个常微分方程的解，方程（3.2.7）在边界条件（3.2.1）下的解为
R（x）=sin（λx）
其中
地球物理数据处理教程
对于i的确定值以λi代入，得R（x）的通解
地球物理数据处理教程
式中：αi为任意系数；m为任意大但非无穷的正整数。
同理，方程（3.2.8）在边界条件（3.2.2）时的通解为
地球物理数据处理教程
式中：
地球物理数据处理教程
βj为任意系数，n为任意大但非无穷的正整数。
方程（3.2.9）在边界条件（3.2.4）时的解为
地球物理数据处理教程
而
地球物理数据处理教程
于是得到
地球物理数据处理教程
即
地球物理数据处理教程
从上式可见，其解呈级数形式，并呈指数衰减。其衰减速度是很快的，在实际计算时m、n均无需取得很大。
按以往的做法，取上半空间的边界条件，则其解为
地球物理数据处理教程
这个解比式（3.2.13）要复杂，因它衰减慢，故项数必须取得较多，加大了计算工作量，而且在异常的边部，出现较大的误差。新边界条件（3.2.1）式和（3.2.2）式的选用，得到了一个更为简便的解（3.2.13），使计算工作量减少了而精度又有较大的提高。
3.2.1.1 当测量在一个平面内进行时
有：
地球物理数据处理教程
式中
地球物理数据处理教程
系数A′ij可用简单积分方法求得，将式（3.2.14）两边同乘以sin（λkx）sin（δLy），并对x和y在延拓得的新区域积分，得
地球物理数据处理教程
令
地球物理数据处理教程
故
地球物理数据处理教程
当x=0和Lx、y=0和Ly时，θ=0和π，ψ=0和π
代入上式的右边，得
地球物理数据处理教程
于是
地球物理数据处理教程
用i代替K，j代替L，并代回Aij得
地球物理数据处理教程
通常只能得到φ（x，y）的离散值φ（xR，yR），所以积分
地球物理数据处理教程
只能用数值积分予以完成。
（1）当测量按矩形网格进行时，若测线数为m1，测点数为n1，点距为D，点线距之比为μ，则有：
地球物理数据处理教程
得
地球物理数据处理教程
其中φ（xu，yv）为第u线第v点上的实测值。
而
地球物理数据处理教程
（2）当测量不按矩形网格进行时，可将各测点按最近距离联成三角形网格，若连接的三角形数目为S0，则显然有
地球物理数据处理教程
其中，Δu为第u个三角形。
可证明，右边的积分近似等于同一积分域（三角形Δu内）的这样一个积分：其被积函数H（u）为常量，且等于φ（x，y）sin（λix）sin（δjy）在三角形三个顶点（α，β，v）上的平均值，即
地球物理数据处理教程
于是得
地球物理数据处理教程
而
地球物理数据处理教程
式（3.2.15）和（3.2.16）有相同的计算精度。
3.2.1.2 测点分布于曲面上
即φ（x，y）沿任意起伏地形曲面测量时，由式（3.2.3）和式（3.2.13）得下列方程组
地球物理数据处理教程
取m、n使m×n=W，由于通常W值很大，从所列方程组解出m×n个系数是困难的，所以用最小二乘法建立法方程的方法求解，为此令 F（A11，A12，…，Aij，…，Am（n-1），Amn）为上列W个方程右边项移往左边后的平方和，即
地球物理数据处理教程
根据极小的条件为
地球物理数据处理教程
由此得出m×n个包含Aij的线性方程组
地球物理数据处理教程
这是m×n阶的方程组，用来解出m×n个系数Aij。
3.2.2 二维情况
对于二维的情况，相应的拉氏方程及边界条件应为
地球物理数据处理教程
同样采用分离变量法，考虑到φ（x）是以离散形式出现的，类似于三维情况，可得
地球物理数据处理教程
3.2.2.1 当z=fL（x）为直线时
即z=z0，φ（x）是在一直线上测得的，有
地球物理数据处理教程
式中：A′i=Ai
对于系数 A′i，可用简单积分的方法求得。用 sin（λkx）（k=1，2，…，m），乘以（3.2.19）式的两边，并对x在延拓后的新区间积分，得
地球物理数据处理教程
令θ=
则dx=
当x=0和Lx时，有θ=0和π，代入右边的积分，得
地球物理数据处理教程
于是得
地球物理数据处理教程
以i代替K，并代回Ai得
地球物理数据处理教程
一般来说，φ（x）呈离散形式φ（xj），则积分
地球物理数据处理教程
只能用数值积分形式完成。
（1）当有相同的测点距D时
地球物理数据处理教程
得
地球物理数据处理教程
式中
地球物理数据处理教程
（2）当测点距不相同时
地球物理数据处理教程
式中Dj=xj+1-xj，为第j段的长度。得
地球物理数据处理教程
3.2.2.2 当φ（xj）是在曲线上测得时
即沿地形起伏曲线测量而得，则有下列方程
地球物理数据处理教程
若取m=n，则可唯一解出Ai，但这样既费时间又无必要，一般取m≪n，用最小二乘法建立法方程求解，令F（A1，A2，…，Am-1，Am）为上述n个方程右端项移往左边后的平方和，即
地球物理数据处理教程
按极小必要条件得
地球物理数据处理教程
得m个含Ai的线性方程组
地球物理数据处理教程
地球物理数据处理教程
得到m个方程的方程组，求解m个系数Ai。
由上可见，从曲线和曲面为上、下半平面和半空间的分界面求拉氏方程的解出发，可以很方便地计算各种位场转换量。因为拉氏方程的解，其每一项均是一个指数函数和两个正弦函数的乘积，其每个因子分别是一个变量的函数，因而是高次可微可积的，所有位场转换的计算无非是一些微分和积分，因此用这种方法可以很容易建立起一个完整的位场转换系统，它包含现有的一切位场转换，还可以进行新的位场量的计算。
3.2.3 有关技术方法
3.2.3.1 影响计算精度的主要因素
（1）在边界条件（3.2.1）和（3.2.2）中，x=0，x=Lx，y=0，y=Ly定义了一组直立边界面，而对于某一高度z=z0而言，边界面上的位场值为零的假设就不成立了，因为异常为零的边界，随着高度的增加，要向外扩展，应为x=αz+α1，x=βz+α2，y=αz+b1，y=βz+b2的斜面为位或场的零边界，如图3.10所示。这样一来，使拉氏方程的解十分复杂和困难，即使求得了一些数字解，也将会有较大的误差，而且不同场源的异常，其零边界随高度向外扩展的规律也不相同，α、β不能准确地确定，只能近似地取值。而取直立边界使问题变得十分简单，当取Lx、Ly比实际边界稍大时，近似的程度也是可行的。
图3.10 零异常边界的图示
（2）测区场值零边界的确定，在逻辑上显然需要首先确定区域背景场，以便从实测值中除去，然而这确是难于严格做到的。实际上，区域背景场定量可靠的求取方法还是一个待解决的问题。取测区之外某一边界的场值为零，看来好像避开了这个问题，其实无非是对所述前提作了一个近似的假设而已。当然，零边界取得足够远，近似程度可以很高。但这却只能保证所述方法的使用条件，并没有将区域背景场去掉。在数据处理过程中，区域背景场总是要在一定的阶段予以去掉的。但由此带来的误差，应归于所使用的确定区域背景场的方法，与这里所介绍的方法系统无关，本书不予讨论。
（3）所得拉氏方程基本级数解的系数，是在空间一定范围内的，并且测量往往来自正常场的实测场。这样，所得级数作为定义的上半空间的位场函数本身就是近似的。在空间内，随着离实测的线和面的距离增大，其近似程度将降低，离得愈远，降低愈多。从实测量经过转换计算而得的量，随着转换中计算过程的增加，其近似程度将降低。计算的过程愈复杂，降低得也愈多。
还应该注意到，当为了减少误差而采取某种措施时，上述（2）和（3）互相矛盾，难于兼顾。从（2）来看，零值边界自然取得离测区边部远些为好。但这样一来，在整个取定的空间范围内，实测区域所占的比例就要减低，从而使位场函数的近似程度减低。反之，则又将使零值边界假定的近似程度降低。在进行具体计算时，只能两者兼顾，求取折衷。
3.2.3.2 实际达到的精度
上述方法系统的误差来源，其中一部分，如计算数学中的近似性，是任何数据处理方法都有的，另一部分则是所述方法系统所特有的。但这绝不是说该法的计算精度差。它的精度是较高的，其任何一种功能均能完全满足数据处理的要求。只要给出足够的计算量，一定可以得到所需的计算精度。
所得的基本解是一个级数，在不超过实测点数目的范围内，显然项数愈高，计算精度也愈高。取级数项数超过测点数是不行的，因为系数要实测数据确定。由确定系数的方程组可知，当测点数少于项数即系数的数目时，会发生方程式少于未知数的情况，其解是不能唯一确定的。实际情况是：为了满足数据处理的要求，所取级数项数可以比测点数小得很多。对一般变化比较平缓的异常数据，级数项数和测点数的比有1：10左右就可保证边部亦有相当高的计算精度。对于梯度大，如强地形的磁测资料等情况，用1：3 左右也够了。再提高计算精度的潜力还相当大，这是问题的一个方面：级数项数取得愈高，精度也愈高。
但这里所说的精度愈高，是仅只对用以确定级数系数的实测数据所在的邻近区域而言的。对比较远的区域，其精度则反而要降低。这是由于级数项数取得愈高，对实测数据的模拟就愈细，因而函数也愈适应于实测数据区域的近旁。其远处的适应能力就反而要有所降低了。这在所述方法系统中是经过多次试算证明了的。在自然界，这样一种规律是有普遍意义的，对局部的过分适应，往往要以对于全局适应性的降低为代价的。
系数最初是由实测数据确定的。所取级数项数合适与否，一个客观存在的检验标准是：当系数已确定，因而位场函数也已确定之后，用此函数反算用以确定系数的实测数据。这不是拟合，而是检验，目的是看一下级数项数取得合适与否。在设定的精度要求下，级数的项数应取得愈少愈好。幸运的是对于具体的数据，在达到一定的精度之后，再增加级数项数，所带来精度上的增加是微小的。例如计算球体模型形成的50×50的规则网，其边部的值在2nT左右，极值为300nT的测区进行计算，取级数项数为100和2500的计算结果的差别大部分在0.1nT左右，全部在1nT以下。原因是取级数项数为100时，误差已在1nT以下了。可见盲目的增加级数项数，实在是大可不必。
3.2.3.3 接图
在重磁数据处理中，所有方法程序在计算时最大的误差都发生在边部，这样使将一个大面积的资料进行分块处理发生困难。因为分块处理后必须进行接图，而要接图则必须边部的误差很小。这对于计算精度是一个非常集中而严峻的问题，特别是接图线划在异常峰值附近时难度更大。
这里介绍的方法用于确定系数的数据所在的空间位置附近，是完全能够接图的。在数据处于同一平面而梯度不是太大时，数据的分块完全可以接图线为界。在地形起伏较大，或异常梯度大的情况下，接图分块也只需从接图线向外超出一、二个点、线就行了，而且在接图线经过异常峰值时亦可以做到。这对系统的精度作了一个综合性的显示，应该说是相当高的。
