地球物理学报  2016, Vol. 59 Issue (1): 240-251   PDF    
引用本文
骆遥, 吴美平. 2016. 位场向下延拓的最小曲率方法. 地球物理学报, 59(1): 240-251, doi: 10.6038/cjg20160120   
Luo Y, Wu M P. 2016. Minimum curvature method for downward continuation of potential field data. Chinese J. Geophysics, 59(1): 240-251, doi: 10.6038/cjg20160120   
位场向下延拓的最小曲率方法
骆遥1,2, 吴美平1    
1. 国防科学技术大学机电工程与自动化学院, 长沙 410073;
2. 中国国土资源航空物探遥感中心, 北京 100083
收稿日期 2015-01-22, 2015-12-09 收修定稿
基金项目: 国家自然科学基金(61174206)、国家高技术研究发展计划(863计划)(2011AA060501,2013AA063901,2013AA063905)和中国地质调查(1212011120189)项目资助.
作者简介: 骆遥,男,1982年生,博士研究生,高级工程师,主要从事地球物理探测等相关技术研究.E-mail:geophy@vip.qq.com
摘要: 针对位场向下延拓的不适定性,我们将位场向下延拓视为向上延拓的反问题,提出以位场最小曲率作为约束条件来求解稳定的下延位场.我们将剖面位场向上延拓表达式用傅里叶矩阵的形式表示,以矩阵乘法形式给出延拓的表达式,同时向待反演的下延位场引入最小曲率约束,得到向下延拓的最小曲率解,并利用正交变换给出了更为简洁的频率域解.随后,利用Kronecker积将上述全部结果拓展至三维位场,给出了三维位场向下延拓的最小曲率解.此外,我们将位场数据的填充、扩充问题与向下延拓问题统筹考虑,提出一种新的向下延拓迭代格式,该算法面向实际资料处理需求、无须预扩充或填补数据.下延迭代时,对原始数据直接向下延拓,而空白部分利用上一次下延位场估计的上延值替代其空白值并对其向下延拓,直至获得最小曲率约束下稳定的向下延拓结果.同时,我们也讨论了利用改进L曲线和广义交叉验证(GCV)计算正则参数最优估计的问题.对理论模型和实际航空重力资料进行了向下延拓检验,处理结果表明位场向下延拓的最小曲率方法解能满足实际位场资料对向下延拓处理的需求,具有较高的下延精度.
关键词: 向下延拓     最小曲率     位场     数据空白     正则化    
Minimum curvature method for downward continuation of potential field data
LUO Yao1,2, WU Mei-Ping1    
1. College of Mechatronics Engineering and Automation, National University of Defense Technology, Changsha 410073, China;
2. China Aero Geophysical Survey and Remote Sensing Center for Land and Resources, Beijing 100083, China
Abstract: Downward continuation calculates the potential field that is closer to the source of anomalies is a powerful but unstable tool in data processing. To obtain a reasonable downward continued result, we formulate the downward continuation process as an inverse problem of upward continuation, and introduce a discretized smoothing continuous curvature spline regularization approach to solve it. The proposed method, which is based on the discrete Fourier transform(DFT) matrix, allows robust smoothing of downward continued field data. In the study, we formulate the upward continuation of 2-D potential field data as matrix multiplication in the form d=Gm represented by DFT matrix. Then, the inverse problem is formulated as minimizing a total objective function consisting of total squared curvature of downward continued field data and data misfit. Because the number of data is often large, calculating the matrix equation is a major computational load of the downward continuation. We greatly simplify and solve it by means of the DFT. Then we extend the method to 3-D potential field by using the Kronecker product, and obtain the minimum curvature solution of downward continued field data in the frequency domain. As the results of downward continuation are strongly influenced by the smoothing parameter, automatic choice of the amount of regularization parameter is carried out using the method of L-curve and generalized cross validation(GCV). In addition, we consider the problems on data missing and grid expansion in the downward continuation of potential field. An iterative robust scheme of downward continuation is then proposed to deal with data expansion or replacing dummy values with interpolated values. To test the performance of the method, we employ the synthetic and real airborne gravity data for downward continuation. Data processing results on both the synthetic and real data confirm that the performance of the proposed technique can satisfy the need of real data processing for downward continuation of the potential field.
Key words: Downward continuation     Minimum curvature     Potential field     Missing data     Regularization    
1 引言

位场(重力场或地磁场)向下延拓可以从数学上拉近观测平面与场源间的距离，增强地质信息(改善信噪比)，不仅在地球物理勘探中具有重要作用，还是利用重力场或地磁场特征进行导航的一项关键技术(胡小平等，2013)，向下延拓技术在未爆军火探测等领域也有重要应用.位场向下延拓是经典的不适定问题(Blakely，1996；管志宁，2005)，常规频率域位场下延强烈的表现出对数据高频部分的振荡，因而难于实际应用.尽管众多学者对位场向下延拓进行了大量细致的研究，但由于天然位场的频带之宽，各种向下延拓方法对不同资料适应性的差异很大.此外，几乎所有涉及频率域位场转换的研究中都假设输入的位场数据是纯粹的数学矩阵或向量，既要数据无空白，又要满足快速傅里叶变换(FFT)对数据长度的要求.实测位场资料很难满足上述条件，必须对资料空白部分进行填充、对数据进行扩充(扩边)，这将直接关系到位场向下延拓的精度(王万银等，2009).数据的填充与扩充是实际位场资料延拓研究中难以回避的问题，需要加以解决.

正则化方法是求解不适定问题较完备且行之有效的方法.针对位场向下延拓的不适定性，我们以位场最小曲率作为约束条件求解稳定的下延位场，同时将位场数据的填充、扩充问题与向下延拓问题统筹考虑，提出一种将多种需求一并处理的向下延拓技术.具体研究中我们将位场向下延拓问题视为向上延拓的反问题，首先从二维位场剖面的向下延拓问题出发，将频率域向上延拓表达式用傅里叶变换矩阵的形式表示；利用矩阵乘法形式给出空间域向上延拓表达式的基础上，推导出向下延拓的频率域表达式及位场向下延拓的迭代解.同时，通过对下延位场引入最小曲率约束，给出了空间域向下延拓的最小曲率解，并最终通过正交变换给出更简洁的频率域解.随后，利用Kronecker积将上述结果拓展至三维位场，并给出了频率域三维位场向下延拓的最小曲率解.最后，我们讨论了频率域延拓中遇到的扩充数据、填补数据空白区等实际问题，给出了一种无须预扩充数据、填补数据空白的位场向下延拓技术，这对解决实际位场资料延拓问题具有重要意义.

2 向下延拓技术方法 2.1 二维位场延拓

位场解析延拓是典型的狄利克莱问题，根据重磁勘探教科书(如Blakely，1996；管志宁，2005；等)，两个高度层的剖面二维位场满足

其中U(x，h)对应高高度层的位场，U(x，0)对应低高度层的位场，高度差h>0为常量.对其进行傅里叶变换，可以得到位场延拓的频率域表达式

其中FT[·]表示傅里叶变换，kx对应频率域波数.实际位场转换中采用傅里叶变换的离散形式即快速傅里叶变换(FFT)，可将傅里叶变换写成矩阵形式

其中 f =[f1，f2，…，fN]T为观测数据，其傅里叶变换 为 F =[F1，F2，…，FN]T， W N={wmn∈ W N｜ 1≤m≤N，

1)/N)

i2=-1)为N阶酉阵，即WN* WN= I N(I N是N阶单位矩阵，上标*表示共轭转置)，则傅里叶逆变换表示为

于是(2)式表达的位场延拓关系可写成矩阵形式，有

其中 U h=[U(1，h)，…，U(n，h)，…，U(N，h)]T，U 0=[U(1，0)，…，U(n，0)，…，U(N，0)]T，U(n，·)对应离散序列在n(n=1，2，…，N)处的位场， Λ N=diag{exp(-2πh｜ kx1 ｜)，…，exp(-2πh｜ kxn ｜)，…，exp(-2πh｜ kxN ｜)}，kxn对应离散序列在n处的波数.若高高度位场 U h为观测位场，计算低高度位场 U 0的问题转化为求解目标函数

的极小化问题，其中‖ · ‖2表示L-2范数.该目标函数最小化的解为

其中上标T表示转置，即

将傅里叶变换矩阵用FT[·]表示，有

(9)式就是位场向下延拓的频率域表达式.(6)式的极小化问题也可采用数值方法求解，如采用超松弛迭代法，有

其中 U0k表示第k次迭代的低高度位场，ω是松弛因子，通常取(0，2).实际上(10)式即为位场向下延拓的迭代解(Xu et al.，2007)，当迭代次数k较大时(10)式与(9)式则在数值上等价，有关迭代法讨论有众多文献参考，如刘东甲等(2009)及姚长利等(2012)的论述，这里不再重复.要获得稳定的位场向下延拓就必须考虑向目标函数引入正则化项.例如，可直接向(7)式中对角阵 Λ NT Λ N中引入微小增量保持稳定，有

则

将傅里叶变换矩阵用FT[·]表示，有

即为通常意义下的位场向下延拓的正则化方法(栾文贵，1983；梁锦文，1989；Cooper，2004；陈生昌和肖鹏飞，2007).上述正则化方法的本质是引入了约束

即(12)或(13)式是目标函数

极小化的解.因此，位场向下延拓能否稳定的关键在于对下延位场施加的约束.实际上(15)式引入的约束仅体现模型最小，未反映位场本身的特性，并非“最优”.类似的，在位场反演中附加到场源的约束条件除模型最小外，还包括对模型的一阶偏导数，但将这类约束条件引入位场向下延拓问题中也不“恰当”，一阶偏导数最小是约束物性分布不存在间断，对位场本身而言仍不充分.

在处理位场资料插值或网格化问题时，Briggs(1974)引入了约束

对离散数据进行网格化——即最小曲率网格化方法，Smith和Wessel(1990)曾对其进行改进.最小曲率网格化方法是目前位场资料网格化的标准方法.因此，利用最小曲率的方式约束位场向下延拓最为恰当，能够继承位场资料的网格化特征.对于剖面位场最小曲率的约束是

将其写成离散形式，有

其中下标N表示方阵阶数，该方阵具体为

为了减少频率域位场转换中Gibbs效应，D N算子采用了周期性边界条件.因此，在最小曲率条件约束下，求解下延位场的目标函数可表示为

其中α>0.对目标函数即(20)式极小化，有

由于 U 0的左乘方阵是对角占优矩阵，可以直接求逆矩阵来解(21)式，也可以采用共轭梯度法求解(21)式.显而易见，随着观测数据的增多，求解该方程将变得极困难，不得不考虑类似反演问题中抽稀观测数据的做法(Foks et al.，2013).目前，航空磁测中单一测区覆盖面积最大达1.14×106 km2(熊盛青等，2001)，面对大数据处理量的实际情况，需考虑避免直接求解超大型方程组.我们考虑对(19)式算子进行正交分解(Strang，1999)，有

其中

故(21)式可记为

其解为

将傅里叶变换矩阵用FT[·]表示，有

其中n=1，2，…，N是U(x，h)经离散傅里叶变换后对应序列中元素的序号.利用(26)式即可通过FFT直接计算在最小曲率条件约束下的下延位场，从而避免直接用(21)式在空间域中进行求解.

2.2 三维位场延拓

类似的，重磁勘探教科书(Blakely，1996；管志宁，2005；等)也给了不同高度层间三维位场的关系，有

对其进行傅里叶变换，可以得到位场延拓的频率域表达式

其中kx、ky对应频率域波数.按照前面的思路，可将(28)式表示成矩阵乘法的关系，但经二维离散傅里叶变换后的位场网格是M×N的矩阵，我们需要将FT[U(x，y，·)]表示成列向量形式，同时借助Kronecker积将二维傅里叶变换表示为矩阵乘法的形式.

设二维矩阵 f ={fmn∈ f ｜ 1≤m≤M，1≤n≤N}，经二维傅里叶变换得到的矩阵为 F ={Fmn∈ F ｜ 1≤m≤M，1≤n≤N}，将其表示为列向量形式，有

其中vec(·)是将矩阵转为列向量的变换.于是离散傅里叶变换的矩阵表达为

其中表示Kronecker积，有

其中矩阵 B 第m行第n列的元素用Bmn表示.显然 W N W M是(M·N)阶酉阵，于是二维离散傅里叶变换逆变换的矩阵表达为

故(28)式可写成矩阵形式为

其中 U h、 U 0分别为高高度层和低高度层的位场网格，并以矩阵形式表示；

kxm表示位场经傅里叶变换后矩阵对应第m行处的波数，kyn表示对应第n列处的波数，显然对角阵 Λ MN主对角线上元素构成的向量实际上是位场经傅里叶变换后矩阵所对应的频率域向上延拓因子的列向量变换.类似的，我们取目标函数

其最小化的解为

即

将傅里叶变换矩阵用FT[·]表示，有

类似剖面位场的处理，我们很容易获得位场向下延拓的正则化解(栾文贵，1983；梁锦文，1989；Cooper，2004；陈生昌和肖鹏飞，2007)或迭代解(Xu et al.，2007)，此处不再重复.我们仅重点讨论附加(16)式约束下的最小曲率解. Δ2U0的离散形式可以表示为(I N D M+ D N I M)vec( U 0)，因此在最小曲率条件约束下求解下延位场的目标函数为

其中

对该目标函数极小化，有

于是，下延位场计算可通过求解该方程获得.类似的物性反演中多采取共轭梯度法，但仍难以处理稍大的数据量，例如处理1024×1024大小的普通网格数据时，(41)式方程组的系数矩阵中元素个数就达1012量级，显然难以承受！出于实际考虑，我们利用频率域方法求解(41)式.

由于 I N= WN* I N W N且 D M= W *M Γ M W M，根据Kronecker积的性质，有

于是

其中

因此(41)式的解为

将傅里叶变换矩阵用FT[·]表示，有

其中m=1，2，…，M和n=1，2，…，N分别是U(x，y，h)经离散傅里叶变换后对应矩阵元素的行序号和列序号.(46)式即最小曲率约束下位场向下延拓的频率域解，该解与(41)式的空间域方法完全等效.

2.3 正则参数选择

前面提出的基于最小曲率的位场向下延拓方法本质上属于正则化反演，正则参数的选择强烈地影响着向下延拓的结果，以往许多文献(栾文贵，1983；梁锦文，1989；Cooper，2004；陈生昌和肖鹏飞，2007)并没有专门对此进行讨论，其正则参数的选取多属经验性的.无论是二维位场还是三维位场其解析延拓都可表达为

其中向量 d 表示观测位场，向量 m 表示下延位场，矩阵 G 表示向上延拓算子.求解下延位场就是求目标函数

的极小值问题，其中矩阵 H 表示前面的最小曲率约束算子(剖面数据时为 D N，网格数据时为 D MN).确定α的方法之一是L曲线法，但L曲线隅角的确定比较复杂，我们建议通过极小化泛函

确定正则参数，其中 m (α)代表不同α下向下延拓位场 m 的最小曲率解.正则参数的最优估计 为

由于φ(α)能直接在频率域中求解，正则参数的确定是很容易的.除L曲线方法外，也可以利用广义交叉验证(GCV)方法确定正则化参数(Golub et al.，1979).广义交叉验证是通过极小化泛函

确定正则参数，其中n是观测位场的数量，Tr(·)表示矩阵的迹，

(51)式的分子(不含n)同(49)式中乘法的第一项因子相同，是下延位场上延至原平面的恢复误差，可以用频率域上延算法快速计算.但对矩阵 A 的迹求解却较为困难，这里我们借鉴Garcia(2010)关于迹的处理，利用频率域与空间域的对应关系进行快速计算.对于剖面二维位场，有

对于网格三维位场，有

于是，通过(53)或(54)式可以快速计算(51)式所需的迹，从而获得对正则参数的最优估计，即

2.4 数据填充与扩充处理

目前几乎所有的频率域位场转换都借助快速傅里叶变换(FFT)实现，而FFT对数据长度是有一定要求的，通常要求剖面数据的点数或网格数据的行、列数为合数，一般要求为以2为基(即2的自然数次幂).插值后的剖面位场或网格化后的位场通常难以满足上述要求，这就要求将剖面数据或网格数据向外扩充至FFT所要求的数据长度，即位场数据的扩充(扩边)问题.我们以FFT中最常见的2n为例，如果位场剖面的长度为1666个点，进行位场转换时就需要将剖面两端各扩充191个数据点，即扩充后数据长度为2048即211；再如位场网格的节点数为365×666，进行位场转换时就需要将网格首、末行的两端分别向外扩充73行和74行，再将扩充后网格的首、末列两端向外各扩充179列，扩充后的网格大小为512×1024即29×210.此外，由于多种因素某些区域缺乏实测位场数据，经过网格化后这些区域也无法形成相应的内插数据进而存在空白，位场转换中还要求将上述空白用合理的数据(称该数据为“假值”)进行填充(替代).填充与扩充中的间断现象将引起Gibbs效应，但许多理论研究者却并没有高度重视该问题.实际上，位场数据中空白值在航空磁测或航空重力测量中广泛存在却是不争的事实.测量中测线方向通常并非南北向或东西向，而是垂直于主地质构造走向，如果采用正网格(网格列方向为南北向)方式进行网格化势必产生大量数据空值(Zhang et al.，2012)，即使采用斜网格(网格列方向为测线方向)方式进行网格化，也会因测区形状的不规则而产生一定数量的空值，有时还因特殊原因造成测区的部分数据缺失，部分历史航磁资料曾出现测区中间空白的情况.在航磁编图中数据空值问题更加突出，例如全国航磁编图中的位场转换中一个极重要的问题就是如何对网格空值进行假值填充处理，航磁编图数据在国境以外存在大量网格空值，即使利用国外磁测资料进行填充，但由于中国藏南地区、中国台湾海峡等尚存的航空磁测空白区，网格空值问题近乎无法避免！

除位场转换处理本身，数据填充与扩充处理是直接影响位场转换精度的重要因素，例如曾小牛等(2011)利用改进的积分迭代法处理重力资料时边界处的误差显著增大，明显是数据扩边不理想所致.由于某些位场转换研究本身就极富挑战性，很难再将填充与扩充处理统筹考虑，多数位场转换研究中仅讨论位场转换本身的精度，对数据空值填充或扩边的问题鲜有涉猎.姚长利等(2003)在低纬度化极时曾指出数据扩充对计算造成的影响不能忽视，但没给出进一步措施.张辉等(2009)结合正则化方法和维纳滤波方法研究波数域延拓算法时曾指出其算法存在边界效应问题，需引入扩边处理，但未对扩边算法进行研究.王万银等(2009)提出了用最小曲率方式对位场转换数据进行空值填补与扩边，其研究表明利用最小曲率方式进行数据扩充优于余弦方法，适于对位场数据扩充和填补空白.我们在前面的讨论中并没有涉及到数据填充与扩充处理，这就需要预先对位场数据中空值进行填补和扩充.似乎这样便解决位场下延中遇到的实际问题，但仔细考虑却存在矛盾，因为本文提出的向下延拓方法本身即基于最小曲率方式，而最小曲率方法却多用于对未知位场的插值处理(Briggs，1974；Smith and Wessel，1990)，这似乎表明没有必要单独对数据进行填充和扩充处理.

基于解决位场数据存在假值与扩边问题的实际需求，结合最小曲率方法的特性，我们将前面所述的位场向下延拓方法进一步扩展.这里我们借鉴了Garcia(2010)对数据空值估计的思想，提出下延处理中对非空数据直接使用，而空白部分则利用上一次下延位场(估计)的上延值代替并进行迭代，具体的计算格式为

其中U0k表示第k次迭代的下延位场；FT[·]表示傅里叶变换，FT-1{·}表示其逆变换，upc(·，h)表示将位场上延h高度的操作；Λ(α，h)表示下延距离为h正则化参数为α的频率域算子，即(26)式或(46)式中的频率域下延因子；γ是权重因子，当观测位场在该节点处空白时取0，否则取1.(56)式的直观解释是下延位场在最小曲率条件约束下满足观测位场延拓关系的同时，保证空白区的上延估计同观测位场连续.具体迭代计算中，下延位场的初始值选零，并满足FFT对位场扩充的数据长度要求.若位场转换中要求边界趋于零，可令待求解位场最佳估计或最佳上延估计在数据边界取零.同时，为了保证在迭代过程中能有效地对数据空值和扩充区域进行估计，要求正则参数在迭代初期要适当大，以充分保证空白处位场的光滑，随着迭代次数的增加要逐步缩小正则参数，直至正则参数接近广义交叉验证(GCV)或L曲线方法对正则参数的预估计，迭代次数以能最佳恢复出观测位场为宜.

3 理论模型计算

为了验证提出的位场下延方法，我们分别就二维位场和三维位场进行理论模型计算.由于剖面位场向下延拓实际应用较少，我们仅展示对(26)式中频率域位场下延算子的理论模型计算部分.选取Pateka等(2012)使用的剖面位场模型数据进行检验，该模型为埋深5000 m、半径500 m的无限延伸水平圆柱产生的重力异常.图 1a和图 1b分别给出了z=0 m(地面)和z=-4000 m(向下为负)处的重力异常理论曲线，图 1c给出使用REGCONT代 码(Pateka et al.，2012)向下延拓得到z=-4000 m 的重力异常曲线.为了检验我们给出的下延算子，我们直接将REGCONT代码中的下延算子修改为(26)式中频率域下延算子，其他不做任何修改(包括正则参数选取的C-norm准则)，图 1d给出了最小曲率方式的向下延拓结果.通过对图 1d和图 1b及图 1c比较，可以看出最小曲率方式能获得较理想的下延结果，但同Tikhonov正则化下延结果相比振荡略显增大，这主要是由于C-norm准则选取的是位场模型最小的情况下的正则参数而非位场曲率最小，C-norm判据下的正则参数仅为8.5×10-10(图 1d)，相比而言极小，故图 1d较图 1c表现出的振荡情况略明显，但图 1d的下延结果可以确认(26)式中频率域下延算子的有效性.

图 1 二维重力模型及其向下延拓结果对比Fig. 1 Comparison of synthetic gravitational field and its downward continuation

我们重点讨论对三维位场的理论模型计算，因为二维位场仅是三维情况的特例.这里采用郭志宏等(2004)给出三维磁性模型作为标准模型，高玉文等(2012)曾利用该模型进行向下延拓检验.该模型为长方体，长宽高分别为8 km、4 km、4 km，中心埋深3 km，磁化强度1 A·m-1，磁化方向为偏角30°倾角50°，地磁场偏角10°倾角60°.我们利用长方体磁场无解析奇点表达式(骆遥和姚长利，2007)计算距离地面1 km高空处的ΔT场作为延拓数据，下延1 km(10倍网格距)计算地面处的磁异常，并以地面处的磁异常正演结果作为评价标准进行检验.图 2a和图 2b分别给出地面处和1 km高空处磁异常的正演结果.利用(46)式计算将图 2b下延至地面处，其中扩边方式直接采用Cooper等编写的taper2d函数，该扩边方式曾在变磁倾角化极(Cooper and Cowan，2005)和Potensoft位场数据处理软件(Arisoy and Dikmen，2011)中应用.正则参数则分别采用改进L曲线方式和GCV方式确定，图 3a给出了GCV曲线和改进的L曲线，从图 3a可知两种方式确定的正则参数曲线形态一致仅幅度不同，确定最优正则参数均为=0.102，这反映了GCV方式和L曲线方式在正则参数确定中的一致性.图 3b是最优正则参数下的向下延拓结果，可以看出图 3b给出的下延结果同图 2a的理论结果非常一致，二者差值均方差为0.6 nT，仅在0 nT等值线的边缘略有畸变，说明了向下延拓的可靠性.

图 2 地面(h=0 km)和高空(h=1 km)处的模型磁异常Fig. 2 Total field anomalies of synthetic model at different observation planes

图 3 正则参数曲线(a)与向下延拓结果(b)Fig. 3 The regularization parameter curves (a) and corresponding evaluated regularized solution of downward continuation (b)

为了进一步检验对数据填充与扩充的处理效果，我们将图 2b中“挖去”三处数据作为模型数据(图 4)，图 4中三处空白分别选取在正异常梯度带 处、正负异常交接处及边部异常平缓处，具有一定的型性.利用(56)式将图 4数据向下延拓1 km，图 5a给出了下延至地面的ΔT等值线(图 5没有将数据空白处或扩充处的下延估计值去掉，故数据范围扩大)，此外将下延后的磁异常重新上延回h=1 km的原平面(图 5b)，进一步比较数据恢复误差.向下延拓时正则参数从2等速率下降至最佳估计=0.102后保持恒定，图 6 给出了该过程使用的正则参数.同时，图 6也给出了迭代过程中恢复观测磁异常的误差变化情况，以说明计算的收敛性.为了进一步说明向下延拓的精度，无论是图 5a向下延拓结果还是图 5b恢复的观测数据都叠加了相应模型数据等值线的彩色填充以供对比.可以看出，图 5a的等值线与理论模型的等值线彩色填充区块叠合严密，仅-80 nT和-60 nT等值线的最南部边缘与色块区域略有错动，将下延结果上延回原平面恢复得到 的磁异常等值线(图 5b)与理论模型的等值线彩色填充区块相叠合则严丝合缝，进一步说明了向下延拓同数据填充、扩充结合所取得的明显试验效果.

图 4 高空(h=1 km)处的模型磁异常(白色区域为三处数据空白区)Fig. 4 Synthetic magnetic anomalies with missing values

图 5 向下延拓结果(a)及其上延恢复的观测异常(b)
其中(a)和(b)对应的理论模型异常数据以等值线的彩色填充形式叠加作为比照.Fig. 5 The downward continuation field (a) and the restored result (b) by continuing (a) to original elevation
All the calculated contour maps added the image layers of corresponding synthetic data in colors.

图 6 迭代过程使用的正则参数及相应迭代步的均方误差Fig. 6 The regularization parameter used in each steps and the corresponding RMS error
4 实际资料处理

为了检验位场最小曲率向下延拓方法在实际数据处理中的实用性，我们采用美国地质勘探局(USGS)2006年在阿富汗实测的航空重力资料进行数据处理检验.该航空重力测量数据经各项处理并最终归算至离地7000 m的高度，形成网格距为1 km的布格重力异常.图 7给出了布格重力异常图(等值线间隔10 mGal)，其中角图中蓝色区块为航空重力测量覆盖范围，红点为1966—1967年地面重力测量点.由图 7可知该测区极不规则且不同区块测线方向不同(测线间距也有所不同)，这些特 征是航空物探作业中的典型特点，但多数向下延拓文献却刻意回避这些特点，而是将类似情况的数据裁剪成特殊的无空值矩形网格进行处理，以避免数据填充或扩充(扩边)问题.我们将图 7的布格重力异常延拓至地面(下延距离7000 m)，由于图 7的网格含有(大量)空白数据，这里采用(56)式给出的迭代格式进行求解(迭代过程中正则参数逐渐缩小并最终等于预估值α=0.65)，图 8给出了最终下延至地面的布格重力异常(等值线间隔10 mGal).可以看出，向下延拓7000 m后得到的重力异常(图 8)没有出现虚假的振荡，在保持原始异常变化趋势的同时增强了局部异常.我们将图 8的布格重力异常同地面布格重力异常进行比较，二者的外符合精度(郭志宏等，2008)为7.5 mGal，反映二者差异较大，但却难以反映延拓处理的精度.USGS为了编图的需要将地面布格异常归算至距地面7000 m的高空，归算后的布格异常同航空重力异常比较，二者的外符合精度也仅为6.4 mGal，可见该地区的地面重力资料和航空重力资料本身符合的就不甚理想，这也是向下延拓至地面的布格异常同地面实测结果相差较大的原因.如果忽略二者间资料本身的差异，向下延拓引起的误差可近视为两种外符合精度之差即1.1 mGal，相对于下延7000 m(7倍网格距)而言，延拓精度是相当可观的.此外，我们将图 8的布格异常上延恢复至原高度比较其与图 7的差异，二者相差小于0.36 mGal的网格点占99%，其差值的均方差仅0.1 mGal，也进一步证明了向下延拓的数值精度.

图 7 航空布格重力异常数据(归算离地高度7 km)Fig. 7 Airborne Bouguer gravity anomaly map at 7 km above ground

图 8 向下延拓至地面的航空布格重力异常Fig. 8 Airborne Bouguer gravity anomaly map continued to the ground
5 结论

我们将位场向下延拓视为向上延拓处理的反问题，并以位场最小曲率条件为正则化约束，反演位场数据向下延拓的最小曲率解.由于使用了傅里叶变换矩阵，我们提出的向下延拓方法既可以在空间域求解，也可在频率域通过快速傅里叶变换直接求解.此外，我们讨论了频率域延拓中遇到的扩充数据、填补数据空白区两项实际问题，给出了无须预扩充数据、填补数据空白的位场向下延拓算法.针对我们提出的频率域向下延拓的最小曲率解，我们分别利用二维重力模型和三维磁异常模型进行了检验，并特别针对三维位场数据中含数据空白的情况进行了专门的检验，无论是向下延拓的数值精度还是对空白区的恢复都取得了良好的理论模型效果.同时，我们选取实测的航空重力资料进行向下延拓处理.该测区极不规则，并含有大块的数据空白区域，甚至个别数据仅由单条测线网格化而来.在未对该资料进行填充、扩充等预处理情况下直接将其向下延拓至地面，处理结果表明该方法具有较高的向下延拓精度，可在实际位场资料处理应用中发挥积极作用.

致谢 图 7使用了美国地质勘探局(USGS)公布的航空重力数据和地面重力数据，在此表示感谢.
参考文献
[1] Arisoy M Ö, Dikmen Ü. 2011. Potensoft:MATLAB-based software for potential field data processing, modeling and mapping. Computers & Geosciences, 37(7):935-942, doi:10.1016/j.cageo.2011.02.008.
[2] Blakely R J. 1996. Potential Theory in Gravity and Magnetic Applications. New York:Cambridge University Press.
[3] Briggs I C. 1974. Machine contouring using minimum curvature. Geophysics, 39(1):39-48, doi:10.1190/1.1440410.
[4] 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.
[5] Cooper G. 2004. The stable downward continuation of potential field data. Exploration Geophysics, 35(4):260-265, doi:10.1071/EG04260.
[6] Cooper G R J, Cowan D R. 2005. Differential reduction to the pole. Computers & Geosciences, 31(8):989-999, doi:10.1016/j.cageo.2005.02.005.
[7] Foks N L, Krahenbuhl R, Li Y G. 2013. Adaptive sampling of potential-field data:A direct approach to compressive inversion. Geophysics, 79(1):IM1-IM9, doi:10.1190/GEO2013-0087.1.
[8] Gao Y W, Luo Y, Wen W. 2012. The compensation method for downward continuation of potential field from horizontal plane and its application. Chinese J. Geophys.(in Chinese), 55(8):2747-2756, doi:10.6038/j.issn.0001-5733.2012.08.026.
[9] Garcia D. 2010. Robust smoothing of gridded data in one and higher dimensions with missing values. Computational Statistics & Data Analysis, 54(4):1167-1178, doi:10.1016/j.csda.2009.09.020.
[10] Golub G H, Heath M, Wahba G. 1979. Generalized cross-validation as a method for choosing a good ridge parameter. Technometrics, 21(2):215-223, doi:10.1080/00401706.1979.10489751.
[11] Guan Z N. 2005. Magnetic Field and Magnetic Exploration(in Chinese). Beijing:Geological Publishing House.
[12] Guo Z H, Guan Z N, Xiong S Q. 2004. Cuboid ΔT and its gradient forward theoretical expressions without analytic odd points. Chinese J. Geophys.(in Chinese), 47(6):1131-1138.
[13] Guo Z H, Xiong S Q, Zhou J X, et al. 2008. The research on quality evaluation method of test repeat lines in airborne gravity survey. Chinese J. Geophys.(in Chinese), 51(5):1538-1543.
[14] Hu X P, Wu M P, Mu H, et al. 2013. Technologies on Underwater Geomagnetic Field Navigation(in Chinese). Beijing:National Defense Industry Press.
[15] Liang J W. 1989. Downward continuation of regularization methods for potential fields. Chinese J. Geophys.(Acta Geophysica Sinica)(in Chinese), 32(5):600-608.
[16] 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.
[17] Luan W G. 1983. The stabilized algorithm of the analytic continuation for the potential field. Chinese J. Geophys.(Acta Geophysica Sinica)(in Chinese), 26(3):263-274.
[18] Luo Y, Yao C L. 2007. Theoretical study on cuboid magnetic field and its gradient expression without analytic singular point. Oil Geophysical Prospecting(in Chinese), 42(6):714-719.
[19] Pateka R, Karcol R, Kunirák D, et al. 2012. REGCONT:A Matlab based program for stable downward continuation of geophysical potential fields using Tikhonov regularization. Computers & Geosciences, 49:278-289, doi:10.1016/j.cageo.2012.06.010.
[20] Smith W H F, Wessel P. 1990. Gridding with continuous curvature splines in tension. Geophysics, 55(3):293-305, doi:10.1190/1.1442837.
[21] Strang G. 1999. The discrete cosine transform. SIAM Review, 41(1):135-147, doi:10.1137/S0036144598336745.
[22] 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 Geophys.(in Chinese), 24(4):1327-1338, doi:10.3969/j.issn.1004-2903.2009.04.022.
[23] Xiong S Q, Zhou F H, Yao Z X, et al. 2001. Aeromagnetic Survey in Central and Western Qinghai-Tibet Plateau(in Chinese). Beijing:Geological Publishing House.
[24] Xu S Z, Yang J Y, Yang C F, et al. 2007. The iteration method for downward continuation of a potential field from a horizontal plane. Geophysical Prospecting, 55(6):883-889, doi:10.1111/j.1365-2478.2007.00634.x.
[25] Yao C L, Guan Z N, Gao D Z, et al. 2003. Reduction-to-the-pole of magnetic anomalies at low latitude with suppression filter. Chinese J. Geophys.(in Chinese), 46(5):690-696.
[26] Yao C L, Li H W, Zheng Y M, et al. 2012. Research on iteration method using in potential field transformations. Chinese J. Geophys.(in Chinese), 55(6):2062-2078, doi:10.6038/j.issn.0001-5733.2012.06.028.
[27] Zhang C, Yao C L, Xie Y M, et al. 2012. Reduction of distortion and improvement of efficiency for gridding of scattered gravity and magnetic data. Applied Geophysics, 9(4):378-390, doi:10.1007/s11770-012-0349-x.
[28] 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 J. Geophys.(in Chinese), 52(4):1107-1113, doi:10.3969/j.issn.0001-5733.2009.04.028.
[29] Zheng X N, Li X H, Liu D Z, et al. 2001. Regularization analysis of integral iteration method and the choice of its optimal step-length. Chinese J. Geophys.(in Chinese), 54(11):2943-2950, doi:10.3969/j.issn.0001-5733.2011.11.024.
[30] 陈生昌, 肖鹏飞. 2007. 位场向下延拓的波数域广义逆算法. 地球物理学, 50(6):1816-1822.
[31] 高玉文, 骆遥, 文武. 2012. 补偿向下延拓方法研究及应用. 地球物理学报, 55(8):2747-2756, doi:10.6038/j.issn.0001-5733.2012.08.026.
[32] 管志宁. 2005. 地磁场与磁力勘探. 北京:地质出版社.
[33] 郭志宏, 管志宁, 熊盛青. 2004. 长方体ΔT场及其梯度场无解析奇点理论表达式. 地球物理学报, 47(6):1131-1138.
[34] 郭志宏, 熊盛青, 周坚鑫等. 2008. 航空重力重复线测试数据质量评价方法研究. 地球物理学报, 51(5):1538-1543.
[35] 胡小平, 吴美平, 穆华等. 2013. 水下地磁导航技术. 北京:国防工业出版社.
[36] 梁锦文. 1989. 位场向下延拓的正则化方法. 地球物理学报, 32(5):600-608.
[37] 刘东甲, 洪天求, 贾志海等. 2009. 位场向下延拓的波数域迭代法及其收敛性. 地球物理学报, 52(6):1599-1605, doi:10.3969/j.issn.0001-5733.2009.06.022.
[38] 栾文贵. 1983. 场位解析延拓的稳定化算法. 地球物理学报, 26(3):263-274.
[39] 骆遥, 姚长利. 2007. 长方体磁场及其梯度无解析奇点表达式理论研究. 石油地球物理勘探, 42(6):714-719.
[40] 王万银, 邱之云, 刘金兰等. 2009. 位场数据处理中的最小曲率扩边和补空方法研究. 地球物理学进展, 24(4):1327-1338, doi:10.3969/j.issn.1004-2903.2009.04.022.
[41] 熊盛青, 周伏洪, 姚正煦等. 2001. 青藏高原中西部航磁调查. 北京:地质出版社.
[42] 姚长利, 管志宁, 高德章等. 2003. 低纬度磁异常化极方法--压制因子法. 地球物理学报, 46(5):690-696.
[43] 姚长利, 李宏伟, 郑元满等. 2012. 重磁位场转换计算中迭代法的综合分析与研究. 地球物理学报, 55(6):2062-2078, doi:10.6038/j.issn.0001-5733.2012.06.028.
[44] 张辉, 陈龙伟, 任治新等. 2009. 位场向下延拓迭代法收敛性分析及稳健向下延拓方法研究. 地球物理学报, 52(4):1107-1113, doi:10.3969/j.issn.0001-5733.2009.04.028.
[45] 曾小牛, 李夕海, 刘代志等. 2011. 积分迭代法的正则性分析及其最优步长的选择. 地球物理学报, 54(11):2943-2950, doi:10.3969/j.issn.0001-5733.2011.11.024.