地球物理学报  2017, Vol. 60 Issue (11): 4212-4220   PDF    
引用本文
张冲, 黄大年, 刘杰. 2017. 重力场向下延拓Milne法. 地球物理学报, 60(11): 4212-4220, doi: 10.6038/cjg20171109.  
Zhang C, Huang D N, Liu J. 2017. Milne method for downward continuation of gravity field. Chinese J. Geophys. (in Chinese), 60(11): 4212-4220, doi: 10.6038/cjg20171109.  
重力场向下延拓Milne法
张冲1,2, 黄大年2, 刘杰3     
1. 中国地质科学院地球深部探测中心, 北京 100037;
2. 吉林大学地球探测科学与技术学院, 长春 130026;
3. 吉林大学仪器科学与电气工程学院, 长春 130026
收稿日期 2016-12-30, 2017-09-25 收修定稿
基金项目: 中国地质调查局地质调查项目(DD20179611)，国家高技术研究发展(863)计划项目(SS2013AA063903)，国家重点研发计划课题(2016YFC0600201)联合资助
第一作者简介: 张冲, 男, 1989年生, 博士, 主要从事重磁场数据处理与解释研究.E-mail:zhangchong12@mails.jlu.edu.cn
通讯作者: 刘杰, 女, 1969年生, 副教授, 主要从事光电测试及地球物理仪器研发工作.E-mail:l_jie@jlu.edu.cn
摘要：重力场向下延拓能够突出局部和浅部的异常信息，分离叠加的异常特征.但是向下延拓通常具有过程不稳定、下延深度小、结果不准确等问题.针对向下延拓所存在的不足，本文利用重力场及其垂向一阶导数，基于辛普森（Simpson）求积公式，推导出重力场向下延拓米尔尼（Milne）公式.将本文向下延拓方法应用于模型数据，向下延拓模型结果及误差曲线表明，相对于向下延拓快速傅里叶变换（FFT）法和积分迭代法，向下延拓Milne法的深度更大，相对误差更小；相对模型值，向下延拓Milne法能够获得稳定且准确的结果.对加拿大乃查科（Nechako）盆地地区实测航空重力数据进行本文方法向下延拓验证，处理结果表明，相对于实测异常，本文方法向下延拓结果能够很好还原实测数据，并且在进一步向下延拓中反映原始异常的趋势，增强局部和细小异常信息.
关键词: Milne公式      重力场下延      垂向导数      不同高度     
Milne method for downward continuation of gravity field
ZHANG Chong1,2, HUANG Da-Nian2, LIU Jie3     
1. SinoProbe Center-China Deep Exploration Center, Chinese Academy of Geological Science, Beijing 100037, China;
2. College of Geo-Exploration Science and Technology, Jilin University, Changchun 130026, China;
3. College of Instrumentation and Electrical Engineering, Jilin University, Changchun 130026, China
Abstract: Downward continuation of gravity field can highlight the anomaly of the local and shallow part of the earth, and can separate the anomaly feature from the overlying field. But downward continuation is usually unstable, with shallow downward continuation depth and inaccurate result. To this end, this paper presents a Milne formula method for downward continuation of gravity field by using gravity field and its first-order vertical derivatives, and based on Simpson quadrature formula. The method is applied to the synthetic model. The synthetic result and the error curve show that compared with the fast Fourier transform downward continuation method and the integral iteration downward continuation method, the new method's depth of downward continuation is larger, the relative error is smaller, and the result is stable and accurate to the model. The method is applied to the downward continuation of the measured airborne gravity data of Nechako basin in Canada. The downward continuation result of this method can restore the original data accurately. And the result can reflect the tendency of the original anomaly and enhance the local and tiny anomaly.
Key words: Milne formula    Gravity downward continuation    Vertical derivative    Different heights    
1 引言

重力场向下延拓能够突出局部、浅部的地质信息(张辉等，2012)，在重力数据处理、解释和辅助导航中具有十分重要作用(高玉文等，2012).但是重力场向下延拓是不适定问题，一直为重力勘探的研究热点.向下延拓空间域插值法(Henderson，1960)的计算复杂、精度不高；常规向下延拓FFT法(Dean，1958)不仅向下延拓因子具有高频放大效应，而且傅里叶变换的离散和截断误差会引起延拓结果发生高频振荡；正则化法(Tikhonov et al., 1968)、广义逆法(陈生昌，2007)、匹配滤波法(Pawlowski，1995)、维纳滤波法(Clarke，2012)等常规向下延拓FFT法的改进方法，虽然提高了向下延拓的稳定性，但向下延拓深度不大(一般不超过5倍点距)；向下延拓积分迭代法(徐世浙，2006)能够实现无噪声重力数据稳定向下延拓，且向下延拓点距大，但其迭代次数较多，导致计算效率降低(陈生昌和肖鹏飞，2007).

针对向下延拓所存在的计算复杂、过程不稳定、向下延拓深度小、结果不准确等不足，本文利用不同高度的重力场及其垂向一阶导数，通过Simpson求积公式，建立重力场向下延拓Milne法.使用模型数据(无噪声和含噪声)对常规向下延拓FFT法、向下延拓积分迭代法和本文向下延拓Milne法进行检验，分析了随着向下延拓深度增加，这三种向下延拓结果与真实值的均方根误差的变化曲线.为了验证本文方法的实际效果，对加拿大Nechako盆地地区实际航空重力数据进行向下延拓.

2 Milne法基本原理

重力场向下延拓可定义为如下褶积型线性积分(Blakely，1995)

(1)

其中，g(ξ, η, z)已知，表示高度为z的观测面上重力场，g(x, y, z+h)未知，表示观测面以下场源以上，深度为z+h重力场，x, y, z, ξ, η表示坐标变量，h表示观测面与向下延拓面之间距离.

在空间域，公式(1)右端可视为g(ξ, η, z)与相乘的二重积分，即g(ξ, η, z)与关于变量x, y的二维褶积：

(2)

其中，*表示褶积运算.上式为空间域向下延拓公式.通过傅里叶变换，空间域褶积关系可转化为波数域乘积关系，因此，对公式(2)等号两边进行傅里叶变换，得到波数域重力场向下延拓公式：

(3)

其中，G(u, v, z)、G(u, v, z+h)分别为g(x, y, z)、g(x, y, z+h)的二维傅里叶变换，e2π(u2+v2)1/2h为 的二维傅里叶变换，(u2+v2)1/2为波数.再对公式(3)进行傅里叶逆变换，得到波数域向下延拓(FFT下延法)公式：

(4)

其中，F-1[]表示傅里叶逆变换.

在观测面以上无源半空间内，重力场具有调和函数性质，利用某一观测面重力值可计算该空间其他面重力值，实现重力场解析延拓.由于地面以下含异常源半空间的重力场不满足调和函数性质，则向下换算时，是不适定问题，解不稳定；又因为实测重力场是有限的，而延拓计算公式中积分域为无限的，需利用已知的有限的实测重力场数据对无限积分进行替代，近似求解未知区域重力场.所以有必要推导基于物理数学性质，适用更广泛的重力场向下延拓方法.重力场在物理属性上不仅是横向连续变化的，还具有垂向连续变化性质(Fedi and Florio, 2002)；在数学原理上，重力场具有水平导数和垂向导数，满足调和函数性质，即重力场微分形式存在.对重力场的垂直方向求导数，得到重力场的微分形式：

(5)

其中，gzz(x, y, z)为重力场的垂向一阶导数，f(z, gz)为已知数值或已知表达式，表示重力场gz在z处的垂向一阶导数，z在区间[a, b]内变化.假设观测面的垂向坐标为z0，z0∈[a, b]，观测面上重力场为

(6)

由等式(5)和(6)构成如下求解重力场值gz(x, y, z)的方程：

(7)

由于重力场gz(x, y, z)满足利普希茨(Lipschitz)条件，方程(7)的解gz(x, y, z)——重力场值存在且唯一，但由于重力场方程(7)的解析解，在实际问题中不实用，因此，求解方程(7)中的重力场主要采用数值解法.将方程(7)中的微分式在区间[zn, zn+1]上积分，由牛顿(Newton)—莱布尼茨(Leibniz)微积分基本公式，得到

(8)

其中，zn=z0+nh, (n=0, 1, 2, …).

等式(8)中，可将gz(x, y, zn+1)视为对gz(x, y, zn)向下延拓h深度后的重力场.所以，通过等式(8)，可以确定向下延拓重力场值gz(x, y, zn+1)，即将等式(8)视为新形式重力场向下延拓基本公式.利用等式(8)求解重力场向下延拓，需要近似计算积分项，Milne格式是微分方程数值解法中行之有效的方法之一(王新民和术洪亮，2005)，利用Simpson求积公式，近似求解等式(8)中积分项，得到重力场向下延拓Milne法.详细推导如下，利用zn-3, zn-1, zn+1三节点的Simpson求积公式，计算等式(8)中，得到

(9)

上式为计算数学中的四阶Milne格式，局部截断误差为，其中，gz(5)(ξ)表示所求区间内某点的5阶导数值.当以(9)式中gz(x, y, zn)为观测重力场gz(x, y, z0)时，即n=0时，则等式(9)化为重力场向下延拓四阶Milne格式：

(10)

其中，h为向下延拓深度，s为计算系数，gz(x, y, z1)为待求的向下延拓深为h的重力场，gz(x, y, z-3)为向上延拓3h高的重力场，gzz(x, y, z-2)为向上延拓的2h高的重力垂向梯度，gzz(x, y, z-1)为向上延拓h高的重力垂向梯度，gzz(x, y, z0)为观测面上的重力垂向梯度.

观测面上重力异常场通常已知，重力异常垂向导数可以直接测得或者通过计算得到，不同高度的重力场及其垂向梯度可采用如下稳定且收敛的向上延拓公式求解，即

(11)

其中，Gz(x, y, z0)、Gzz(x, y, z0)分别表示gz(x, y, z0)、gzz(x, y, z0)的傅里叶变换.

将公式(11)带入等式(10)，实现重力场向下延拓四阶Milne法的计算.根据求解微分方程多步法的Milne格式具有不同阶数，重力场向下延拓公式(10)也存在不同阶公式：

(12)

式中，β0、β1、…, βn是常数系数，其取值可以参考数值计算公式(Butcher, 2008; Fornberg, 1998).gz(x, y, z0+h)表示向下延拓h深的重力场，z0表示观测面的纵坐标.gz(x, y, z0-(n+1)h)表示向上延拓(n+1)h高度的重力场，gzz(x, y, z0)表示观测面上重力场垂向梯度，gzz(x, y, z0-h)和gzz(x, y, z0-nh)表示向上延拓h和nh高度的重力场垂向梯度.实际应用该公式时，应选有限项公式计算.例如，为便于说明，本文选择等式(10)重力场向下延拓四阶Milne法，即β-1=0，α0=α1=α2=0时，等式(12)的具体化.

3 模型试验

为说明本文向下延拓方法效果，在模型试验中，给出本文四阶Milne法与常规FFT法、积分迭代法的重力场向下延拓结果比较.理论模型由两个顶底界面埋深不同、大小不一的长方体组合而成.模型效果如图 1所示：长方体1(红色标记)x、y和z方向边长分别为10 m、20 m和20 m，中心点坐标为(60, 80, 23)m，其顶界面距观测面为13 m，底界面距观测面33 m；长方体2(蓝色标记)x、y和z方向边长分别为40 m、10 m和20 m，中心点坐标为(60, 50, 25)m，其顶界面距观测面为15 m，底界面距观测面35 m.两个长方体密度相同，与围岩密度差均为0.9 kg·m-3.观测面上点、线数均为130，网格间距为1 m.模型试验计算中，不同高度的重力垂向一阶导数是将观测面上理论值向上延拓至不同高度得到的.延拓前需要进行扩边处理，防止延拓结果出现边界效应.但是扩边方法种类繁多，对于不同数据类型，不同扩边方法对延拓结果影响较大，引入扩边方法不利于比较延拓方法本身的好坏(王万银，2009; 骆遥, 2016).因此，为显明延拓方法本身的优缺点，本文向下(向上)延拓所涉及三种方法均未进行扩边、去噪滤波以及正则化处理.

图 1 地下模型空间分布 Fig. 1 The position of underground rectangular model

不含噪声时，FFT法(图 2e)、积分迭代法(图 2g：迭代步长为1，迭代收敛误差为0.0001 mGal，模型数据的精度也为0.0001 mGal，迭代次数为4次)和本文方法(图 2j)均能对重力异常向下延拓2倍点距.FFT法的向下延拓结果具有强烈边界损失，积分迭代法和本文方法的向下延拓边界损失小，向下延拓结果均稳定且收敛.同时对比向下延拓2倍点距的理论异常值(图 2c)，积分迭代法和本文方法的向下延拓结果都十分准确，说明了本文方法进行向下延拓的有效性.含噪声时，FFT法(图 2f)和积分迭代法(图 2h：迭代步长为1，迭代收敛误差为0.001 mGal，迭代次数为50次，建议最大次数.)结果发散，不能获得向下延拓2倍点距结果，对比含噪声向下延拓FFT法和积分迭代法结果，本文方法(图 2i)对噪声不敏感，向下延拓相对稳定，得到了相对明显的向下延拓结果.无论是否存在噪声干扰，本文方法均可获得较为稳定且准确的向下延拓结果.

图 2 向下延拓2倍点距(2 m)结果(图中色标单位为mGal) (a)不含噪声观测面理论计算重力异常值；(b)含5%噪声观测面理论计算重力异常值；(c)不含噪声向下理论计算重力异常值；(d)含噪声向下理论计算重力异常值；(e)利用FFT法对(a)向下延拓结果；(f)利用FFT法对(b)向下延拓结果；(g)利用积分迭代法对(a)向下延拓结果；(h)利用积分迭代法对(b)向下延拓结果；(i)利用本文方法对(a)向下延拓结果；(j)利用本文方法对(b)向下延拓结果. Fig. 2 The synthetic result of downward continuation with 2 m (a) The theoretical gravity anomalies over the ground without noises; (b) The theoretical gravity anomalies over the ground with 5% noises; (c) The theoretical gravity anomalies of 2m; (d) The theoretical gravity anomalies of 2m with 5% noises; (e) The downward continuation anomalies by FFT method with the results of (a); (f) The downward continuation anomalies of by FFT method with the results of (b); (g) The downward continuation anomalies by Integral-iteration method with the results of (a); (h) The downward continuation anomalies by Integral-iteration method with the results of (b); (i) The downward continuation anomalies by this method with the results of (a); (j) The downward continuation anomalies by this method with the results of (b).

为进一步阐明本文方法的优越性，向下延拓深度增大至10倍点距.在无噪声情况下，FFT法(图 3e)向下延拓结果发散，不能获得向下延拓结果，积分迭代法(图 3g：迭代步长为1，迭代收敛误差为0.0001 mGal，迭代次数为50次.)和本文方法(图 3j)的结果均稳定且收敛，可以延拓至10倍点距乃至更大点距，但相对于积分迭代法结果，本文方法边界损失更小、计算结果更接近理论异常(图 3c).在噪声干扰下，FFT法(图 3f)和积分迭代法(图 3h：迭代步长为1，迭代收敛误差为0.001 mGal，迭代次数为50次.)均发散，不能得到向下延拓结果，而本文方法(图 3i)仍有很好效果，在向下延拓深度大时，本文方法对噪声依然不是特别敏感，能够获得不错延拓结果.

图 3 向下延拓10倍点距(10 m)结果(图中色标单位为mGal) (a)不含噪声观测面理论计算重力异常值；(b)含5%噪声观测面理论计算重力异常值；(c)不含噪声向下理论计算重力异常值；(d)含噪声向下理论计算重力异常值；(e)利用FFT法对(a)向下延拓结果；(f)利用FFT法对(b)向下延拓结果；(g)利用积分迭代法对(a)向下延拓结果；(h)利用积分迭代法对(b)向下延拓结果；(i)利用本文方法对(a)向下延拓结果；(j)利用本文方法对(b)向下延拓结果. Fig. 3 The synthetic result of downward continuation with 10 m (a) The theoretical gravity anomalies over the ground without noises; (b) The theoretical gravity anomalies over the ground with 5% noises; (c) The theoretical gravity anomalies of 10m; (d) The theoretical gravity anomalies of 10m with 5% noises; (e) The downward continuation anomalies by FFT method with the results of (a); (f) The downward continuation anomalies by FFT method with the results of (b); (g) The downward continuation anomalies by Integral-iteration method with the results of (a); (h) The downward continuation anomalies by Integral-iteration method with the results of (b); (i) The downward continuation anomalies by this method with the results of (a); (j) The downward continuation anomalies by this method with the results of (b).

图 4a是无噪声数据FFT法和本文方法的结果与真实值的均方误差随下延点距的对数函数变化曲线，图 4b是无噪声数据积分迭代法和本文方法的结果与真实值的均方误差随向下延拓点距的变化曲线.通过对比该误差曲线可知，向下延拓相同点距时，FFT法的均方误差随向下延拓深度增加指数增大，而本文方法的均方误差小很多，同时，当向下延拓深度增大到5倍点距之后，本文方法的均方误差小于积分迭代法.均方误差曲线说明本文方法的结果更准确，更接近真实异常值.

图 4 误差曲线 (a) FFT法(红色线)和本文方法(绿色线)向下延拓结果与真实结果的均方误差对数函数随向下延拓点距的变化曲线; (b)积分迭代法(蓝色线)和本文方法(绿色线)向下延拓结果与真实结果的均方误差随向下延拓点距的变化曲线. Fig. 4 The RMS error curve (a) The logarithmic function of RMS variation with the downward continuation depth of FFT and this method; (b) The RMS variation with the downward continuation depth of integral-iteration and this method.
4 实际应用

为了检验所提出的向下延拓Milne法在实际数据处理的实用性，本文对加拿大Nechako盆地北部地区实测航空重力数据进行向下延拓试验.加拿大Nechako盆地地区位于加拿大西南部的大不列颠哥伦比亚省的中部，是该省八个盆地中最大的一个，覆盖面积约为69000 km2.地表以第三纪始新世早期到第三纪晚期的火山熔岩为主，盆内富产石油和天然气(Calvert et al., 2011；Kim et al., 2014).

该重力观测资料的原始区域的平均经度约为-126.22°，最小经度约为-95.65°，最大经度约为-151.33°.飞机飞行高度为2000 m.所选用的重力数据仅为该测区观测资料的一部分，其网格间距为400 m×400 m，网格点数为179×179，航空布格重力异常以及垂向一阶导数均为实测资料.实测航空布格重力异常结果为图 5a所示，将其向上延拓2000 m(5倍点距)，作为已知重力异常如图 5b.将2000 m高度重力异常向上延拓8000 m，再将2000 m高度的重力异常垂向一阶导数分别向上延拓至4000 m和6000 m.利用本文向下延拓Milne法(公式(10))，对图 5b向下延拓2000 m(5倍点距)，得到结果图 5c.延拓后结果(图 5c)与原始异常(图 5a)的误差结果如图 5d.通过向下延拓5倍点距的结果(图 5d)可知，本文向下延拓Milne法稳定且收敛.在进行向下延拓的过程中未使用任何扩边方法，从误差图像可以看出本文向下延拓Milne法的边界效应不明显，对比图 5a和图 5c可以看出本文延拓结果与实测重力异常相吻合，原本在更大高度位置图 5b中已经不可见的异常，在图 5c中均有很好、很准确的体现.

图 5 加拿大Nechako盆地重力数据向下延拓比较结果 (a)该地区布格重力异常结果；(b)是(a)向上延拓5倍点距结果；(c)是由(b)经过本文方法向下延拓5倍点距结果；(d)是(c)与(a)的误差. Fig. 5 The comparison result of downward continuation over Nechako in Canada (a) The gravity Bouguer anomaly over this area; (b) The result of upward continuation of 5 times the interval; (c) The downward continuation of (b)′s result of 5 times the interval; (d) The error between (c) and (a).

我们再将图 5a的实测布格重力异常向下延拓至2000 m(5倍点距)、4000 m(10倍点距)，图 6a和图 6b分别给出了相应向下延拓结果.可以看出，向下延拓2000 m(5倍点距)、4000 m(10倍点距)后得到的重力异常没有出现强烈的边界损失，没有发生剧烈震荡，增强了局部异常的同时基本保持原始异常变化趋势.

图 6 加拿大Nechako盆地重力数据实际向下延拓结果 (a)该地区布格重力异常经过本文方法向下延拓5倍点距结果；(b)该地区布格重力异常经过本文方法向下延拓10倍点距结果. Fig. 6 The result of downward continuation over Nechako in Canada (a) The downward continuation result of 5 times the interval of the gravity Bouguer anomaly over this area; (b) The downward continuation result of 10 times the interval of the gravity Bouguer anomaly over this area.
5 结论

本文通过引入不同高度的垂向一阶导数，利用Milne公式，给出重力场向下延拓Milne法.由于垂向一阶导数容易计算或实测，其向上延拓结果稳定收敛，且导数的引入会增加细节信息，所以使向下延拓结果更准确.利用模型数据，将本文方法与FFT法和积分迭代法进行对比验证，均未使用任何扩边方法，FFT法结果出现强烈的边界损失，而本文方法结果边界效应不明显.FFT法向下延拓2倍点距已经发散，而本文方法可以向下延拓至10倍点距甚至异常体上顶面，表明本文方法向下延拓深度大，且有很好的稳定性和收敛性.对本文方法与积分迭代法做比较研究，虽然对不含噪声数据，积分迭代法结果与本文方法结果相当，但含噪声时，积分迭代法结果发散，本文方法结果稳定，说明本文方法具有很好的抗噪能力，同时相对于需要反复迭代和对噪声具有累加作用的积分迭代法，本文方法向下延拓的计算更直接、更稳定且结果更准确.本文方法所得到向下延拓结果对比真实重力异常值，表明本文方法，能够很好刻画并反映真实异常值.分析FFT法、积分迭代法和本文方法的向下延拓重力场值与真实值的均方误差随着向下延拓深度的变化关系.本文方法的均方误差远远小于FFT法，当向下延拓深度大于5倍点距之后，本文方法的均方误差小于积分迭代法，且随着迭代深度的增加，误差曲线变大的不明显，表明本文方法具有很好的精确度.为进一步阐释本文方法实用效果，将本文方法应用到实际数据向下延拓，得到了稳定准确的向下延拓结果，边界效应不明显.本文方法计算向下延拓场与实际数据之间误差小，与原始异常变化趋势相吻合，很好的恢复刻画了向下延拓场中的微小异常，表明了本文方法具有稳定、可靠且准确实际应用性.

致谢

感谢加拿大自然资源部(Natural Resources Canada)提供乃查科盆地(Nechako Basin)实测航空重力数据.感谢两位匿名审稿专家的意见以及地球物理学报主编的建议和编辑的修订.

参考文献
Blakely R J. 1995. Potential Theory in Gravity and Magnetic Applications. Cambridge: Cambridge University Press.
Butcher J C. 2008. Numerical Methods for Ordinary Differential Equations. Chichester: John Wiley & Sons.
Calvert A J, Hayward N E, Spratt J E, et al. 2011. Seismic reflection constraints on upper crustal structures in the volcanic-covered central Nechako basin, British Columbia. Canadian Journal of Earth Sciences, 48(6): 1021-1037. DOI:10.1139/e10-097
Chen S C, Xiao P F. 2007. Wavenumber domain generalized inverse algorithm for potential field downward continuation. Chinese Journal of Geophysics, 50(6): 1816-1822.
Clarke G K C. 2012. Optimum second-derivative and downward-continuation filters. Geophysics, 34(3): 424-437.
Dean W C. 1958. Frequency analysis for gravity and magnetic interpretation. Geophysics, 23(1): 97-127. DOI:10.1190/1.1438457
Fedi M, Florio G. 2002. A stable downward continuation by using the ISVD method. Geophysical Journal International, 151(1): 146-156. DOI:10.1046/j.1365-246X.2002.01767.x
Fornberg B. 1998. Calculation of weights in finite difference formulas. SIAM Review, 40(3): 685-691. DOI:10.1137/S0036144596322507
Gao Y W, Luo Y, Wen W. 2012. The compensation method for downward continuation of potential field from horizontal plane and its application. Chinese Journal of Geophysics, 55(8): 2747-2756. DOI:10.6038/j.issn.0001-5733.2012.08.026
Henderson R G. 1960. A comprehensive system of automatic computation in magnetic and gravity interpretation. Geophysics, 25(3): 569-585. DOI:10.1190/1.1438736
Kim H S, Cassidy J F, Dosso S E, et al. 2014. Mapping crustal structure of the Nechako-Chilcotin plateau using teleseismic receiver function analysis. Canadian Journal of Earth Sciences, 51(4): 407-417. DOI:10.1139/cjes-2013-0147
Luo Y, Wu M P. 2016. Minimum curvature method for downward continuation of potential field data. Chinese Journal of Geophysics, 59(1): 240-251. DOI:10.6038/cjg20160120
Pawlowski R S. 1995. Preferential continuation for potential-field anomaly enhancement. Geophysics, 60(2): 390-398. DOI:10.1190/1.1443775
Tikhonov A N, Glasko V B, Litvinenko O K, et al. 1968. Analytic continuation of a potential in the direction of disturbing masses by the regularization method. Izv., Earth Physics (in Russian), 12: 30-48.
Wang W Y, Qiu Z Y, Liu J L, et al. 2009. The research to the extending edge and interpolation based on the minimum curvature method in potential field data processing. Progress in Geophysics, 24(4): 1327-1338. DOI:10.3969/j.issn.1004-2903.2009.04.022
Wang X M, Shu H L. 2005. Engineering Mathematics: Numerical Methods. Beijing: Higher Education Press.
Xu S Z. 2006. The integral-iteration method for continuation of potential fields. Chinese Journal of Geophysics, 49(4): 1176-1182.
Zhang H, Chen L W, Ren Z X, et al. 2009. Analysis on convergence of iteration method for potential fields downward continuation and research on robust downward continuation method. Chinese Journal of Geophysics, 52(4): 1107-1113.
陈生昌, 肖鹏飞. 2007. 位场向下延拓的波数域广义逆算法. 地球物理学报, 50(6): 1816–1822.
高玉文, 骆遥, 文武. 2012. 补偿向下延拓方法研究及应用. 地球物理学报, 55(8): 2747–2756. DOI:10.6038/j.issn.0001-5733.2012.08.026
骆遥, 吴美平. 2016. 位场向下延拓的最小曲率方法. 地球物理学报, 59(1): 240–251. DOI:10.6038/cjg20160120
王万银, 邱之云, 刘金兰, 等. 2009. 位场数据处理中的最小曲率扩边和补空方法研究. 地球物理学进展, 24(4): 1327–1338. DOI:10.3969/j.issn.1004-2903.2009.04.022
王新民, 术洪亮. 2005. 工程数学:计算方法. 北京: 高等教育出版社.
徐世浙. 2006. 位场延拓的积分-迭代法. 地球物理学报, 49(4): 1176–1182.
张辉, 陈龙伟, 任治新, 等. 2009. 位场向下延拓迭代法收敛性分析及稳健向下延拓方法研究. 地球物理学报, 52(4): 1107–1113.