2012, Vol. 55
2. 中国国土资源航空物探遥感中心,北京 100083
2. China Aero Geophysical Survey and Remote Sensing Center for Land and Resources, Beijing 100083, China
位场数据处理与解释中，向下延拓是多种数据处理以及反演解释中需使用的重要技术手段.例如，在重力界面反演中需要进行多次的向下延拓；此外，向下延拓能够突出局部异常，突出浅部的地质信息.但是，位场的向下延拓同向上延拓相反，是地球物理中典型的不适定问题，频率域的向下延拓因子具有极强的高频放大作用；此外，傅里叶变换的离散及其截断误差，也会引起延拓结果的高频振荡.为了克服向下延拓的不稳定性，通常针对向下延拓因子的高频放大特性对向下延拓算子加低通类滤波的方法保证向下延拓的稳定性[1].众多学者提出各种改善延拓稳定性的方法，如正则化方法[2-3]，Wiener滤波方法[4]，多尺度边缘约束法[5]，匹配组合滤波方法[6-7]，补偿与逐次向下延拓方法[8]，以及泰勒级数展开方法[9]和广义逆方法[10-11]等.这些改进方法使向下延拓稳定性得到一定的提高，但随着向下延拓的跨度增大，向下延拓深度仍不能超过点距的3～5 倍.最近，徐世浙[12-14]系统研究了位场向下延拓迭代法，在无噪声的模型研究中将延拓的距离提高到20 倍点距实现了极大跨度的向下延拓.
尽管，徐世浙空间域迭代法能较好地解决位场向下延拓不稳定性这一问题，实现大跨度的向下延拓，但陈生昌等[11]认为由于需要迭代的次数较多，每次迭代均需要进行位场的向上延拓致使计算效率较低.陈生昌给出了一种频率域下延因子，其延拓结果很大程度上取决于阻尼因子，当阻尼因子适当放大时延拓精度会迅速下降.本文借鉴空间域向下延拓迭代法的思想[15]，将位场向下延拓的波数域广义逆算法[11]与迭代法[13]结合，形成一种新的向下延拓补偿方法，通过逐次补偿的办法实现了位场稳定的向下延拓，同时将这种补偿因子引入三维重力密度界面反演中，取得了较好的效果.
2 补偿原理及方法利用傅里叶变换矩阵的正交对称特性及广义逆原理，陈生昌等将向下延拓的褶积型线性积分方程变换为主对角阵形式的下延奇异值，并向奇异值中引入阻尼，给出了波数域傅里叶变换广义逆向下延拓算子[11]:
|
(1) |
其中，h（h＞0）是向下延拓的距离；u、v分别为水平方向上的波数；α（α≥0）是阻尼因子.当α=0 时，上述广义逆向下延拓算子将退化为普通的频率域下延因子.实际上，陈生昌等推导的结果可以认为是Li和Oldeuburg[16]频率域位场反演中不考虑模型偏导数的一种特解.
通过（1）式广义逆向下延拓算子可以较为稳定地实现位场的向下延拓，但实际计算中向下延拓的距离与延拓精度却是相互矛盾的.当阻尼因子的取值过大时，将会压低位场延拓数据的精度（幅值）、降低延拓数据的精度；当的取值过小或时，对于高波数数据或大跨度的延拓会产生不稳定性.为了实现大跨度的向下延拓并保持延拓结果的稳定，我们引入一个较大的阻尼因子，由此带来的延拓精度较低的问题，则采取补偿的方法进行弥补.
在补偿过程中，将（1）式计算出的下延结果向上延拓至原观测平面，用观测面上的实测值与计算值的差值进行向下延拓，对下延平面上的位场场值进行补偿.如此反复进行补偿，直至观测面上的实测值与向上延拓计算值的差值可以忽略.上述补偿过程类似于迭代法，其具体补偿形式为
|
(2) |
其中gB是已知观测平面B上的位；gAn是待求的向下延拓平面A上的位，n是补偿次数；D（g）是对位场g进行向下延拓；U（g）是对位场g进行向上延拓.实际补偿延拓中，向上延拓采用普通的频率域延拓因子，向下延拓通过式（1）广义逆向下延拓算子进行频率域下延.
与迭代法不同，补偿法考虑了迭代初始过程中初始位场gA0更为精确的给出，实际上，如果令将向下延拓进一步近似为D（gB）=gB，则有：
|
(3) |
即演变为向下延拓的迭代算法，可见迭代法可以认为是补偿法的一种特殊情况.由于补偿过程中考虑了迭代的初值，并通过向下延拓来进行较为准确的补偿，所以在适当的阻尼条件下，该补偿方法将具有更快的收敛速度，能够通过较少的补偿次数获得较高的延拓精度.
3 补偿法的频率域响应因子 3.1 频率响应因子尽管上述补偿方法可以通过较少的补偿次数实现精度较高的向下延拓，但由于每次补偿过程需要进行两次频率域延拓，计算效率较低，为此我们从频率域空间对上述补偿方法进行讨论.
将补偿公式（2）两端进行傅里叶变换，有
|
(4) |
其中，HD 是式中定义的广义逆向下延拓算子；${H_{\rm{U}}}{\rm{ = }}{{\rm{e}}^{ - 2\pi h\sqrt {{u^2} + {v^2}} }}$是频率域上延因子；GB是已知观测平面B上位gB的频谱；GAn是待求的向下延拓平面A上位gAn的频谱；n是补偿次数.
设第n次补偿的频率域响应因子为Hn，则有
|
(5) |
代入（4）式，利用频率域中向上延拓因子和广义逆向下延拓算子的虚部为零的特性，可以得到：
|
(6) |
再将具体的解析延拓因子代入（6）式即可以得到具体的补偿法向下延拓的频率域响应因子：
|
(7) |
（7）式中物理量的具体含义与（1）式广义逆向下延拓算子完全相同，不再重复.
上面得到的下延频率域补偿响应因子，是通过逐次补偿来获取大跨度向下延拓的精确频率响应，其原理与侯重初[17]提出的补偿圆滑滤波方法在补偿思想以及具体给出形式上都是非常一致的，实际上侯重初后来提出的经过补偿的向下延拓滤波算子[8]是从频率域角度进行了一次补偿处理，在不考虑侯重初使用的向下延拓的滤波算子与本文采用的广义逆向下延拓算子具体形式时，选取（6）式进行一次补偿，就可得出侯重初所给出向下延拓滤波算子[8]是本文频率域响应的一阶近似这一结论，而本文则根据补偿原理进行了更多次的频率域补偿，因此这种频率域补偿响应因子是对侯重初向下延拓滤波算子的一种发展.
3.2 频率特性补偿法作为一种实用化的方法技术，从理论上给予认识并揭示其数学物理本质，往往是有意义的.为此本文尝试从补偿法或迭代法本身及其傅里叶空间的频率特性两个角度来进行讨论.式（2）的补偿公式实际上是一种广义反演，观测位场gB可认为是数据空间，下延位场gA则可认为是模型空间，于是
|
(8) |
其中，U是向上延拓运算，是从模型空间到数据空间的映射关系.于是，观测位场可以通过向上延拓这种“正演”方式，用下延位场进行计算，而求解下延位场就是通过观测位场来进行“反演”求解，而这种反演则是不适定的.式（2）的补偿法首先通过一种近似的稳定向下延拓给出初始模型，反演过程则是通过对模型的不断进行校正，可视为试错法反演，而每一次模型的校正值是根据式（8）正演拟合得到的.从反演过程来看，补偿法与迭代法的区别仅在于初始模型的给出及校正值的计算.实际计算表明经过一定次数的迭代后，正演拟合的位场（上延位场）总能趋近于观测位场，而对下延位场的校正值则由拟合误差决定，经多次迭代计算出的下延位场确实趋近于某一位场，迭代法是事实收敛的.当选用补偿法进行下延时，由于阻尼的作用式（1）进行下延校正的过程中能一定程度地压制高频成分---甚至较上延操作更能整体地压制位场的频谱，所以补偿法也是事实收敛的.同其他的反演方法一样，上述反演的收敛性是可以很容易的理解的，但其严格的证明却是比较困难的.
对于从频率域空间讨论迭代法收敛性的问题，张辉等[18]认为其收敛于FFT 法下延的理论解，但FFT 法下延出的位场往往发散[13]，而迭代法计算表明计算结果能够逼近于某一稳定的位场，实际上空间域的迭代还隐含了低通滤波的操作.我们认为避开空间域收敛性的问题，直接讨论其频率域的滤波特性更具有实际意义.为此，进一步考虑普通的下延频率因子H、张辉等给出的下延因子Rn、陈生昌等给出的下延因子HD 即式（1）以及本文补偿法的下延因子Hn即式（7）的频率响应特征，图 1给出了各种下延因子的频率响应情况（其中Rn、HD 及Hn的频谱放大程度被调节的近似一致）.可以看出普通的下延频率因子对频谱指数型的放大，而其他改进的因子能对高频成分进行一定的压制，显然Rn对高频部分压制效果不如HD，HD 与本文的响应因子则非常类似，但本文的响应因子在对某一低频段逼近H的同时，较HD 更能压制高频成分，即低频段能更精确地逼近普通下延因子而对高频成分能更为有效地进行压制因而也更稳定.当然，由于FFT 给出的频谱是复数空间，很难对比观测位场和下延位场的频谱差异，这需要借助其它数学物理方法进行更深入的分析，给出迭代法的真实频率响应特征.
|
图 1 各种下延因子的频率响应特征 Fig. 1 The characteristic of different downward continuation factor in frequency domain |
考虑从下延频率域补偿响应因子中提取出了隐含的低通滤波因子H* ，令
|
(9) |
于是
|
(10) |
提取的低通滤波因子H* 其表达式非常复杂，其构造的低通滤波器涉及下延距离h、阻尼因子α及补偿次数n，用频率响应曲线更能揭示其滤波实质（以下讨论中向下延拓的数据其点距均归一化为1，既Δx=1m）.如图 2所示，随着补偿次数n的增加，得到一组低通滤波曲线族，其中α=0.1、h=5Δx时，曲线H* 沿波数向前移动，每条曲线都具有一段全通区域后快速下降并趋近于零，补偿次数越多，通过的低通滤波因子的频率成分越多，这与侯重初构建的补偿圆滑滤波相类似.其次，考虑下延距离h的影响，如图 3 所示，n=200、α=0.5 时，随着下延距离h的加大，曲线H* 沿波数向低频区域移动，且曲线的陡度在增加，即随着下延深度的增大，通过的低通滤波因子的频率成分越少.最后，将补偿次数固定，取n=100、h=5Δx时，如图 4所示，随着阻尼因子的增大，曲线H* 沿波数向低频区域移动，即α 越大，通过的低通滤波因子的频率成分越少.不难看出补偿法所构成的低通滤波器H* 具有理想低通滤波器的特性，通过对下延距离h、阻尼因子α 及补偿次数n的调节，能构造出满足不同频率滤波特性的低通滤波器，可以说H* 是一种很理想的补偿圆滑滤波因子，一定程度上解释了补偿法的收敛问题，这是补偿圆滑滤波理论运用的又一发展.
|
图 2 不同补偿次数下的频率响应特征图中自左向右的曲线其补偿次数依次为1、10、50、100、200、500、1000、2000 Fig. 2 The characteristic of the low-pass filter with different compensation in frequency domain. The compensation numbers from left to right are, in order: 1、10、50、100、200、500、1000、2000 |
|
图 3 不同下延深度的频率响应特征 图中自左向右的曲线其下延深度依次为20、15、10、5、4、3、2、1 Fig. 3 The characteristic of the low-pass filter with different depths of downward continuation in frequency domain. The depths of downward continuation from left to right are， in order： 20、15、10、5、4、3、2、1 |
|
图 4 不同阻尼因子条件下的频率响应特征图中自左向右的曲线其阻尼因子依次为10、3、1、0.3、0.1、0.03、0.01、0.003 Fig. 4 The characteristic of the low-pass filter with different damping factor in frequency domain. The damping factors from left to right are, in order： 10、3、1、0.3、0. 1、0. 03、0. 01、0. 003 |
为了检验位场向下延拓的补偿方法及其频率响应的准确性和实用效果，本文采用三维理论模型数据进行数值试验.我们选取了郭志宏等[19]提出的磁性长方体无解析奇点模型作为标准参照，具体的模型设置和网格参数与其完全一致.正演计算距离地面2km 高空处的ΔT场作为延拓数据，下延2km（20倍点线距）延拓出地面处的ΔT场，以地面处的ΔT场正演计算结果作为评价标准，进行检验.图 5a和图 5b分别给出地面处和2km 高空处ΔT场的正演结果.
|
图 5 磁性长方体AT场理论模型 (a)是地面处场；（b)是距地面2 km高空处场；网格间距100 m. Fig. 5 Total field anomalies of cuboid model (a) The total field data at the height of 0 km； (b) The total field data at the height of 2 km. The grid interval is 100 m, respectively. |
采用空间域和频率域的补偿法对图 5b给出的距地面2km 高空处场向下延拓至0 m 高度，其中阻尼因子分别取0.01 和0.05.当阻尼因子选取较小时，图 6a未经补偿的延拓结果与图 5a的理论模型结果相比无论从幅度还是形态仍有一定差距，这主要因为在大跨度的下延中，阻尼为了保持稳定下延会压低位场延拓数据的幅度，降低精度.而图 6b给出的经18次补偿的向下延拓结果与图 5a的理论模型相比较无论从形态上还是从具体的数值上都是非常一致的，二者的均方误差不超过0.05nT.当响应因子中的阻尼取0.05时，图 6c给出未经补偿的向下延拓结果，图 6d给出了经38 次补偿的频率域向下延拓结果.对比图 6a和图 6c未经补偿直接采用广义逆向下延拓算子的延拓结果可以发现，当阻尼因子加大后向下延拓结果将更加稳定，例如图 6a中-20nT 的等值线略有畸变；但加大阻尼因子的同时会降低延拓结果的精度，表现为图 6c异常的幅度明显低于图 6a，与图 5a理论模型结果最大幅度相差竟高达110nT 以上.当向下延拓跨度过大时，为了保持延拓的稳定，适当放大阻尼因子会大大降低广义逆向下延拓算子向下延拓的精度.但是图 6d通过频率域补偿获得的大跨度向下延拓结果则明显增强了延拓结果的精度，异常的幅度和形态都有较大的恢复，可以看出除底部因边界造成0nT 等值线略有不同外，下延结果非常稳定且具有较高的精度.可见给出的频率域补偿响应因子在解决向下延拓跨度和精度之间具有调和作用，能够用于位场资料向下延拓.
|
图 6 向下延拓的效果检验 (a) = 0.01，未经补偿的向下延拓结果；（b) = 0.01，空间域18次补偿的向下延拓结果；(c) = 0.05,未经补偿的向下延拓结果；（d) = 0.05,频率域38次补偿的向下延拓结果. Fig. 6 The examining results of downward continuation (a) α = 0. 01，Total field anomalies by downward continuation without compensation； (b) α=0. 01，with 18th continued compensation in spacing domain； (c) α = 0. 05，Total field anomalies by downward continuation without compensation； (d)α=0. 05，with 38th continued compensation in frequency domain. |
为了检验向下延拓的补偿方法及其频率响应对实际资料的处理效果，对南非某地的实测重力异常数据进行延拓处理.该重力资料处于-26.1°E--26.9°E、-26.9°N--26.1°N 范围，网格间距0.004°，如图 7a所示.将该网格数据利用频率域位场转换向上延拓20倍网格距（约9km），得到图 7b的重力异常资料.可以看出由于延拓高度极大，上延后的重力资料异常极其平缓，将图 7b作为原始资料用位场向下延拓的补偿方法及其频率响应分别进行大跨度的向下延拓，以检验对实际资料的处理效果.
|
图 7 南非某地实际重力资料（a)和向上延拓20倍网格距的上延结果（b) Fig. 7 Contour map of gravity data in South Africa (a) and result of upward continuation 20 times interval of grid (b) |
分别采用位场向下延拓的补偿方法及其频率响应对图 7b进行20倍网格距的向下延拓，图 8a和图 8b分别给出空间域和频率域补偿的向下延拓结果.可以看出，二者延拓获得的位场与图 7a非常吻合，仅仅是细微的局部异常无法体现，位场的形态和幅度都获得了极大程度的恢复，可见用图 7b极其平缓的上延位场资料能较好地恢复出原始位场资料，说明位场向下延拓的补偿方法及其频率响应具有较高的计算精度，能够用于实际位场资料的稳定向下延拓.当然，图 7b上延位场过程中无疑会损失图 7a中原始资料中的部分高频异常的信息，这部分信息是无法通过任何数学方法所弥补的，向下延拓是无法延拓出这部分异常信息的，这一点在应用补偿法向下延拓时要尤为注意.
|
图 8 空间域补偿法的下延结果（a)和频率域响应因子下延的结果（b) Fig. 8 Gravity data by downward continuation with compensation in spacing domain (a) and result by downward continuation with compensation in frequency domain (b) |
重力异常中含有密度分布和密度分界面的深度及起伏的重要信息，根据重力资料进行密度界面的反演是位场处理解释的重要内容之一[20].计算地壳底界面即莫霍面及密度界面分层，对于研究区域构造及探测矿产资源具有重要的意义.与反射地震方法相比，重力勘探研究莫霍面的起伏和深度变化具有成本低、适用范围广的优势.国内外众多学者提出了许多密度界面反演的方法，如，Bott[21]提出了根据剩余重力异常计算二维沉积盆地深度的方法，刘元龙等[22]提出用压缩质面法反演地壳构造，Oldenburg[23]根据Parker[24]公式提出了频率域的密度界面的迭代反演方法.这些方法极大地丰富了重力资料反演方法理论，其中频率域方法以其较快的运算速度而被广泛应用，但该方法中包含向下延拓算子.为了限制高频振荡，通常在反演过程中加入低通滤波，此处将本文的下延频率域补偿响应因子引入界面反演中，改进算法的稳定性.
参考Bott和Oldenburg方法，给出一种三维重力界面反演的迭代形式：
|
(11) |
式中，hi（x，y）是第i次迭代时（x，y）处的密度界面的平均深度z0 的起伏情况；u、v分别为水平方向上的波数，Δg是由密度界面起伏引起的重力异常；Δgi表示起伏界面hi（x，y）产生的重力异常；G为万有引力常数；ρ 为界面密度差；F-1[]表示傅里叶正变换；⊗表示褶积.经过逐次迭代逼近，使计算场与实测场之间的差越来越小，但由于下延因子
为了检验实际效果，利用法国某地区的重力资料[25-26]进行莫霍面起伏的反演研究，如图 9所示，该区域范围是300km×300km，重力异常幅度约为40mGal.反演过程中，莫霍面上下层的密度差选用0.4g/cm3，平均深度为30km.
|
图 9 法国布列尼塔地区莫霍面重力异常(等值线间距2 mGal) Fig. 9 Bouguer anomaly map attributed to Moho interface inBtittany (France) The contour interval s 2 mGal，and data from Gomez-Ortiz and Agarwal (2005) |
图 10a是由法国布列尼塔地区莫霍面重力异常反演得到的基地起伏情况，莫霍面最大深度为32.3km，最小深度为27.6km.反演过程中仅通过4次迭代，便使得计算场与实际场之间的均方根误差降为0.75mGal.与以往反演结果[26]相比较，莫霍面起伏形态基本一致，仅是在幅值有些偏差，最大差值为0.2km，75%以上的数值差异小于0.08km（图 10b）.
|
图 10 重力反演的法国布列尼塔地区莫霍面深度a和与以往反演结果的差值(b) Fig. 10 Moho depth map of inversion in Brittany (France) (a) and difference of Moho depth between Fig. 10a and the inversion result before (b) |
为了进一步验证反演方法的效果，将得到的莫霍面深度结果进行Parker算法正演，如图 11所示.图 11a是由基底界面起伏得到的重力异常，图 11b为正演得到的重力异常与布列尼塔地区莫霍面重力异常的差值.可以看到，图 11a与图 4 之间非常吻合，无论从形态上还是从具体的数值上都是十分一致，误差范围在-1.45～ -0.04mGal.通过引入稳定性较好的向下延拓因子，一定程度上解决了下延计算不稳定问题，达到了较理想的效果.
|
图 11 莫霍面深度起伏引起的重力异常（)和与实际重力异常的差值(b) Fig. 11 Gravity anomaly map computed from Moho relief (a) and difference between anomaly of Fig. 11a and original anomaly presented in Fig. 9 (b) |
基于空间域向下延拓的迭代法以及补偿圆滑滤波的思想，将位场向下延拓的波数域广义逆算法与迭代法相互结合，通过逐次补偿的办法实现了位场数据稳定的向下延拓.同时，在频率域空间讨论了所提出的位场向下延拓的补偿方法，并给出了这种延拓方法的频率域响应因子.对比理论模型和实际位场资料的下延结果表明，位场向下延拓的补偿方法及其频率响应在一定程度上能够实现位场资料的稳定下延，且具有较高的延拓精度.通过将该补偿方法引入重力密度界面反演中，取得了较为理想的效果.但值得注意的是补偿法的本质仍是迭代法，在位场资料向下延拓处理时存在一定的应用条件，只有可下延位场才具备下延条件[27]，最近姚长利等[28]对重磁位场转换计算中迭代法的综合分析仍适用于对补偿法向下延拓的分析.
致谢感谢南非Witwatersrand 大学地学系Gordon R.J. Cooper教授提供了南非重力资料，感谢David Gomez Ortiz博士提供的法国重力资料.
| [1] | 北京师范大学数学系, 国家地质总局一五〇工程方法组. 用二维富氏变换进行重磁位场的转换. 北京师范大学学报 , 1978(2): 56–69. Mathematics Department of Beijing Normal University, 150 Engineering Group of General Administration of geology of the People's Republic of China. Potential field transform of gravity and magnetic data by Fourier transform. Journal of Beijing Normal University (in Chinese) (in Chinese) , 1978(2): 56-69. |
| [2] | 栾文贵. 场位解析延拓的稳定化算法. 地球物理学报 , 1985, 28(5): 537–543. Luan W G. The stabilized algorithm of the analytic continuation for the potential field. Chinese J. Geophys (in Chinese) , 1985, 28(5): 537-543. |
| [3] | 梁锦文. 位场向下延拓的正则化方法. 地球物理学报 , 1989, 32(5): 600–608. Liang J W. Downward continuation of regularization methods for potential fields. Chinese J. Geophys (in Chinese) , 1989, 32(5): 600-608. |
| [4] | Pawlowski R S. Preferential continuation for potential-field anomaly enhancement. Geophysics , 1995, 60(2): 390-398. DOI:10.1190/1.1443775 |
| [5] | 宁津生, 汪海洪, 罗志才. 基于多尺度边缘约束的重力场信号的向下延拓. 地球物理学报 , 2005, 48(1): 63–68. Ning J S, Wang H H, Luo Z C. Application of the multiscale edges in the downward continuation of gravity signal. Chinese J. Geophys (in Chinese) , 2005, 48(1): 63-68. |
| [6] | 王延忠, 熊光楚. 位场向下延拓组合滤波器的设计和应用. 地球物理学报 , 1985, 28(5): 537–543. Wang Y Z, Xiong G C. A combined filter used for the downward continuation of the potential field. Chinese J. Geophys (in Chinese) , 1985, 28(5): 537-543. |
| [7] | 毛小平, 吴蓉元, 曲赞. 频率域位场下延的振荡机制及消除方法. 石油地球物理勘探 , 1998, 33(2): 230–237. Mao X P, Wu R Y, Qu Z. Oscillation mechanism of potential field downward continuation in frequency domain and its elimination. Oil Geophysical Prospecting (in Chinese) , 1998, 33(2): 230-237. |
| [8] | 侯重初. 位场的频率域向下延拓方法. 物探与化探 , 1982, 6(1): 33–40. Hou Z C. Potential field downward continuation in frequency domain. Geophysical and Geochemical Exploration (in Chinese) , 1982, 6(1): 33-40. |
| [9] | Fedi M, Florio G. A stable downward continuation by using the ISVD method. Geophys. J. Int , 2002, 151(1): 146-156. DOI:10.1046/j.1365-246X.2002.01767.x |
| [10] | 杨文采. 用于位场数据处理的广义反演技术. 地球物理学报 , 1986, 29(3): 283–291. Yang W C. A generalized inversion technique for potential field data processing. Chinese J. Geophys (in Chinese) , 1986, 29(3): 283-291. |
| [11] | 陈生昌, 肖鹏飞. 位场向下延拓的波数域广义逆算法. 地球物理学报 , 2007, 50(6): 1819–1822. Chen S C, Xiao P F. Wavenumber domain generalized inverse algorithm for potential field downward continuation. Chinese J. Geophys (in Chinese) , 2007, 50(6): 1819-1822. |
| [12] | 徐世浙. 位场延拓的积分—迭代法. 地球物理学报 , 2006, 49(4): 1176–1182. Xu SZ. The integral-iteration method for continuation of potential fields. Chinese J. Geophys (in Chinese) , 2006, 49(4): 1176-1182. |
| [13] | 徐世浙. 迭代法与FFT法位场向下延拓效果的比较. 地球物理学报 , 2007, 50(1): 285–289. Xu S Z. A comparison of effects between the iteration method and FFT for downward continuation of potential fields. Chinese J. Geophys (in Chinese) , 2007, 50(1): 285-289. |
| [14] | Xu S Z, Yang J Y, Yang C F, et al. The iteration method for downward continuation of a potential field from a horizontal plane. Geophysical Prospecting , 2007, 55(6): 883-889. DOI:10.1111/gpr.2007.55.issue-6 |
| [15] | 王谦身. 重力学. 北京: 地震出版社, 2003 . Wang Q S. Gravitology (in Chinese). Beijing: Seismological Press, 2003 . |
| [16] | Li Y, Oldeuburg D W. Stable reduction to the pole at the magnetic equator. Geophysics , 2001, 66(2): 571-578. DOI:10.1190/1.1444948 |
| [17] | 侯重初. 补偿圆滑滤波方法. 石油物探 , 1981, 20(2): 22–29. Hou Z C. Filtering of smooth compensation. Geophysical Prospecting for Petroleum (in Chinese) (in Chinese) , 1981, 20(2): 22-29. |
| [18] | 张辉, 陈龙伟, 任治新, 等. 位场向下延拓迭代法收敛性分析及稳健向下延拓方法研究. 地球物理学报 , 2009, 52(4): 1107–1113. Zhang H, Chen L W, Ren Z X, et al. Analysis on convergence of iteration method for potential fields downward continuation and research on robust downward continuation method. Chinese J. Geophys (in Chinese) , 2009, 52(4): 1107-1113. |
| [19] | 郭志宏, 管志宁, 熊盛青. 长方体ΔT场及其梯度场无解析奇点理论表达式. 地球物理学报 , 2004, 47(6): 1131–1138. Guo Z H, Guan Z N, Xiong S Q. Cuboid ΔT and its gradient forward theoretical expressions without analytic odd points. Chinese J. Geophy (in Chinese) , 2004, 47(6): 1131-1138. |
| [20] | Granser H. Nonlinear inversion of gravity data using the Schmidt-Lichtenstein approach. Geophysics , 1987, 52(1): 88-93. DOI:10.1190/1.1442243 |
| [21] | Bott P. The use of rapid digital computing methods for direct gravity interpretation of sedimentary basins. Geophysical Journal of the Royal Astronomical Society , 1960, 3(1): 63-67. DOI:10.1111/gji.1960.3.issue-1 |
| [22] | 刘元龙, 王谦身. 用压缩质面法反演重力资料以估算地壳构造. 地球物理学报 , 1997, 20(1): 59–69. Liu Y L, Wang Q S. Inversion of gravity data by use of a method of "compressed mass plane" to estimate crustal structure. Chinese J. Geophys (in Chinese) , 1997, 20(1): 59-69. |
| [23] | Oldenburg D W. The inversion and interpretation of gravity anomalies. Geophysics , 1974, 39(4): 526-536. DOI:10.1190/1.1440444 |
| [24] | Parker R L. The rapid calculation of potential anomalies. Geophy Geophysical Journal of the Royal Astronomical Society , 1973, 31(4): 447-455. DOI:10.1111/j.1365-246X.1973.tb06513.x |
| [25] | Lefort J P, Agarwal B N P. Gravity and geomorphological evidence for a large crustal bulge cutting across Brittany (France): a tectonic response to the closure of the Bay of Biscay. Tectonophysics , 2000, 323(3-4): 149-162. DOI:10.1016/S0040-1951(00)00103-7 |
| [26] | Gomez O D, Agarwal B N P. 3DINVER. Computers and Geosciences , 2005, 31(4): 513-520. DOI:10.1016/j.cageo.2004.11.004 |
| [27] | 骆遥. 位场迭代法向下延拓的地球物理含义——以可下延异常逐次分离过程为例. 地球物理学进展 , 2011, 26(4): 1197–1200. Luo Y. Geophysical implications of the iteration method for downward continuation of potential field: a case by computing regional anomaly can be continued downward. Progress in Geophysics| (in Chinese) (in Chinese) , 2011, 26(4): 1197-1200. |
| [28] | 姚长利, 李宏伟, 郑元满, 等. 重磁位场转换计算中迭代法的综合分析与研究. 地球物理学报 , 2012, 55(6): 2062–2078. Yao C L, Li H W, Zheng Y M, et al. Research on iteration method using in potential field transformations. Chinese J. Geophys (in Chinese) , 2012, 55(6): 2062-2078. |