（地球探测与信息技术专业论文）三维磁场曲化平研究与编程（边界单元法）.pdf

桂林工学院硕士学位论文 摘要 在实际的野外磁法勘探中，三维起伏地形不可避免，地形影响使三维磁场等 值线变得复杂，或者使异常曲线特征点变得不明显，直接影响和降低了磁法勘探 的地质效果。因此。为突出有用的磁异常信息。弄清楚地形对三维磁场异常曲线 的影响，并采取有效的措施即所谓地形改正的办法是山区开展磁法勘探急需解决 的实际问题。本文应用边界单元法对三维磁场进行三维任意地形改正。 在此过程中，应用边界单元法数学理论基础，导出了利用边界单元法进行三 维磁场曲化平的计算表达式。根据需要编制了列主元高斯消元法、三维磁场正演 程序、三维磁场曲化平软件。由于边界单元法中的系数矩阵是对角占优的正定矩 阵，因此采用列主元高斯消元法求解方程组。由于系数矩阵非奇异，而且是正定 正定矩阵，用列主元高斯消元法具有速度快及良好的数值稳定性。编制的正演程 序可以模拟多个规则六面磁性体在任意三维地形上产生的磁场，采用对话框界面， 清晰明了，其主要作用是为曲化平计算提供模型数据和进行误差分析对比。通过 与模拟实验结果对比，证明是正确的。三维磁场曲化平软件利用v i s u a lc + + 6 0 编制，软件编制中采用了数据库进行输入输出处理，大大减少了对原始数据文件 格式的依赖，降低了前期数据处理的难度，提高了效率。为了兼顾误差和速度。 软件中边界单元法采用三角单元剖分，然后用一次插值进行单元积分计算，非奇 异积分采用高斯7 点积分，而奇异积分采用的自己推导的解析解，实算表明是正 确的。此外三维磁场曲化平软件采用了w i n d o w n s 界面，运算时能实时给出相关信 息。软件运行时，使用者能时刻了解运算的进展情况，及时处理问题。为了提高 三维磁场曲化平软件的实用性，在软件中加入了处理边界效应的程序代码。实算 表明，对边界效应的处理是有效的，经过处理后的边界点延拓值误差减小到1 0 左右，网格中间点的延拓值误差减小到5 左右。 经过对多个模型的曲化平计算，将曲化平延拓后的平面等值线与平面理论等 值线对比，证明曲化平延拓能有效的恢复磁场在水平平面上的形态特点；而通过 对比边界改正前和改正后的边界剖面线及平面等值线，证明了边改方法的正确性 和有效性。所编的正演及曲化平软件经多次应用，运行稳定可靠，说明编制的程 序是正确的。 关键词：三维曲化平、三维磁场、边界改正、v c + + 、边界单元法 桂林工学院硕士学位论文 a b s t l 认c t i nt h ea c t u a lo p e nc o u n t r ym a g n e t i cp r o s p e c t i n g ，t h et h r e ed i m e n s i o n a l a c c i d e n t e dr e li e fi si n e v i t a b l e ，t h et e r r a i ni n f l u e n c ec a u s e st h et h r e e d i m e n s i o n a lm a g n e t i cf i e l de q u i v a l e n t1 i n et ob e c 伽ec o p l e x ，o rc a u s e s t h eu n u s u a lc u r v ec h a r a c t e r i s t i cp o i n tt ob e c o m en o to b v i o u s ， d i r e c t l y a f f e c t e da n dr e d u c e dt h em a g n e t i cp r o s p e c t i n gg e o l o g i c a le f f e c t t h e r e f o r e ， f o rt h ep r o m i n e n tu s e f u lm a g n e t i s mu n u s u a li n f o r m a t i o n ，c l a r i f i e st h e t e r r a i nt ot h et h r e ed i m e n s i o n a lm g n e t i cf i e l du n u s u a lc u r v ei n f l u e n c e ， a n da d o p t st h em e a n sw h i c ht h ee f f e c t i v em e a s u r en a m e l ys o c a l l e dt e r r a i n c o r r e c t si st h em o u n t a i n o u sa r e ad e v e l o p m e n tm a g n e t i cp r o s p e c t i n gu r g e n t n e e ds o l u t i o na c t u a lp r o b l e m t h i sa r t i c l ea p p l i e st h eb o u n d a r yc e l l t e c h n i q u e 胃a yt oc a r r yo nt h et h r e ed i m e n s i o n a lr a n d o mt e r r a i nc o r r e c t i o n t ot h et h r e ed i m e n s i o n a lm a g n e t i cf i e l d i nt h i sp r o c e s s ， a p p l i e stb o u n d a r yc e l lt e c h n i q u e w a yo p e r a t o rs t u d y r a t i o n a l e ， d e r i v e dh a sc a r r i e do nt h r e ed i m e n s i o n a lm a g n e t i cf i e l d c o n v e s i o no fc u r v e d s u r f a c ed a t a i n t ob o r i z o n t a l 一p l a n ed a t au s i n gt h e b o u n d a r yc e llt e c h n i q u e a yt h ec o 皿p u t a t i o ne x p r e s s i o n a c c o r d i n gt on e e d e d t oe s t a b li s ht h er o wp r i n c i p a le l e 加e n tg a u s se l i m i n a t i o n ， t h et h r e e d i m e n s i o n a lm a g n e ti cf ie l dh a sb e e nd e v e l o p i n gt h ep r o c e d u r e ，t h r e e d i m e n s i o n a lm a g n e t i cf i e l do n v e s i o no fc u r v e d s u r f a c ed a t ai n t o b o r i z o n t a l p l a n ec l a t at h es o f t w a r e b e c a u s ei nt h eb o u n d a r yc e ll t e c h n i q u e w a yc o e f f i c i e n t 髓t r i xi st h eo p p o s i t ea n 9 1 eo c c u p i e ss u p e r i o r l y i sd e c i d i n gt h em a t r i x ， t h e r e f o r eu s e st h er o wp r i n c i p a le l e m e n tg a u s s e 1i m i n a t i o ns o l u t i o ns y s t 伽o fe q u a t i o n s b e c a u s ec o e f f i c i e n tm a t r i x n o n s t r a n g e ，m o r e o v e ri st h ec o r r e c t i o nd e c i d e st h em a t r i x ，h a st h es p e e d w i t ht h er o wp r i n c i p a le l e m e n tg a u s se l i m i n a t i o nq u i c ka n dt h eg o o dv a l u e s t a b i l i t y t h ee s t a b l i s h e n ti sd e v e l o p i n gt h ep r o c e d u r et ob ep o s s i b l e t os i m u l a t et h em a n yr u l es i xm a g n e t i c st h em a g n e t i cf i e l d - h i c hp r o d u c e s w i11 f u ll yi nt h et h r e ed i m e n s i o n a lt e r r a i n ，u s e st h ed i a l o gb o xc o n t a c t s u r f a c e ， t h ec l e a rp e r s p i c u i t y ， i t sm a i nf u n c t i o ni sc a l c u l a t e sf o r o n v e s i o no fc u r v e d s u r f a c ed a t ai n t ob o r i z o n t a l p l a n ed a t ap r o v i d e st h e m o d e ld i a t aa n dc a c i r i e so nt h ee r r o ra n a l y s i sc o n t r a s t t h r o u g h2 u l dt h e m o d e l li n gr e s u l tc o n t r a s t ， t h ep r o o fi sc o r r e c t t h r e ed i m e n s i o n a l m a g n e t i cf i e l dc o n v e s i o no fc u r v e d s u r f a c ed a t ai n t ob o r i z o n t a l 一p l a n ed a t a v 桂林工学院硕士学位论文 t h es o f t w a r eu s e sv is u a lc + + 6 oe s t a b lis h m e n t 。 i nt h es o f t w a r e e s t a b l i s h m e n tu s e dt h ed a t a b a s et oc a r r yo ni n p u to u t p u tp r o c e s s i n g ， r e d u c e dg r e a t1yt ot h ep r i m a r yd a t ad o c u m e n t f o r m a td e p e n d e n c e ， r e d u c e d t h ep r e l i m i n a r yd a t ap r o c e s s i n gd i f f i c u l t y ， e n h a n c e dt h ee f f i c i e n c y i n o r d e rt og i v ed u a la t t e n t i o nt ot h ee r r o ra n dt h es p e e d ， i nt h es o f t 霄a r e t h eb o u n d a r yc e l lt e c h n i q u e w a yu s e st r i a n g l eu n i tc u t t i n gi nh a l f ， t h e n u s e s a ni n t e r p 0 1 a t i o nt oc a r r yo nt h eu n i ti n t e g r a lc o m p u t a t i o n ， t h e n o n s i n g u l a ri n t e g r a lu s e st h eg a u s s7i n t e g r a l s ，b u tt h es i n g u l a ri n t e g r a l u s e so - ni n f e rt h ea n a l y t i cs o l u t i o n ，t h er e a l i t yc a l c u l a t e di n d i c a t e di s c o r r e c t i na d d it i o nt h r e ed i m e n s i o n a lm a g n e t i cf i e l dc o n v e s i o no f c u r v e d s u r f a c ed a t ai n t ob o r iz o n t a l 一p l a n ed a t at h es o f t w a r eh a su s e dt h e w i n d o w n sc o n t a c ts u r f a c e ，w h e nt h eo p e r a t i o nc a nr e a l t i m eg i v et h er e l a t e d i 力f o r m a t i o n r h e nt h es o f t w a r em o v e s ，t h eu s e rc a nt h et i eu n d e r s t a n d i n g o p e r a t i o np r o g r e s ss i t u a t i o n ， d e a lw i t ht h ei s s u ep r o 叩t l y i no r d e rt o e n h a n c et h r e ed i m e n s i o n a lm a g n e t i cf i e l dc o n v e s i o no fc u r v e d s u r f a c ed a t a i n t ob o r i z o n t a l 一p l a n ed a t at h es o f t w a r et h eu s a b il i t y ， j o i n e di nt h e s o f t w a r eh a sp r o c e s s e dt h eb o u n d a r ye f f e c tt h ep r o c e d u r ec o d e t h er e a l i t y c a l c u l a t e si n d i c a t e dt h a t ，t ob o u n d a r ye f f e c tp r o c e s s i n gi se f f e c t i v e ， a f t e rt h ep r o c e s sp r o c e s s i n gb o u n d a r yp o i n ta n a l y t i cc o n t i n u a t i o nv a l u e e r r o rr e d u c e st o1 0 a b o u t ，m e a s u r e dt h en e ti n t e r m e d i a t ep o i n tt h e a n a l y t i cc o n t i n u a t i o nv a l u ee r r o rr e d u c e st o5 a b o u t p a s s e st h r o u g ht o t h em a n y o d e lt u n ee v e nc o 叩u t a t i o n ，a f t e ro n v e s i o no fc u r v e d s u r f a c ed a t a i n t ob o r i z o n t a l p l a n ed a t aa n a l y t i cc o n t i n u a t i o np l a n ee q u i v a l e n t1 i n ea n d p l a n et h e o r ye q u i v a l e n tl i n ec o n t r a s t ， t h ec a r dm i n gd y n a s t yd r 锄ae v e n a n a l y t i cc o n t i n u a t i o nc a nt h ee f f e c t i v er e s t o r a t i o nm a g n e t i cf i e l di n l a t e r a lp l a n em o r p h 0 1 0 9 i c a lp r o p e r t y b e f o r eb u ta f t e rc o n t r a s t st h e b o u n d a r yt oc o r r e c ta n dt h ec o r r e c t i o nb o u n d a r yp r o f i l el i n ea n dt h ep l a n e e q u i v a l e n t1 i n e ， h a dp r o v e ns i d ec h a n g e st h em e t h o dt h ea c c u r a c ya n dt h e v a l i d i t y a r r a n g e si sd e v e l o p i n g 锄dc o n v e s i o no fc u r v e d s u r f a c ed a t ai n t o b o r i z o n t a l 一p l a n ec i a t at h es o f t w a r ea f t e ra p p li e sm a n yti m e s ，m o v e m e n t s t a b l er e l i a b l e ，t h es h o w i n ge s t a b l i s h m e n tp r o c e d u r ei sc o r r e c t k e y - d r d ：t h r e ed i m e n s i o n a lc o n v e s i o no fc u r v e d s u r f a c ed a t ai n t o b o r i z o n t a 卜p l a n ed a t a ， t h r e ed i m e n s i o n a lm a g n e t i cf i e l d ， b o u n d a r y c o r r e c t i o n ，v c + + ，b o u n d a r yc e l lt e c h n i q u e w a y 桂林工学院硕士学位论文 研究生学位论文独创性声明和版权使用授权说明 独创性声明 本人声明：所呈交的论文是我个人在张小路研究员指导下进行的研究工作及 取得的研究成果。尽我所知，除了文中特别加以标注和致谢的地方外，论文中不 包含他人已经发表或撰写过的研究成果，也不包含为获得桂林工学院或其它教育 机构的学位或证书而使用过的材料。对论文的完成提供过帮助的有关人员已在论 文中作了明确的说明并致以了谢意。 学位论文作者( 签字) ： 签字日期： 版权使用授权说明 本人完全了解桂林工学院关于收集、保存、使用学位论文的规定，即：按照 学校要求提交学位论文的印刷本和电子版本；学校有权保存学位论文的印刷本和 电子版，并提供目录检索与阅览服务；学校可以采用影印、缩印、数字化或其它 复制手段保存论文；在不以赢利为目的前提下，学校可以公布论文的部分或全部 内容。( 保密论文在解密后遵守此规定) 学位论文作者( 签字) ： 指导教师签字： 签字日期： 桂林工学院硕士学位论文 第一章引言 对山区磁测资料曲化平一直是重磁位场理论中的个重要问题，目前已由二维发展 到三维，但三维曲化平虽经多方研究，目前仍处发展之中。起伏地形对磁异常形态特征 的影响主要有以下几种现象：1 地形起伏时，磁异常极值将发生偏移；2 由于地形影响， 会使对称变得不对称，使不对称异常变为对称。3 地形影响还会使简单异常变为复杂形 态。4 地形起伏变化，会使无负值的异常出现负值，有负值的异常在一定的范围内变得 无负值，或负值更大。以上讨论说明：当地形起伏时，沿地形线观测的磁异常与沿水平 面观测的磁异常相比，不论是异常的强度还是形态，都有很大变化。若将起伏地形上的 磁异常当作水平观测面的情况来处理，会产生很大的误差，以致发生解释错误。克服地 形影响的办法之一就是将起伏地形上的异常化为水平地形上的异常既“曲化平。 对二度体磁异常的曲化平方法从手算角度看主要有两类：1 、利用复变函数中的柯 西积分的量板法及保角映射计算方法n 1 ，利用保角映射原理n 1 的电模拟方法。此类方 法主要由2 0 世纪5 0 年代的黄树棠d 1 ( 1 9 6 6 ) 、雷林源1 ( 1 9 6 5 ) 等人进行了系统研究。 其中保角映射计算方法的基本原理是：在z 平面上位场的狄利希莱问题经过一个解析变 换( 映射) 后，在w 平面上仍为位场的狄利希莱问题。这样，就有可能把复平面( z 平 面) 上较复杂的边界值问题化为另一个复平面( w 平面) 上较简单的边界值问题，也就 是说可以把z 平面上曲面延拓问题变为w 平面的平面延拓问题，然后再把w 平面上延拓 结果搬回到z 平面上，就完成了曲化平的工作。由于此方法要求用一函数来近似地描绘 地形线，而有些不规则的地形不易描绘，并且其变换函数难以得到，计算又较麻烦，上 述方法的使用就很困难，此时较好的方法是用保角映射原理的电模拟方法，原理与上面 方法相同，所不同的是不用计算方法，而是用电模拟方法求模拟出地形线，水平延拓线 的映射；量板法5 “1 是一种曾经广泛使用的方法，但由于使用繁琐、效率和精度低， 现在已很少有人使用。2 、利用一个特定的函数来近似地表示地形上磁测曲线，如罗朗 级数法。6 0 7 0 年代，开始运用数值计算方法，二度体磁异常的曲化平及三度体磁异常 的曲化平方法都得到了极大发展。出现许多新的方法，主要有以下几种：1 、等效源偶 层法或等效源法7 l 盯。e i 玎il a 9 1 ( 1 9 7 3 ) 、b h a t t a c h a r y y a 与c h e n ( 1 9 7 7 ) 1 盯、n a k a t s u k a ( 1 9 8 1 ) 1 1 1 、h a n s e n 与m i y a z a k i ( 1 9 8 4 ) 1 2 1 ，以及周熙襄( 1 9 7 9 ) 、杜维本n 3 1 ( 1 9 8 2 ) 、 陈钟琦“钉( 1 9 8 3 ) 、侯重初1 6 1 ( 1 9 8 5 ) 、管志宁埔1 ( 1 9 8 5 ) 、王忠敏m ( 1 9 8 6 ) 、许连喜 n 盯( 1 9 8 8 ) 、王万银1 ( 1 9 9 1 ) 等利用位场理论将起伏地形的重磁异常换算单层或偶层 桂林工学院硕士学位论文 形式三维等效源，再由此等效源计算空间某一平面上的重磁场达到曲化平的目的，这类 方法通常称为等效源偶层法或等效源法，用等效源法作磁异常的曲化平方法在计算时不 会损失边界的测点信息，必要时还可以将异常适当外推，计算结果表明在一定的条件下 在场的变换上精度较高，并对起伏较大的观测面进行曲化平的处理有一定的效果，但也 存在根据磁场实际特征求得合适的等效源及其分布也是一个不易解决的问题；3 、级数 类方法，h a n a d e r s o nc o r d e l l 剐( 1 9 7 1 ) 、何玉辉2 ( 1 9 8 5 ) 、梁锦文阮1 ( 1 9 8 8 ) 、王 万银n 3 剐( 1 9 9 1 ) 等则把重磁场展成级数，求出级数的系数，再计算空间任意平面的 重磁场，由此来实现重磁资料的曲化平。级数类方法较多，如调和级数法2 钉，此方法 在地形较复杂条件下，效果较差，异常易出现震荡现象，为克服这一问题，进行正则化 处理可以起到一定的效果。4 、样条函数法，王硕儒、黎春华协1 ( 1 9 9 1 ) 利用样条函数 来处理积分方程中的导数的计算，也称样条函数法。一般选择双三次样条函数，利用样 条函数求导的观测场的垂向二阶和三阶导数，最后用泰勒级数计算出某一平面上的近似 值，达到曲化平的目的，此方法精度较高，速度也较快，主要的问题是延拓的高度有限， 误差随高度而增加。5 、边界单元法，有限单元法：程振炎瞄1 ( 1 9 8 1 ) 、徐世浙伽剐阳1 ( 1 9 8 5 ) 利用边界元法、有限单元法开创了磁场曲化平新局面，用有限单元法进行磁场 曲化平时精度很高，但计算时数据量大，对计算机的要求高以及边界条件不好确定等问 题，而边界单元法进行计算时，数据量较小，速度快，原理简单，精度也较高，但边界 误差较大，此外当测点数据多时，对计算机的要求也较高。此外，国外关于曲化平的方 法3 1 1 也很多，除了上面介绍的方法外还有一些简易而有效的方法，如棋盘法( c h e s s b o a r d ) ，由美国地调局l c o r d e l l 脚1 提出，优点是可计算大批量数据，速度快，缺点是 理论不严格；位函数拟合面延拓法由美国地调局y j s g r a u r o h 提出，此方法采用逐次 剥离的迭代算法，不断改善拟合度，直到满足给定的精度为止，此法的适应性强且计算 中可调节算子提高计算速度，但计算耗时较长。 磁场曲化平能有效的消除地形对磁异常的影响，恢复了磁场在平面上的真实形态， 有助于解释人员运用通常的解释方法，正确的认识地下磁性体的分布。目前，曲化平方 法很多，取得了很好的效果。例如：等效偶层法是目前理论上比较完善的曲化平方法， 但是无法处理大面积、大数据量的实际资料；同样，有限单元法也是一种非常有效的方 法，但有和等效偶层法相同的缺点。此外还有许多方法，它们有的计算速度快，但误差 大，有的速度慢但误差小；因此寻找一种又快又精的方法是大家的共同目标。通过多年 的研究，目前曲化平方法已由模型模拟发展到数值模拟，并随着计算机的发展而飞速发 2 桂林工学院硕士学位论文 展。理论方面由于边界单元法和有限单元法的应用为曲化平方法的发展提供了强有力的 手段。此外，各种对曲化平方法的改进仍在进行中，主要是在速度、精度方面进行改进， 相信不久曲化平方法将在磁法勘探中广泛使用。 本文在前人工作的基础上，为了进一步完善磁场曲化平的边界单元方法主要作了以 下研究工作： 1 改进了用边界单元法实现三维磁场曲化平的方法。利用三角单元自动剖分 网格极大的减少了原始数据的输入工作量；利用一次线性插值提高了计算 精度。 2 利用拉格朗日二维插值函数及差分对三维磁场边界混合向外插值，实现用 。外推法 对边界单元法边界效应进行混合改正的目的，大大降低了边界 点及内部各点的相对误差。 3 边界单元法中的奇异积分是编程中的难点，为了提高了计算精度和速度， 本文对三维边界单元法中一次插值时的奇异积分做了推导，求出了其解析 解，简化并加快了奇异积分的计算，计算结果表明完全正确。 4 原始数据的输入和结果的输出采用数据库技术，减轻了原始数据的整理和 输入工作，降低了对初始数据文件格式的要求。同时，将结果和原始数据 放在同一数据库中，实现了输出和输入的统一，整理结果时方便简洁。 5 采用v c + + 编制了w i n d o w s 程序界面的三维磁场曲化平软件和三维正演程 序。曲化平软件简洁、易于操作和维护，运行时提供多重反馈信息，使计 算过程清晰明了。编制的三维正演程序可以计算任意多个均匀磁化三度体 在任意位置时在空间任意地点磁场值( 包括各分量及总值) ，为曲化平的研 究提供了强有力的保证。 3 桂林工学院硕士学位论文 第二章三维磁场曲化平的理论基础 磁场是位场，位场延拓在磁测资料的解释中有重要的应用。根据地表实测的磁场值计 算上部空间的磁场值，称为向上延拓；计算下部空间的场值，称为向下延拓3 3 1 “m 轴。 利用延拓的方法可以将起伏地形上的磁场值向上延拓到水平平面上，从而达到曲化平的 目的。对于位场，边值问题是非常重要的，磁场的边值问题如下。 2 1 三维地形模拟与其边值问题 一三维地形模拟 在地球物理中，场源一般位于地下，在地表上部，没有场源，位函数u 满足拉普 拉斯方程2 蚰 v ：“：生生+ 生竺+ 生兰：o 在无源空间 ( 2 1 - 1 ) a x 2 a y 。8 z 2 对于磁力勘探来说，在地面测量的是场，不是位。场是位在某一方向的负梯度，但是经换 算3 可知，对于磁场磁感应强度b 及其分量同样有 vz b ：生旦+ 生旦+ 坐旦：o 在无源空间 ( 2 1 2 ) 8 x t 8 y t a z 2 因此，我们可以将场强看作位，都用u 表示。为了能确定边界条件，必须先建立一个三维 模拟地形，并且给定地形上的u 。为此我们给出如下模拟地形，在一起伏的有限的区域 q 内，地表如图2 1 所示。 4 桂林工学院硕士学位论文 各测量点在网格交点处，场源位于地表下，地表上无场源，以rs 代表地表边界，rs 向周围延伸至足够远，用r 二表示无穷远。有了模型后，我们就可以确定边界条件了。 二边值问题3 8 1 地形上部区域q ，由地表边界rs 及无穷远边界l 所包围。三维位场u 向上延拓的 边值问题是： vz 材：塑+ 塑+ 塑；o q ， 8 x 2 a y l 8 z 2 ( 2 1 3 ) “= o ，y ，z ) n ( 2 1 4 ) 甜2 素 l ( 2 l 5 ) 其中f ( x ，y ，z ) 是已知函数；在r 帕上，位场与场源的距离r 一次方成反比，也可以假 定与距离的二次方成反比。已证明啪1 ，在l 上的场值对q 中的场值的影响可以忽略m 1 。 三边界积分方程嚣1 对于起伏地形可取空间任意一点的基本解为：缈= 1 ( 4 万，) ，其中r 为空间任意 点p 到边界rs 上任意点的距离，将u 和缈代入三维格林公式有： l ( 以v2 9 一缈v2 州q = k 。( 撑署一伊d r ( 2 l6 ) 对( 2 1 6 ) 式的左侧积分有： 上( 扰v2 缈一缈v2 越n = 上( 一口万( p ) 一缈o 油舱 2 一上 万( p ) 施 国o = 一吉甜p 4 覆 对( 2 1 6 ) 式的右侧积分有： 鼻。c 甜等一缈筹= j r 。t 詈掣一去掣 在r 上，有r r 所以 靠。( ”等一缈等= 。 ( 2 1 7 ) 由此可得到公式： 5 桂林工学院硕士学位论文 等“，= j r ，r ”掣+ 署寺 亿m ， 当p 点在rs 上时，彩。为p 点对区域q 的立体角，用边界单元法解该积分方程可求得 锄锄。当p 点在q 内部时，有缈。= 4 万；将边界l 上的u 和加锄代入上式就可求得 区域q 内任意一点p 的场值 讨论：由于rs 和r 。构成一个封闭的空间区域，因此rs 和r 。之间必定有分界线。 因为分界线上的点不一定有r r ，所以为了使( 2 1 7 ) 式成立，u 和锄锄应该为零， 否则( 2 1 7 ) 式将有误差。实际工作中，勘测面积是有限的，往往不能满足上述要求， 既在现实中，r 。上的积分一般不为零。因此由此产生的误差将是不可避免的。但一般 情况下情况下，对测区中间的测点的误差较小。在后面的误差分析中将详细讨论。 2 2 三维磁场拉普拉斯方程的基本解3 7 1 拉普拉斯方程的基本解是在将边值问题转变为积分方程过程中一个非常重要的工 具，为此本节将对其基本解进行说明。 2 2 1 基本解的定义 线性微分算子l 对函数u 进行某种微分运算，构成一个微分方程 l ( u ) = o 用q ( x ，y ，z ) 表示动点，p ( x o ，y o ，z o ) 表示固定点。若某函数1 i r 的微分运算，得到一个以p 点为中心的负6 函数，即 l ( v ) = 一6 则称1 i r 为l ( u ) = 0 的基本解。 2 2 1 三维拉普拉斯方程的基本解 三维拉普拉斯方程 v ：堕+ 塑+ 生：o 8 x t a y t 8 z l 的基本解是 l 伊2 石万 ( 2 2 1 ) ( q ，p ) 经过l ( 2 2 2 ) ( 2 2 3 ) ( 2 2 4 ) 式中r 是三维空间中某点p 到任意点q 的距离( 见图2 ) 。基本解的求解过程如下。 6 周2 2 以p 点为坐标原点，将方程 v 2 缈= 一万( p ) ( 2 2 5 ) 写成球坐标形式，讨论它的球对称解，此时，9 只与r 有关。当r o 时，万( p ) = 0 。 此时( 2 2 5 ) 可写成 v 2 驴= 吉昙( ，2 警) - - 却) ：0 ( 2 2 6 ) 上式积分后可得其通解 矽：兰+ 曰 ( 2 2 7 ) 式中，a ，b 是积分常数。显然，伊满足 v 2 矽= ，r = 0 ( 2 2 8 ) 此外，对全部区域q ( 如图2 2 所不) 硐 l v 2 碑= 一l 万( p 闺一1 在图2 中作以p 点为中心得小球e ，根据( 2 2 6 ) 和( 2 2 7 ) 式有 l v2 徊q = 一上v 2 缈d q2p v 缈擒 = 垂。( v 九订= 垂。争一垂。争 一等协一乒4 矿一4 蒯一一7 孵12 7 彻r 一叶翮 由( 2 2 9 ) 和( 2 2 1 0 ) 知a = l ( 4 ) ， 1 回= 一 7 4 万r ( 2 2 9 ) ( 2 2 1 0 ) 再令b = 0 得三维拉普拉斯方程得基本解为 ( 2 2 1 1 ) 7 桂林工学院硕士学位论文 因为磁场的磁感应强度在无源三维空间满足式( 2 2 3 ) ，因此基本解也是磁场在无源三 维空间的基本解，故在( 2 1 6 ) 式到( 2 1 7 ) 式的推导中引用了该基本解。 8 桂林工学院硕士学位论文 第三章利用边界单元法进行三维磁场曲化平的方法 由第二章的公式( 2 1 7 ) 式可知，如果能求解该积分方程就能得到区域q 内任意一 点p 的场值从而用延拓的方法实现曲化平的目的。求解该方程的数值方法可以用有 限单元法，也可用边界单元法等。由于边界单元法简单可靠，所占计算机内存较小，速 度快，因此本文采用边界单元法。 3 1 边界单元法的基本方法 边界单元法的基本方法是：用三角单元对地表边界rs 进行剖分( 见图3 1 ) ，使每 个单元的地形近似平面，然后在三角单元内用线性等参单元进行一次插值3 钉阳1 。其中 剖分方法和过程比较重要，好的剖分方法能简化后续计算。 纽距 9 桂林工学院硕士学位论文 自动剖分对于编程计算非常重要，将节点号按上图排列，能实现自动剖分。每个三 角单元含有三个节点( 如图3 1 ) ，其排列顺序为逆时针方向。例如：1 号单元的节点编 号顺序为( 1 ，7 ，2 ) ；又如5 0 号单元的节点编号为( 3 0 ，3 5 ，3 6 ) 。按此排序后就可进 行网格自动剖分，具体过程如下： 第一步：按图3 1 排列测量数据( 包括t 和测点的三维坐标，三维坐标可以是设 计坐标也可是g p s 坐标) ，既一条测线接一条测线的输入数据。 第二步：自动剖分。为了能找出剖分过程的规律，我们先将上例图3 1 中的部分单 元的节点号排列如下： 第1 线和第2 线间的单元第2 线和第3 线间的单元 单元号节点号单元号节点号 11721 171 38 22781 281 31 4 32831 381 49 43891 491 41 5 5394 1 591 51 0 6491 0 1 61 01 51 6 741 051 7l o1 61 1 851 01 11 8l l1 61 7 951 161 9l l1 71 2 1 061 11 22 01 21 71 8 从上面的数据中我们可以发现： 1 、第一条测线和第二条测线间的单元的第一个节点号的排列规律是：设单元号为 m ，其第一个节点的编号设为n 1 。当m 为奇数时其第一个节点的编号n 1 等于其前一个 单元的第一个节点的编号，例如第三单元的第一个节点的编号为“2 就等于第二单元 的第一个节点的编号“2 一；当m 为偶数时其第一个节点的编号n l 等于其前一个单元的 第一个节点的编号加“1 一，例如第二单元的第一个节点的编号为“2 刀就等于第一单元 的第一个节点的编号“l + 1 ，既奇数不变，偶数加l 。 2 、第一条测线和第二条测线间的单元的第二个节点号的排列规律是：当m 为奇数 时，其第二个节点的编号n 2 等于同单元第一个节点编号n l 加上每条测线的测点总数， 既n 2 = n 1 + ( 测点总数测线) 。例如：第五单元，其第一个节点编号为3 ，而每条测线 l o 桂林工学院硕士学位论文 的测点数为6 ，所以第五单元的第二个节点编号为3 + 6 = 9 ；当m 为偶数时， n 2 = n l + ( 测点总数测线) 一1 。 3 、第一条测线和第二条测线间的单元的第三个节点号的排列规律是：当m 为奇数 时其第三个节点的编号n 3 = n 1 + 1 ；当m 为偶数时，n 3 = n 1 + ( 测点总数测线) 。 4 、第二条测线和第三条测线间的单元的第一个节点号的排列规律同第一条测线和 第二条测线问的单元的第一个节点号的排列规律，只是将结果加上( 测点总数测线) 即可。如第l l 单元的第一个节点号是“7 ，就是将第一单元的第一个节点号“1 加上 ( 测点总数测线= 6 ) 得到的。至于其余节点号类推即可。 由上我们就可以通过计算机编程实现自动剖分，所需输入的数据主要是每条测线上 的总测点数( 注：上述规律要求每条测线上的总测点数相同) 及全部的测点总数。此外， 剖分时并不要求测线是直线，测点能按逆时针排列即可，这对实际测量的计算是有利的， 不需重新网格化，可以利用原始数据进行计算。 关于单元内的插值运算之所以选择一次插值主要是考虑到如下两点：l 、利于剖分， 因为此时的节点为三角单元的顶点；2 、提高精度，在徐世浙教授关于边界元三维延拓 的论述中，举例说明了零次插值时计算误差约为1 1 ，因此使用一次插值能完全满足 精度要求；3 、保持运算速度，随着插值次数的增高将增大计算机内存的需求和增加运 算时间。应用一次插值能合理兼顾速度与精度的关系。下面是插值运算的具体过程。 首先，通过将形函数与单元上各节点的值相乘得到单元内各点的值。 甜= ：，甜r ( 3 1 1 ) 要= ，( 署) ， ( 3 l2 ) a 刀 一 、a 刀7 x ：y ，z ， ( 3 1 3 ) y = ，y ， ( 3 1 4 ) z = ，z ， 1 = j ，k ，m ( 3 1 5 ) 其中u 是( x ，y ，z ) 点的值，n j ，m ，n - 是形函数，其计算公式为： nj = 丢 xxkx m yyk y 嘛 zz t z m 小扰 xxm yym zzm 桂林工学院硕士学位论文 七= 吉 x j x k x yj yky z j z k z l 1 ，= 一 6 x j ) i z j 其次，将( 3 1 1 ) ( 3 1 2 ) 代入( 2 1 7 ) 式得单元积分： 甜等碧dr ： 4 万，。 一。”。 署d r = 州等) ，a 刀4 万，二一“、a 刀 7 l = j ，k ，m 式中 厶= 上，争 最后将各单元的积分相加得： d a = k n | 争 x 孟 y 七 z 孟 毛j r 材掣d r2 ；兀 ，。莓凡甜， 莩i 筹一手幽( 等) - 李州驯， x 历 ym z 埘 ( 3 1 6 ) ( 3 1 7 ) ( 3 1 8 ) ( 3 1 9 ) ( 3 1 1 0 ) 将( 3 1 1 0 ) 式和( 3 1 1 1 ) 式代入( 2 1 8 ) 式可得到由n 个节点构成的线性方程组： ( 国4 万一，) 甜= d ( a 甜a 甩) ( 3 1 1 2 ) 其中 = d i a g ( ( i ) 。) ，f = ( f 。j ) ，d = ( d ；j ) ，u = ( u ；j ) 7 ，a 甜归刀= 【( a 甜归刀) ，】丁 i ，j = 1 ，2 ，3 ，n 求解( 3 1 1 2 ) 式得到锄锄，然后将u 和a l | 锄代入由( 2 1 8 ) 式构成的方程组 “p = ，铭+ d ( a 甜a 刀) 既可得地形上部任意点的甜p 。 1