2. 地下信息探测技术与仪器教育部重点实验室(中国地质大学, 北京), 北京 100083
2. Key Laboratory of Geo-detection(China University of Geosciences, Beijing), Ministry of Education, Beijing 100083, China
位场延拓在重力、磁测等位场资料的处理和解释中具有广泛用途.因为向下延拓平面更接近于场源,场在场源附近的特性对于反演有重要意义,所以对位场向下延拓方法的研究一直是重磁数据处理方法研究中的热点之一.方法研究沿空间域和频率域两个方向进行,但数学原理相同,都是由外部狄里希莱问题求出向上延拓的严格解,再导出向下延拓的解,并根据计算方法的不同,发展了多种多样向下延拓的具体方法(Huestis et al.,1979,1998;Fedi et al.,2002;梁锦文,1989;王兴涛等,2004;徐世浙,2006;陈生昌和肖鹏飞,2007;刘东甲等,2009;闫辉等,2010).但现有的延拓方法在向下延拓接近场源时,都会产生剧烈震荡,计算都十分不稳定.杨文采(1997)曾从数学上证明 “位场向下延拓的最大深度是接近最浅的奇点,再深时解不存在”.因此目前所有的向下延拓方法,理论上都不能越过场源.
位场向下延拓不能越过场源的根本原因,是因为位场在场源奇点上趋于无穷大,位场函数在有源空间里不是解析函数.针对这一问题,本文在复数域里定义了一个新函数,该函数在场源奇点上是有限值,并在包含场源的空间里解析,解决了位场在有源空间不解析给向下延拓造成的困扰.提出了一套新的延拓方法,可以越过场源,不会产生震荡,在理论上实现了位场在有源空间的解析延拓.
2 方法原理二维水平圆柱体的复磁场强度F的复数表达式(吴功建等,1980;张秋生,1985):
F的实部和虚部分别是磁场的垂直分量Z和水平分量H.
二维水平圆柱体的复磁位V的复数表达式:
(1)、(2)式中M为水平圆柱体的复磁矩;
ζ0=x0+iy0为 水平圆柱体柱心复坐标;
z=x+iy 为空间任意点的复坐标(见图 1).
定义一个新函数R:
以(1)、(2)式代入R:
因为复磁位V和复场强F在无场源区是解析函数,解析函数的商也是解析函数(梁昆淼,1965;张秋生,1985),所以R在无场源区域是解析的.在场源处z=ζ0,R=z-ζ0=0 是有限值.很容易证明R满足柯西-黎曼方程.
则:
所以R在含有场源的有限大空间里(z≠∞)解析,可以使用复变函数论中的柯西公式根据解析函数在边界l上的函数值f(l),计算出在闭区域内的函数值f(z),实现解析延拓(梁昆淼,1965;张秋生,1985).
设水平圆柱的中心坐标 ζ0=20+20i,已知其在半径50的圆环边界上的R函数值:
按照以上推导,根据柯西公式由边界上的值计算圆内任意z点的R值:
实际计算得到R函数值,绘出R的等值线如图 2.
(1)图 2结果表示,可以根据边界值计算出圆内所有点上R的函数值,即R函数可以在有源空间内解析延拓.R在圆柱体中心的 ζ0 点有极小值,R等值线为包围场源中心 ζ0 的同心圆.随着远离场源,R逐渐增大.
(2)由(4)式知,R等于空间观测点z与场源奇点 ζ0 之间的距离,在复平面上,R的长度可用两点复坐标之差来计算.R等值线图形是描述空间观测点与场源奇点长度关系的几何图形,R具有几何属性.
(3)因R是由复磁位V和复磁场强度F相除计算得出的,故R可以反映V和F的位场属性,比如R在空间任意点的函数值,应反映所有场源产生的场在该点的叠加效应.
(4)R的几何属性和反映的位场属性在以下的运算中都要使用,很难说清这个具有双重属性函数的物理意义.借用爱因斯坦在广义相对论中把引力和空间弯曲结合在一起的思想,我们把R叫做PFGR函数(Potential field geometric representation).
3 延拓计算方法实际工作中只能测量得到地面或空中的磁场强度,不能测出地下闭曲线上的磁场值.要使用柯西公式沿闭曲线积分进行解析延拓,需将直角坐标系中的地面直线映射到闭合圆周上才能进行.
3.1 进行三次保角映射第一次映射通过公式(8)将Z平面映射为W1平面:
第二次映射通过公式(9)将W1平面映射为W2平面:
第三次映射通过公式(10)将W2平面映射为W3平面:
三次映射坐标对应关系见图 3.
三次映射可合并为一个统一的正映射公式:
至此已把Z平面映射为W平面,Z平面上的有限长线段AC映射为单位圆周,根据解析函数论中的边界对应原理,Z平面上的A→B→C→D→A包围的三角形区域映射为W1,W2,W3上相应点包围的区域,在W(W3)平面上映射为单位圆的内域.
按照(11)式计算出Z平面AC线段上测点Z(x,0),与W平面上单位圆周上对应点ei φ.x与φ 的对应关系见图 4.
由图 4见映射前直线上均匀分布的测点x,映射为单位圆周上总体均匀分布的 φ 点,这对使用柯西公式进行数值积分是有利的.
由正映射公式(11)可推出由W到Z的反映射公式:
3.2 构建PFGR函数根据地面上的实测场值构建单位圆周上的PFGR函数.已知地面上的V和F,则地面上的R函数为
所以
经三次保角映射后,Z平面上的Z(x,0)点映射为W平面单位圆周上的Wz点,Wz=eiφ. Z平面上的 ζ0 点映射为W平面上的 Wζ0 点(见图 5):
则Z平面上的距离R映射为W平面上的距离WR:
WR是单位圆边界上的PFGR函数值.
3.3 解析延拓已知单位圆边界上的PFGR函数值,则W平面单位圆内任意点m的PFGR函数值可按照柯西公式计算得出
以上计算的流程框图如下(图 6):
设定水平圆柱的圆心坐标 x0=y0=-100,ζ0=-100-100i,2M=1000cgs,代入(1)、(2)式,其理论复磁位V和复场强F为
根据地面上的V和F,采取上述方法向下半空间解析延拓,延拓结果见图 7.
由图 7见,延拓计算很稳定,延拓出的场源奇点等于设定的理论值,延拓计算精度很高. PFGR形态为包含场源圆心的同心圆,与磁位V的形态类似.
4 任意形状二度体的解析延拓4.1 延拓方法原理物体的磁化是其内部磁偶极子定向排列的结果,磁性体的外部磁场则是其内部磁偶极子场的叠加,二维时则是偶极线磁场的叠加.因为水平圆柱体的外部磁场等效于偶极线的磁场,故任意形状磁性体的磁场就可用无限多个水平圆柱体磁场的叠加来等效表示.
式中,Δζi=z-ζi,ζi 为场源内奇点(水平圆柱体柱心)的坐标.
式中 Δ ζi2=(z-ζi)2.
考虑到均匀磁化时各个水平圆柱体的磁矩 Mi 相等:
则
(24)式中R是一系列(z-ζi)的组合,只要任意一个 Δζi=z-ζi=0则R=0,即任意形状的磁场源,当观测点坐标等于场源内任意一个奇点的坐标时,PFGR函数都等于0.而当所有的 Δζi≠0 时,|R|≠0. PFGR等值线表现为在所有的场源奇点上,R取极小值,不在场源奇点上,R值增大,任意形状磁性体的PFGR等值线性质与单个水平圆柱体场源时相同.
分析(24)式可见,式中所有自变量 Δζi=z-ζi 表示的都是测点至场源的距离,故R和单一圆柱体一样,表现的也是距离参数,只是很难求出如 x-ζi 那样的解析表达式.我们采用数值方法来计算.
4.2 任意形状二度体解析延拓的计算方法在每个测点x下方放一个等效圆柱体,规定该点的磁场只是由这一个圆柱体产生的,则该圆柱体的柱心位置就可由(14)式计算:
N个测点下有N个圆柱体,可以算出N个 ζi,求出各x测点对应的 ζi,就可以按照单一圆柱体源同样的方法进行解析延拓了.
这些圆柱体在测线下方很近的位置排列成一个密集的链条,只要测点足够密,这个链条产生的磁场就与磁偶层等效.垂直磁化磁偶层面上的磁场是垂直方向的(实测磁场可先化极),故规定x点的磁场只受其下方一个圆柱体的作用,在一定精度下是合理的.
求出的 ζi 是否正确需要验证. 理论上x轴各点的磁场应等于所有等效圆柱体场的叠加,设其磁位为V1,场强为F1,每个圆柱体的磁矩等于常数K,则:
将计算出的等效源磁场V1和F1和实测磁场V和F进行比较,如果二者在允许误差范围内,则认为求出的 ζi 是正确的,可以继续进行延拓;如果不相等,延拓的误差可能增大.
4.3 验证实例1——倾斜板状体的延拓设定倾斜板状体的上顶坐标为(50,-100),下底坐标为(100,-150),复磁矩2M=1000 cgs,先计算出其理论复磁位和复场强,再采用上述延拓方法进行延拓,结果如图 8.
在图 8中,图a的实线表示倾斜板状体的理论磁场V和F的模,标记×表示由等效源计算出的磁场V1和F1的模,由图可见二者吻合得很好.说明计算出等效源的ζ(i)是正确的.用这些ζ(i)进行解析延拓,延拓顺利地越过场源,PFGR等值线确定的极小值范围与倾斜板状体重合,等值线圆滑,没有出现震荡,延拓精度较高.
4.4 验证实例2——三角棱柱体的延拓设定三角棱柱体顶点坐标分别为(80,-100)、(60,-120)和(100,-120),复磁矩2M=1000 cgs,计算出理论磁场V和F,延拓结果如图 9.
在图 9中,图a的实线表示三角棱柱体的理论磁场V和F的模,标记×表示由等效源计算出的磁场V1和F1的模,由图可见,等效源计算出的磁场与三角棱柱体理论磁场吻合很好,向下延拓计算稳定,没有出现震荡,PFGR等值线确定的奇点位置与三角棱柱体重合,延拓精度较高.
5 延拓影响因素分析5.1 保角映射时如何选择测线长度的大小由图 3知,测线长度AC=2a的大小决定着向下延拓三角形的深度和边长.因为PFGR函数 R=z-ζ,如果a过大,R中的z相应也很大.前已证明R只在包含 ζ 的有限大空间里解析,z过大,ζ 提供的有用信息相对减小,延拓误差会增大;如果a过小,实测剖面长度太短,磁场参数V和F形态不完整,不能给出场源全部信息,误差也会增大.试算结果表明,以测点距离为单位1,a选择500较好.
5.2 测量误差的影响实际工作中磁场参数V和F都存在测量误差,我们对一个组合长方体模型的V和F添加10%的随机噪声,延拓结果见图 10.延拓结果表明,延拓方法的抗干扰能力较好,计算稳定,没有出现震荡,延拓出的奇点位置也很准确.
本文推导计算中用的是复磁位和复场强,验证实例中使用的是垂直磁化的复磁场.实际工作中一般只测量磁场强度ΔT,或重力场强的垂直分量Δg,这些参数都可以在频率域中通过快速傅里叶变换(FFT)互相转换,简单易行,本文不赘述.
6 结 论(1)针对现有方法向下延拓出现困难的根本原因是位场在有源空间里不解析,本文定义了一个新函数PFGR,其值等于复磁位和复磁场强度的商.理论推导和模型试验表明,PFGR函数在包含场源的有限大空间里解析,在场源奇点处为有限值,可以用柯西公式对PFGR函数进行解析延拓.
(2)利用PFGR函数的几何特性和反映的位场特性,提出了一套新的向下延拓的数值计算方法,经过理论模型试算表明,对于任意形状的二度体,可以实现在包含场源的下半空间的解析延拓,PFGR等值线圈定了场源奇点.计算过程稳定,没有产生震荡,计算速度很快.
(3)模型试算表明,延拓方法具有较高的精度,对测量中出现的随机误差也有较好的抗干扰能力.
致 谢 衷心感谢董树文研究员、姚长利教授、孟小红教授和张永谦博士在本文研究工作中所给予的支持.[1] | Chen S C, Xiao P F. 2007. Wavenumber domain generalized inverse algorithm for potential field downward continuation. Chinese J. Geophys. (in Chinese), 50(6):1816-1822. |
[2] | Fedi M, Florio G. 2002. A Stable downward continuation by using the isvd method. Geophysical Journal International, 151(1):146-156. |
[3] | Huestis S P, Parker R L. 1979. Upward and downward continuation as inverse problems. Geophysical Journal International, 57(1):171-188. |
[4] | Liang J W. 1989. Downward continuation of regularization methods for potential fields. Chinese J. Geophys.(in Chinese), 32(5):600-608. |
[5] | Liang K M. 1965. Methods of Mathematical Physics (in Chinese). Beijing:Higher Education Press. |
[6] | Liu D J, Hong T Q, Jia Z H, et al. 2009. Wave number domain iteration method for downward continuation of potential fields and its convergence. Chinese J. Geophys. (in Chinese), 52:(6):1599-1605,doi:10.3969/j.issn.0001-5733.2009.06.022. |
[7] | Huestis S P. 1998. The continuation inverse problem revisited. Geophysical Journal International, 133(3):705-712. |
[8] | Wang X T. 2004. A comparison of different downward continuation methods for airborne gravity data. Chinese J. Geophys.(in Chinese), 47(6):1017-1022. |
[9] | Wu G J. 1980. Applied Geophysics-Magnetic Course (in Chinese). Beijing:Geological Publishing House. |
[10] | Xu S Z. 2006. The integral-iteration method for continuation of potential fields. Chinese J. Geophys. (in Chinese), 49(4):1176-1182. |
[11] | Xue Q F. 1978. Field Theories (in Chinese). Beijing:Geological Publishing House. |
[12] | Yan H, Xiao C H, Zhang Z Y, et al. 2010. Iterative method for continuation of three-component magnetic field. Chinese Journal of Computational Physics. (in Chinese), 27(05):705-710. |
[13] | Yang W C. 1997. Inverse theories and methods of geophysical problems (in Chinese). Beijing:Geological Publishing House. |
[14] | Zhang Q S. 1985. Field Theory (in Chinese). Beijing:Geological Publishing House. |
[15] | 陈生昌, 肖鹏飞. 2007. 位场向下延拓的波数域广义逆算法. 地球物理学报, 50(06):1816-1822. |
[16] | 梁锦文. 1989. 位场向下延拓的正则化方法. 地球物理学报, 32(5):600-608. |
[17] | 梁昆淼.1965. 数学物理方法.北京:高等教育出版社. |
[18] | 刘东甲,洪天求,贾志海等. 2009. 位场向下延拓的波数域迭代法及其收敛性. 地球物理学报, 52:(6):1599-1605,doi:10.3969/j.issn.0001-5733.2009.06.022. |
[19] | 王兴涛. 2004. 航空重力测量数据向下延拓方法比较. 地球物理学报, 47(6):1017-1022. |
[20] | 吴功建.1980. 应用地球物理学-磁法教程.北京:地质出版社. |
[21] | 徐世浙.2006. 位场延拓的积分-迭代法.地球物理学报, 49(4):1176-1182. |
[22] | 薛琴访.1978. 场论.北京:地质出版社. |
[23] | 闫辉,肖昌汉,张朝阳等. 2010. 三分量磁场延拓的递推算法. 计算物理, 27(05):705-710. |
[24] | 杨文采.1997. 地球物理反演的理论与方法.北京:地质出版社. |
[25] | 张秋生.1985. 场论.北京:地质出版社. |