駱遙， 吳美平

1 國防科學技術大學 機電工程與自動化學院， 長沙　410073 2 中國國土資源航空物探遙感中心， 北京　100083

?

位場向下延拓的最小曲率方法

駱遙1, 2， 吳美平1

1 國防科學技術大學 機電工程與自動化學院， 長沙410073 2 中國國土資源航空物探遙感中心， 北京100083

摘要針對位場向下延拓的不適定性，我們將位場向下延拓視為向上延拓的反問題，提出以位場最小曲率作為約束條件來求解穩(wěn)定的下延位場.我們將剖面位場向上延拓表達式用傅里葉矩陣的形式表示，以矩陣乘法形式給出延拓的表達式，同時向待反演的下延位場引入最小曲率約束，得到向下延拓的最小曲率解，并利用正交變換給出了更為簡潔的頻率域解.隨后，利用Kronecker積將上述全部結果拓展至三維位場，給出了三維位場向下延拓的最小曲率解.此外，我們將位場數(shù)據(jù)的填充、擴充問題與向下延拓問題統(tǒng)籌考慮，提出一種新的向下延拓迭代格式，該算法面向實際資料處理需求、無須預擴充或填補數(shù)據(jù).下延迭代時，對原始數(shù)據(jù)直接向下延拓，而空白部分利用上一次下延位場估計的上延值替代其空白值并對其向下延拓，直至獲得最小曲率約束下穩(wěn)定的向下延拓結果.同時，我們也討論了利用改進L曲線和廣義交叉驗證(GCV)計算正則參數(shù)最優(yōu)估計的問題.對理論模型和實際航空重力資料進行了向下延拓檢驗，處理結果表明位場向下延拓的最小曲率方法解能滿足實際位場資料對向下延拓處理的需求，具有較高的下延精度.

關鍵詞向下延拓； 最小曲率； 位場； 數(shù)據(jù)空白； 正則化

1引言

位場(重力場或地磁場)向下延拓可以從數(shù)學上拉近觀測平面與場源間的距離，增強地質信息(改善信噪比)，不僅在地球物理勘探中具有重要作用，還是利用重力場或地磁場特征進行導航的一項關鍵技術(胡小平等，2013)，向下延拓技術在未爆軍火探測等領域也有重要應用.位場向下延拓是經典的不適定問題(Blakely，1996；管志寧，2005)，常規(guī)頻率域位場下延強烈的表現(xiàn)出對數(shù)據(jù)高頻部分的振蕩，因而難于實際應用.盡管眾多學者對位場向下延拓進行了大量細致的研究，但由于天然位場的頻帶之寬，各種向下延拓方法對不同資料適應性的差異很大.此外，幾乎所有涉及頻率域位場轉換的研究中都假設輸入的位場數(shù)據(jù)是純粹的數(shù)學矩陣或向量，既要數(shù)據(jù)無空白，又要滿足快速傅里葉變換(FFT)對數(shù)據(jù)長度的要求.實測位場資料很難滿足上述條件，必須對資料空白部分進行填充、對數(shù)據(jù)進行擴充(擴邊)，這將直接關系到位場向下延拓的精度(王萬銀等，2009).數(shù)據(jù)的填充與擴充是實際位場資料延拓研究中難以回避的問題，需要加以解決.

正則化方法是求解不適定問題較完備且行之有效的方法.針對位場向下延拓的不適定性，我們以位場最小曲率作為約束條件求解穩(wěn)定的下延位場，同時將位場數(shù)據(jù)的填充、擴充問題與向下延拓問題統(tǒng)籌考慮，提出一種將多種需求一并處理的向下延拓技術.具體研究中我們將位場向下延拓問題視為向上延拓的反問題，首先從二維位場剖面的向下延拓問題出發(fā)，將頻率域向上延拓表達式用傅里葉變換矩陣的形式表示；利用矩陣乘法形式給出空間域向上延拓表達式的基礎上，推導出向下延拓的頻率域表達式及位場向下延拓的迭代解.同時，通過對下延位場引入最小曲率約束，給出了空間域向下延拓的最小曲率解，并最終通過正交變換給出更簡潔的頻率域解.隨后，利用Kronecker積將上述結果拓展至三維位場，并給出了頻率域三維位場向下延拓的最小曲率解.最后，我們討論了頻率域延拓中遇到的擴充數(shù)據(jù)、填補數(shù)據(jù)空白區(qū)等實際問題，給出了一種無須預擴充數(shù)據(jù)、填補數(shù)據(jù)空白的位場向下延拓技術，這對解決實際位場資料延拓問題具有重要意義.

2向下延拓技術方法

2.1二維位場延拓

位場解析延拓是典型的狄利克萊問題，根據(jù)重磁勘探教科書(如Blakely，1996；管志寧，2005；等)，兩個高度層的剖面二維位場滿足

(1)

其中U(x,h)對應高高度層的位場，U(x,0)對應低高度層的位場，高度差h>0為常量.對其進行傅里葉變換，可以得到位場延拓的頻率域表達式

FT[U(x,h)]=exp(-2πh|kx|)FT[U(x,0)],

(2)

(3)

(4)

于是(2)式表達的位場延拓關系可寫成矩陣形式，有

(5)

其中Uh=[U(1,h),…,U(n,h),…,U(N,h)]T，U0=[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處的波數(shù).若高高度位場Uh為觀測位場，計算低高度位場U0的問題轉化為求解目標函數(shù)

(6)

的極小化問題，其中‖·‖2表示L-2范數(shù).該目標函數(shù)最小化的解為

(7)

其中上標T表示轉置，即

(8)

將傅里葉變換矩陣用FT[·]表示，有

FT[U(x,0)]=exp(2πh|kx|)FT[U(x,h)].(9)

(9)

式就是位場向下延拓的頻率域表達式.(6)式的極小化問題也可采用數(shù)值方法求解，如采用超松弛迭代法，有

(10)

(11)

則

(12)將傅里葉變換矩陣用FT[·]表示，有

(13)即為通常意義下的位場向下延拓的正則化方法(欒文貴，1983；梁錦文，1989；Cooper，2004；陳生昌和肖鵬飛，2007).上述正則化方法的本質是引入了約束

(14)

即(12)或(13)式是目標函數(shù)

(15)

極小化的解.因此，位場向下延拓能否穩(wěn)定的關鍵在于對下延位場施加的約束.實際上(15)式引入的約束僅體現(xiàn)模型最小，未反映位場本身的特性，并非“最優(yōu)”.類似的，在位場反演中附加到場源的約束條件除模型最小外，還包括對模型的一階偏導數(shù)，但將這類約束條件引入位場向下延拓問題中也不“恰當”，一階偏導數(shù)最小是約束物性分布不存在間斷，對位場本身而言仍不充分.

在處理位場資料插值或網(wǎng)格化問題時，Briggs(1974)引入了約束

(16)對離散數(shù)據(jù)進行網(wǎng)格化——即最小曲率網(wǎng)格化方法，Smith和Wessel(1990)曾對其進行改進.最小曲率網(wǎng)格化方法是目前位場資料網(wǎng)格化的標準方法.因此，利用最小曲率的方式約束位場向下延拓最為恰當，能夠繼承位場資料的網(wǎng)格化特征.對于剖面位場最小曲率的約束是

(17)

將其寫成離散形式，有

(18)

其中下標N表示方陣階數(shù)，該方陣具體為

(19)

為了減少頻率域位場轉換中Gibbs效應，DN算子采用了周期性邊界條件.因此，在最小曲率條件約束下，求解下延位場的目標函數(shù)可表示為

在選擇耕作方法時，綜合考慮當季作物的效應及其對土壤肥力的影響。例如，免耕法具有提高播種質量、保持土壤水分、減少土壤侵蝕、降低生產成本等優(yōu)點。但是采用免耕法，肥料只能施于土壤表層，肥料的利用率低，特別是有機肥料增肥土壤的作用大大降低。多點生產實踐證明，不管什么類型的土壤，長期免耕必然造成土壤耕層變淺，影響作物根系下扎，不利于作物抗倒伏。要推廣科學施肥技術，優(yōu)化肥料運籌，改進施肥方法，發(fā)揮養(yǎng)分協(xié)同作用，提高肥料利用率，提升農作物產量。

(20)

其中α>0.對目標函數(shù)即(20)式極小化，有

(21)

由于U0的左乘方陣是對角占優(yōu)矩陣，可以直接求逆矩陣來解(21)式，也可以采用共軛梯度法求解(21)式.顯而易見，隨著觀測數(shù)據(jù)的增多，求解該方程將變得極困難，不得不考慮類似反演問題中抽稀觀測數(shù)據(jù)的做法(Foksetal., 2013).目前，航空磁測中單一測區(qū)覆蓋面積最大達1.14×106km2(熊盛青等，2001)，面對大數(shù)據(jù)處理量的實際情況，需考慮避免直接求解超大型方程組.我們考慮對(19)式算子進行正交分解(Strang, 1999)，有

(22)

其中

(23)

故(21)式可記為

(24)

其解為

(25)

FT[U(x,0)]=

×FT[U(x,h)],

(26)

其中n=1,2,…,N是U(x,h)經離散傅里葉變換后對應序列中元素的序號.利用(26)式即可通過FFT直接計算在最小曲率條件約束下的下延位場，從而避免直接用(21)式在空間域中進行求解.

2.2三維位場延拓

類似的，重磁勘探教科書(Blakely，1996；管志寧，2005；等)也給了不同高度層間三維位場的關系，有U(x,y,h)=

對其進行傅里葉變換，可以得到位場延拓的頻率域表達式

FT[U(x,y,h)]=

(28)

其中kx、ky對應頻率域波數(shù).按照前面的思路，可將(28)式表示成矩陣乘法的關系，但經二維離散傅里葉變換后的位場網(wǎng)格是M×N的矩陣，我們需要將FT[U(x,y,·)]表示成列向量形式，同時借助Kronecker積將二維傅里葉變換表示為矩陣乘法的形式.

設二維矩陣f={fm n∈f|1≤m≤M,1≤n≤N}，經二維傅里葉變換得到的矩陣為F={Fm n∈F|1≤m≤M,1≤n≤N}，將其表示為列向量形式，有vec(f)=[f11,f21,…,fM1,f12,f22,…,fM2,…,

f1N,f2N,…,fMN]T,

(29)

vec(F)=[F11,F21,…,FM1,F12,F22,…,FM2,…,

F1N,F2N,…,FMN]T,

(30)

其中vec(·)是將矩陣轉為列向量的變換.于是離散傅里葉變換的矩陣表達為

(31)

其中?表示Kronecker積，有

(32)

(33)

故(28)式可寫成矩陣形式為

(34)

其中Uh、U0分別為高高度層和低高度層的位場網(wǎng)格，并以矩陣形式表示；

…,

kxm表示位場經傅里葉變換后矩陣對應第m行處的波數(shù)，kyn表示對應第n列處的波數(shù)，顯然對角陣ΛMN主對角線上元素構成的向量實際上是位場經傅里葉變換后矩陣所對應的頻率域向上延拓因子的列向量變換.類似的，我們取目標函數(shù)

?WM)vec(U0)‖2,

(35)

其最小化的解為

(36)

即

?WMvec(Uh),

(37)

將傅里葉變換矩陣用FT[·]表示，有

(38)

?WMvec(U0)‖2+α‖DMNvec(U0)‖2,

(39)

其中

(40)

對該目標函數(shù)極小化，有

(41)

于是，下延位場計算可通過求解該方程獲得.類似的物性反演中多采取共軛梯度法，但仍難以處理稍大的數(shù)據(jù)量，例如處理1024×1024大小的普通網(wǎng)格數(shù)據(jù)時，(41)式方程組的系數(shù)矩陣中元素個數(shù)就達1012量級，顯然難以承受！出于實際考慮，我們利用頻率域方法求解(41)式.

(42)

于是

(43)

其中

(44)

因此(41)式的解為

(45)

將傅里葉變換矩陣用FT[·]表示，有

(46)

其中m=1,2,…,M和n=1,2,…,N分別是U(x,y,h)經離散傅里葉變換后對應矩陣元素的行序號和列序號.(46)式即最小曲率約束下位場向下延拓的頻率域解，該解與(41)式的空間域方法完全等效.

2.3正則參數(shù)選擇

前面提出的基于最小曲率的位場向下延拓方法本質上屬于正則化反演，正則參數(shù)的選擇強烈地影響著向下延拓的結果，以往許多文獻(欒文貴，1983；梁錦文，1989；Cooper，2004；陳生昌和肖鵬飛，2007)并沒有專門對此進行討論，其正則參數(shù)的選取多屬經驗性的.無論是二維位場還是三維位場其解析延拓都可表達為

(47)

其中向量d表示觀測位場，向量m表示下延位場，矩陣G表示向上延拓算子.求解下延位場就是求目標函數(shù)

(48)

的極小值問題，其中矩陣H表示前面的最小曲率約束算子(剖面數(shù)據(jù)時為DN，網(wǎng)格數(shù)據(jù)時為DMN).確定α的方法之一是L曲線法，但L曲線隅角的確定比較復雜，我們建議通過極小化泛函

(50)

由于φ(α)能直接在頻率域中求解，正則參數(shù)的確定是很容易的.除L曲線方法外，也可以利用廣義交叉驗證(GCV)方法確定正則化參數(shù)(Golub et al., 1979).廣義交叉驗證是通過極小化泛函

GCV(α)=n‖d-Gm(α)‖2/[n-Tr(A(α))]2,

(51)

確定正則參數(shù)，其中n是觀測位場的數(shù)量，Tr(·)表示矩陣的跡，

(52)

(51)式的分子(不含n)同(49)式中乘法的第一項因子相同，是下延位場上延至原平面的恢復誤差，可以用頻率域上延算法快速計算.但對矩陣A的跡求解卻較為困難，這里我們借鑒Garcia(2010)關于跡的處理，利用頻率域與空間域的對應關系進行快速計算.對于剖面二維位場，有

(53)

對于網(wǎng)格三維位場，有

(54)

于是，通過(53)或(54)式可以快速計算(51)式所需的跡，從而獲得對正則參數(shù)的最優(yōu)估計，即

(55)

2.4數(shù)據(jù)填充與擴充處理

目前幾乎所有的頻率域位場轉換都借助快速傅里葉變換(FFT)實現(xiàn)，而FFT對數(shù)據(jù)長度是有一定要求的，通常要求剖面數(shù)據(jù)的點數(shù)或網(wǎng)格數(shù)據(jù)的行、列數(shù)為合數(shù)，一般要求為以2為基(即2的自然數(shù)次冪).插值后的剖面位場或網(wǎng)格化后的位場通常難以滿足上述要求，這就要求將剖面數(shù)據(jù)或網(wǎng)格數(shù)據(jù)向外擴充至FFT所要求的數(shù)據(jù)長度，即位場數(shù)據(jù)的擴充(擴邊)問題.我們以FFT中最常見的2n為例，如果位場剖面的長度為1666個點，進行位場轉換時就需要將剖面兩端各擴充191個數(shù)據(jù)點，即擴充后數(shù)據(jù)長度為2048即211；再如位場網(wǎng)格的節(jié)點數(shù)為365×666，進行位場轉換時就需要將網(wǎng)格首、末行的兩端分別向外擴充73行和74行，再將擴充后網(wǎng)格的首、末列兩端向外各擴充179列，擴充后的網(wǎng)格大小為512×1024即29×210.此外，由于多種因素某些區(qū)域缺乏實測位場數(shù)據(jù)，經過網(wǎng)格化后這些區(qū)域也無法形成相應的內插數(shù)據(jù)進而存在空白，位場轉換中還要求將上述空白用合理的數(shù)據(jù)(稱該數(shù)據(jù)為“假值”)進行填充(替代).填充與擴充中的間斷現(xiàn)象將引起Gibbs效應，但許多理論研究者卻并沒有高度重視該問題.實際上，位場數(shù)據(jù)中空白值在航空磁測或航空重力測量中廣泛存在卻是不爭的事實.測量中測線方向通常并非南北向或東西向，而是垂直于主地質構造走向，如果采用正網(wǎng)格(網(wǎng)格列方向為南北向)方式進行網(wǎng)格化勢必產生大量數(shù)據(jù)空值(Zhang et al., 2012)，即使采用斜網(wǎng)格(網(wǎng)格列方向為測線方向)方式進行網(wǎng)格化，也會因測區(qū)形狀的不規(guī)則而產生一定數(shù)量的空值，有時還因特殊原因造成測區(qū)的部分數(shù)據(jù)缺失，部分歷史航磁資料曾出現(xiàn)測區(qū)中間空白的情況.在航磁編圖中數(shù)據(jù)空值問題更加突出，例如全國航磁編圖中的位場轉換中一個極重要的問題就是如何對網(wǎng)格空值進行假值填充處理，航磁編圖數(shù)據(jù)在國境以外存在大量網(wǎng)格空值，即使利用國外磁測資料進行填充，但由于中國藏南地區(qū)、中國臺灣海峽等尚存的航空磁測空白區(qū)，網(wǎng)格空值問題近乎無法避免！

除位場轉換處理本身，數(shù)據(jù)填充與擴充處理是直接影響位場轉換精度的重要因素，例如曾小牛等(2011)利用改進的積分迭代法處理重力資料時邊界處的誤差顯著增大，明顯是數(shù)據(jù)擴邊不理想所致.由于某些位場轉換研究本身就極富挑戰(zhàn)性，很難再將填充與擴充處理統(tǒng)籌考慮，多數(shù)位場轉換研究中僅討論位場轉換本身的精度，對數(shù)據(jù)空值填充或擴邊的問題鮮有涉獵.姚長利等(2003)在低緯度化極時曾指出數(shù)據(jù)擴充對計算造成的影響不能忽視，但沒給出進一步措施.張輝等(2009)結合正則化方法和維納濾波方法研究波數(shù)域延拓算法時曾指出其算法存在邊界效應問題，需引入擴邊處理，但未對擴邊算法進行研究.王萬銀等(2009)提出了用最小曲率方式對位場轉換數(shù)據(jù)進行空值填補與擴邊，其研究表明利用最小曲率方式進行數(shù)據(jù)擴充優(yōu)于余弦方法，適于對位場數(shù)據(jù)擴充和填補空白.我們在前面的討論中并沒有涉及到數(shù)據(jù)填充與擴充處理，這就需要預先對位場數(shù)據(jù)中空值進行填補和擴充.似乎這樣便解決位場下延中遇到的實際問題，但仔細考慮卻存在矛盾，因為本文提出的向下延拓方法本身即基于最小曲率方式，而最小曲率方法卻多用于對未知位場的插值處理(Briggs, 1974；Smith and Wessel, 1990)，這似乎表明沒有必要單獨對數(shù)據(jù)進行填充和擴充處理.

基于解決位場數(shù)據(jù)存在假值與擴邊問題的實際需求，結合最小曲率方法的特性，我們將前面所述的位場向下延拓方法進一步擴展.這里我們借鑒了Garcia(2010)對數(shù)據(jù)空值估計的思想，提出下延處理中對非空數(shù)據(jù)直接使用，而空白部分則利用上一次下延位場(估計)的上延值代替并進行迭代，具體的計算格式為

(56)

3理論模型計算

圖1　二維重力模型及其向下延拓結果對比

圖2　地面(h=0 km)和高空(h=1 km)處的模型磁異常

圖3　正則參數(shù)曲線(a)與向下延拓結果(b)

圖4　高空(h=1 km)處的模型磁異常(白色區(qū)域為三處數(shù)據(jù)空白區(qū))

圖6　迭代過程使用的正則參數(shù)及相應迭代步的均方誤差

圖5　向下延拓結果(a)及其上延恢復的觀測異常(b)其中(a)和(b)對應的理論模型異常數(shù)據(jù)以等值線的彩色填充形式疊加作為比照.

4實際資料處理

為了檢驗位場最小曲率向下延拓方法在實際數(shù)據(jù)處理中的實用性，我們采用美國地質勘探局(USGS)2006年在阿富汗實測的航空重力資料進行數(shù)據(jù)處理檢驗.該航空重力測量數(shù)據(jù)經各項處理并最終歸算至離地7000 m的高度，形成網(wǎng)格距為1 km的布格重力異常.圖7給出了布格重力異常圖(等值線間隔10 mGal)，其中角圖中藍色區(qū)塊為航空重力測量覆蓋范圍，紅點為1966—1967年地面重力測量點.由圖7可知該測區(qū)極不規(guī)則且不同區(qū)塊測線方向不同(測線間距也有所不同)，這些特征是航空物探作業(yè)中的典型特點，但多數(shù)向下延拓文獻卻刻意回避這些特點，而是將類似情況的數(shù)據(jù)裁剪成特殊的無空值矩形網(wǎng)格進行處理，以避免數(shù)據(jù)填充或擴充(擴邊)問題.我們將圖7的布格重力異常延拓至地面(下延距離7000 m)，由于圖7的網(wǎng)格含有(大量)空白數(shù)據(jù)，這里采用(56)式給出的迭代格式進行求解(迭代過程中正則參數(shù)逐漸縮小并最終等于預估值α=0.65)，圖8給出了最終下延至地面的布格重力異常(等值線間隔10 mGal).可以看出，向下延拓7000 m后得到的重力異常(圖8)沒有出現(xiàn)虛假的振蕩，在保持原始異常變化趨勢的同時增強了局部異常.我們將圖8的布格重力異常同地面布格重力異常進行比較，二者的外符合精度(郭志宏等, 2008)為7.5 mGal，反映二者差異較大，但卻難以反映延拓處理的精度.USGS為了編圖的需要將地面布格異常歸算至距地面7000 m的高空，歸算后的布格異常同航空重力異常比較，二者的外符合精度也僅為6.4 mGal，可見該地區(qū)的地面重力資料和航空重力資料本身符合的就不甚理想，這也是向下延拓至地面的布格異常同地面實測結果相差較大的原因.如果忽略二者間資料本身的差異，向下延拓引起的誤差可近視為兩種外符合精度之差即1.1 mGal，相對于下延7000 m(7倍網(wǎng)格距)而言，延拓精度是相當可觀的.此外，我們將圖8的布格異常上延恢復至原高度比較其與圖7的差異，二者相差小于0.36 mGal的網(wǎng)格點占99%，其差值的均方差僅0.1 mGal，也進一步證明了向下延拓的數(shù)值精度.

圖7　航空布格重力異常數(shù)據(jù)(歸算離地高度7 km)

圖8　向下延拓至地面的航空布格重力異常

5結論

我們將位場向下延拓視為向上延拓處理的反問題，并以位場最小曲率條件為正則化約束，反演位場數(shù)據(jù)向下延拓的最小曲率解.由于使用了傅里葉變換矩陣，我們提出的向下延拓方法既可以在空間域求解，也可在頻率域通過快速傅里葉變換直接求解.此外，我們討論了頻率域延拓中遇到的擴充數(shù)據(jù)、填補數(shù)據(jù)空白區(qū)兩項實際問題，給出了無須預擴充數(shù)據(jù)、填補數(shù)據(jù)空白的位場向下延拓算法.針對我們提出的頻率域向下延拓的最小曲率解，我們分別利用二維重力模型和三維磁異常模型進行了檢驗，并特別針對三維位場數(shù)據(jù)中含數(shù)據(jù)空白的情況進行了專門的檢驗，無論是向下延拓的數(shù)值精度還是對空白區(qū)的恢復都取得了良好的理論模型效果.同時，我們選取實測的航空重力資料進行向下延拓處理.該測區(qū)極不規(guī)則，并含有大塊的數(shù)據(jù)空白區(qū)域，甚至個別數(shù)據(jù)僅由單條測線網(wǎng)格化而來.在未對該資料進行填充、擴充等預處理情況下直接將其向下延拓至地面，處理結果表明該方法具有較高的向下延拓精度，可在實際位場資料處理應用中發(fā)揮積極作用.

致謝圖7使用了美國地質勘探局(USGS)公布的航空重力數(shù)據(jù)和地面重力數(shù)據(jù)，在此表示感謝.

References

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.Blakely R J. 1996. Potential Theory in Gravity and Magnetic Applications. New York: Cambridge University Press.

Briggs I C. 1974. Machine contouring using minimum curvature.Geophysics, 39(1): 39-48, doi: 10.1190/1.1440410.

Chen S C, Xiao P F. 2007. Wavenumber domain generalized inverse algorithm for potential field downward continuation.ChineseJ.Geophys. (in Chinese), 50(6): 1816-1822.

Cooper G. 2004. The stable downward continuation of potential field data.ExplorationGeophysics, 35(4)：260-265, doi: 10.1071/EG04260.

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.

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.

Gao Y W, Luo Y, Wen W. 2012. The compensation method for downward continuation of potential field from horizontal plane and its application.ChineseJ.Geophys. (in Chinese), 55(8): 2747-2756, doi: 10.6038/j.issn.0001-5733.2012.08.026.

Garcia D. 2010. Robust smoothing of gridded data in one and higher dimensions with missing values.ComputationalStatistics&DataAnalysis, 54(4): 1167-1178, doi: 10.1016/j.csda.2009.09.020.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.

Guan Z N. 2005. Magnetic Field and Magnetic Exploration (in Chinese). Beijing: Geological Publishing House.

Guo Z H, Guan Z N, Xiong S Q. 2004. Cuboid ΔTand its gradient forward theoretical expressions without analytic odd points.ChineseJ.Geophys. (in Chinese), 47(6): 1131-1138.

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.ChineseJ.Geophys. (in Chinese), 51(5): 1538-1543.

Hu X P, Wu M P, Mu H, et al. 2013. Technologies on Underwater Geomagnetic Field Navigation (in Chinese). Beijing: National Defense Industry Press.

Liang J W. 1989. Downward continuation of regularization methods for potential fields.ChineseJ.Geophys. (ActaGeophysicaSinica) (in Chinese), 32(5): 600-608.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.ChineseJ.Geophys. (in Chinese), 52(6): 1599-1605, doi: 10.3969/j.issn.0001-5733.2009.06.022.Luan W G. 1983. The stabilized algorithm of the analytic continuation for the potential field.ChineseJ.Geophys. (ActaGeophysicaSinica) (in Chinese), 26(3): 263-274.Luo Y, Yao C L. 2007. Theoretical study on cuboid magnetic field and its gradient expression without analytic singular point.OilGeophysicalProspecting(in Chinese), 42(6): 714-719.

Smith W H F, Wessel P. 1990. Gridding with continuous curvature splines in tension.Geophysics, 55(3): 293-305, doi: 10.1190/1.1442837.

Strang G. 1999. The discrete cosine transform.SIAMReview, 41(1): 135-147, doi: 10.1137/S0036144598336745.

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.ProgressinGeophys. (in Chinese), 24(4): 1327-1338, doi: 10.3969/j.issn.1004-2903.2009.04.022.

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.

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.GeophysicalProspecting, 55(6): 883-889, doi: 10.1111/j.1365-2478.2007.00634.x.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.ChineseJ.Geophys. (in Chinese), 46(5): 690-696.

Yao C L, Li H W, Zheng Y M, et al. 2012. Research on iteration method using in potential field transformations.ChineseJ.Geophys. (in Chinese), 55(6): 2062-2078, doi: 10.6038/j.issn.0001-5733.2012.06.028.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.AppliedGeophysics, 9(4): 378-390, doi: 10.1007/s11770-012-0349-x. 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.ChineseJ.Geophys. (in Chinese), 52(4): 1107-1113, doi: 10.3969/j.issn.0001-5733.2009.04.028.

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.ChineseJ.Geophys. (in Chinese), 54(11): 2943-2950, doi: 10.3969/j.issn.0001-5733.2011.11.024.

附中文參考文獻

陳生昌, 肖鵬飛. 2007. 位場向下延拓的波數(shù)域廣義逆算法. 地球物理學, 50(6): 1816-1822.

高玉文, 駱遙, 文武. 2012. 補償向下延拓方法研究及應用. 地球物理學報, 55(8): 2747-2756, doi: 10.6038/j.issn.0001-5733.2012.08.026.

管志寧. 2005. 地磁場與磁力勘探. 北京: 地質出版社.

郭志宏, 管志寧, 熊盛青. 2004. 長方體ΔT場及其梯度場無解析奇點理論表達式. 地球物理學報, 47(6): 1131-1138.

郭志宏, 熊盛青, 周堅鑫等. 2008. 航空重力重復線測試數(shù)據(jù)質量評價方法研究. 地球物理學報, 51(5): 1538-1543.

胡小平, 吳美平, 穆華等. 2013. 水下地磁導航技術. 北京: 國防工業(yè)出版社.

梁錦文. 1989. 位場向下延拓的正則化方法. 地球物理學報, 32(5): 600-608.

劉東甲, 洪天求, 賈志海等. 2009. 位場向下延拓的波數(shù)域迭代法及其收斂性. 地球物理學報, 52(6): 1599-1605, doi: 10.3969/j.issn.0001-5733.2009.06.022.

欒文貴. 1983. 場位解析延拓的穩(wěn)定化算法. 地球物理學報, 26(3): 263-274.

駱遙, 姚長利. 2007. 長方體磁場及其梯度無解析奇點表達式理論研究. 石油地球物理勘探, 42(6): 714-719.

王萬銀, 邱之云, 劉金蘭等. 2009. 位場數(shù)據(jù)處理中的最小曲率擴邊和補空方法研究. 地球物理學進展, 24(4): 1327-1338, doi: 10.3969/j.issn.1004-2903.2009.04.022.

熊盛青, 周伏洪, 姚正煦等. 2001. 青藏高原中西部航磁調查. 北京: 地質出版社.

姚長利, 管志寧, 高德章等. 2003. 低緯度磁異?；瘶O方法——壓制因子法. 地球物理學報, 46(5): 690-696.

姚長利, 李宏偉, 鄭元滿等. 2012. 重磁位場轉換計算中迭代法的綜合分析與研究. 地球物理學報, 55(6): 2062-2078, doi: 10.6038/j.issn.0001-5733.2012.06.028.

張輝, 陳龍偉, 任治新等. 2009. 位場向下延拓迭代法收斂性分析及穩(wěn)健向下延拓方法研究. 地球物理學報, 52(4): 1107-1113, doi: 10.3969/j.issn.0001-5733.2009.04.028.

曾小牛, 李夕海, 劉代志等. 2011. 積分迭代法的正則性分析及其最優(yōu)步長的選擇. 地球物理學報, 54(11): 2943-2950, doi: 10.3969/j.issn.0001-5733.2011.11.024.

(本文編輯何燕)

基金項目國家自然科學基金(61174206)、國家高技術研究發(fā)展計劃(863計劃)(2011AA060501，2013AA063901，2013AA063905)和中國地質調查(1212011120189)項目資助．

作者簡介駱遙，男，1982年生，博士研究生，高級工程師，主要從事地球物理探測等相關技術研究．E-mail：geophy@vip.qq.com

doi:10.6038/cjg20160120 中圖分類號P631

收稿日期2015-01-22，2015-12-09收修定稿

Minimum curvature method for downward continuation of potential field data

LUO Yao1, 2, WU Mei-Ping1

1CollegeofMechatronicsEngineeringandAutomation,NationalUniversityofDefenseTechnology,Changsha410073,China2ChinaAeroGeophysicalSurveyandRemoteSensingCenterforLandandResources,Beijing100083,China

AbstractDownward 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.

KeywordsDownward continuation; Minimum curvature; Potential field; Missing data; Regularization

駱遙，吳美平. 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.ChineseJ.Geophys. (in Chinese),59(1):240-251,doi:10.6038/cjg20160120.