基于曲波變換的航空電磁數(shù)據(jù)調(diào)平方法研究

高玲琦， 殷長春* ， 王寧， 劉云鶴， 蘇揚， 熊彬

1 吉林大學地球探測科學與技術(shù)學院， 長春 130026 2 桂林理工大學地球科學學院， 桂林 541006

0 引言

航空電磁法是基于直升機、固定翼飛機等移動搭載平臺的地球物理方法.它利用不接地回線進行電磁發(fā)射并通過磁傳感器接收電磁信號，實現(xiàn)對地下目體的探測.航空電磁具有效率高、成本低的特點，特別適合高山、沙漠、湖泊沼澤、森林覆蓋等地形地質(zhì)條件復雜地區(qū)的資源探測.目前航空電磁法已在礦產(chǎn)資源、地下水和地熱、環(huán)境工程、城市地下空間、地質(zhì)災害監(jiān)測等領(lǐng)域發(fā)揮了重要的作用(Reed，1981；Wolfgram and Golden，2001；Haas，2004；Smith et al.，2006；Viezzoli et al.，2010；Wiederhold et al.，2010；Minsley et al.，2012).

由于采用飛機平臺搭載電磁設備，航空電磁數(shù)據(jù)中會產(chǎn)生很多地面電磁法遇不到的噪聲和誤差.最典型包括由飛機飛行狀態(tài)不同導致電磁系統(tǒng)發(fā)生變化，從而產(chǎn)生的調(diào)平誤差.調(diào)平誤差與隨機噪聲存在本質(zhì)區(qū)別，調(diào)平誤差存在特定的物理成因，不具隨機性，因此需要采取特殊方法對其進行處理.Huang和Fraser(1999)將調(diào)平誤差分為測區(qū)調(diào)平誤差、測線調(diào)平誤差以及短波長調(diào)平誤差.測區(qū)調(diào)平誤差主要由不同架次之間電磁系統(tǒng)物理條件變化引起，在數(shù)據(jù)中表現(xiàn)為區(qū)塊狀的異常抬升.測線調(diào)平誤差是由不同測線間由于溫度、壓力、風速和風向等條件發(fā)生變化，導致儀器狀態(tài)發(fā)生變化而產(chǎn)生的(Huang，2008；Beiki et al.，2010；Siemon,2009).航空電磁(特別是直升機航空電磁)通常采用相鄰測線來回飛行的“S”形觀測方式，因此某一條測線上電磁系統(tǒng)面向太陽(或者迎風)，相鄰測線上則背對太陽(或者避風)，造成相鄰測線間系統(tǒng)狀態(tài)發(fā)生變化，進而在電磁信號中產(chǎn)生條帶狀調(diào)平誤差(稱為窗簾效應).另外，由儀器噪聲、飛行速度突然改變及局部熱變化等因素還可引起短波長調(diào)平誤差，在數(shù)據(jù)中表現(xiàn)為局部干擾(Huang，2008).

航空電磁數(shù)據(jù)調(diào)平旨在消除調(diào)平誤差，改善數(shù)據(jù)質(zhì)量.傳統(tǒng)的位場數(shù)據(jù)調(diào)平大多基于位場隨高度變化緩慢的特點，采用切割線調(diào)平(Yarger et al.，1978；Bandy et al.，1990；Minty，1991；Luyendyk，1997；Ferraccioli et al.，1998).然而，航空電磁信號對飛行高度非常敏感，在切割線與測線交叉處，難以保證飛行高度一致，因此切割線調(diào)平無法進行電磁數(shù)據(jù)調(diào)平.Huang和Fraser(1999)提出了一種針對視電阻率數(shù)據(jù)的自動調(diào)平方法.該方法通過一維與二維數(shù)據(jù)非線性相關(guān)濾波器的組合實現(xiàn)數(shù)據(jù)調(diào)平.Fedi和Florio(2003)提出了一種基于小波變換的調(diào)平技術(shù)，然而該方法受到小波變換僅有三個方向特征的限制.Mauring和Kihle(2006)提出了基于移動微分中值濾波器的調(diào)平方法，該方法可以應用于非平行分布的測線數(shù)據(jù)調(diào)平.Huang(2008)提出了基于測線數(shù)據(jù)相關(guān)性的調(diào)平方法.該方法基于相鄰測線間信號具有較大的相關(guān)性，而調(diào)平誤差不具有相關(guān)性的特點，通過最小二乘擬合相鄰測線的數(shù)據(jù)差，實現(xiàn)航空電磁數(shù)據(jù)調(diào)平.Beiki等(2010)提出了基于一維和二維滑動窗口多項式擬合的航空電磁和航空磁測數(shù)據(jù)調(diào)平方法.Zhang等(2018)提出了基于主成分分析的調(diào)平方法并將其應用于航空電磁數(shù)據(jù)與航磁數(shù)據(jù)的調(diào)平.然而，上述方法均存在一定的局限性.其中，基于一維及二維濾波的調(diào)平方法本質(zhì)上對數(shù)據(jù)整體進行了平滑處理，容易導致高頻細節(jié)的丟失，而基于誤差多項式擬合的方法由于其階數(shù)限制，常常不能完美地擬合誤差，造成調(diào)平過度或殘余.

為改善上述方法的局限性，本文提出一種基于曲波變換的調(diào)平方法，通過曲波變換的多尺度和多方向特征，調(diào)平誤差得以在特定尺度和方向內(nèi)被精確擬合并去除，從而避免多項式在擬合誤差時的不準確性.另一方面，由于大多數(shù)尺度及方向未經(jīng)調(diào)平處理，故數(shù)據(jù)畸變小，且不容易導致高頻細節(jié)的丟失.在曲波變換等多尺度分析方法中，閾值法通過對曲波系數(shù)設置閾值實現(xiàn)數(shù)據(jù)去噪和數(shù)據(jù)重構(gòu)等.Donoho和Johnstone(1994)首先提出了用于小波變換重構(gòu)的硬閾值方法，之后，Donoho(1995)又提出了一種用于噪聲去除的軟閾值方法.目前兩種方法均獲得廣泛應用.本文實現(xiàn)航空電磁數(shù)據(jù)調(diào)平，我們首先對飛行高度進行歸一化處理(Huang，2008)，在此基礎(chǔ)上對數(shù)據(jù)進行多尺度分析，通過設置閾值進行方向濾波有效去除數(shù)據(jù)中的調(diào)平誤差.這種新調(diào)平方法與基于曲波變換數(shù)據(jù)去噪的不同之處在于：數(shù)據(jù)去噪僅需對信號進行多尺度分析，通過對高階尺度設置閾值即可去除隨機噪聲(王寧等，2020)；然而，調(diào)平誤差具有方向性，但不具有隨機性，因此需要在特定方向和尺度上對其進行識別，進而通過設置閾值法去除調(diào)平誤差.下面我們首先介紹曲波變換及其多尺度和多方向性，進而通過對調(diào)平誤差進行分析，討論如何通過在不同尺度和不同方向上設置閾值對調(diào)平誤差進行去除.最后，我們通過對合成數(shù)據(jù)以及愛爾蘭Waterford地區(qū)頻率域航空電磁數(shù)據(jù)的調(diào)平驗證本文調(diào)平方法的有效性.

1 基于曲波變換的數(shù)據(jù)調(diào)平方法

1.1 曲波變換

曲波變換是在小波變換基礎(chǔ)上發(fā)展而來，同時繼承了小波變換的多尺度性質(zhì).除此之外，曲波變換還具有多方向性，這使得我們能夠利用曲波變換更好地表征二維或更高維的信號.本文將利用這些曲波變換特性實現(xiàn)航空電磁數(shù)據(jù)調(diào)平.為了方便理解，我們將首先介紹曲波變換的基本理論和離散算法.由于航空電磁數(shù)據(jù)是離散的，我們主要描述離散形式的曲波變換.

1.1.1 波數(shù)域中的笛卡爾窗函數(shù)

在二維波數(shù)域中，離散曲波變換通過笛卡爾窗函數(shù)平滑地提取不同尺度和不同方向上的數(shù)據(jù)信息(Candès et al.，2006).笛卡爾窗函數(shù)由徑向窗與角度窗一起共同構(gòu)建.其中，徑向窗是一系列同心的正方形，可定義為

(1)

(2)

其中，ω=(ω1,ω2)為二維波數(shù)域中的變量，j為尺度參數(shù)，Φ由兩個一維低通濾波窗口的乘積定義，即：

Φj(ω1,ω2)=φ(2-jω1)φ(2-jω2)，

(3)

二維波數(shù)域是由所有徑向窗支撐而成，因此：

(4)

離散曲波變換中的角度窗Vj(ω)呈放射狀，可定義為

Vj(ω)=V(2?j/2_|ω2/ω1),

(5)

其中，|_j/2_|表示對j/2取整.在對徑向窗與角度窗分別進行定義后，離散曲波變換的笛卡爾窗函數(shù)可表示為

(6)

進而，我們引入等間隔的斜率：tanθl∶=l·2-?j/2_|,l=-2-?j/2_|,…,2?j/2_|-1，則(6)式可表示為

(7)

其中，剪切矩陣Sθ可表示為

(8)

圖1 離散曲波變換定義的楔形窗(參考Ying et al.，2005)Fig.1 Wedge-shaped window defined by discrete curvelet transform (Refer to Ying et al., 2005)

圖1展示了二維波數(shù)域中的窗函數(shù).其中，最內(nèi)層代表粗尺度.該尺度內(nèi)的曲波不具有方向性(相當于小波).在粗尺度之外，徑向窗呈同心的正方形狀，而角度窗呈圖中虛線所示的放射狀.圖中灰色的楔形區(qū)域由徑向窗與角度窗共同支撐，它包含了特定尺度和方向上的數(shù)據(jù)信息.

1.1.2 空間域中的曲波函數(shù)

離散曲波函數(shù)為連續(xù)曲波函數(shù)的近似，它保留了連續(xù)曲波函數(shù)的所有性質(zhì)(Ying et al.，2005).為了介紹離散曲波函數(shù)，我們首先定義連續(xù)曲波函數(shù).在連續(xù)曲波變換理論中，母曲波φj(x)由楔形窗的逆傅里葉變換得到，進而所有尺度為2-j的曲波均可由母曲波經(jīng)位移和旋轉(zhuǎn)得到(Candès et al.，2006).

(9)

其中，x=(x1,x2)為空間域自變量，旋轉(zhuǎn)矩陣定義為

(10)

由上面的討論，我們可以通過二維信號與曲波函數(shù)的內(nèi)積計算曲波系數(shù)，即：

(11)

通過二維傅里葉變換，我們也可以在波數(shù)域中計算曲波系數(shù)，即：

(12)

其中，Uj為連續(xù)曲波變換中的窗函數(shù)，可由φj的傅里葉變換得到.

在定義了連續(xù)曲波函數(shù)后，我們可以利用與(9)式相似的形式定義離散曲波函數(shù)：

(13)

其中，b為離散變量，可表示為

(14)

圖2展示了四個不同的曲波函數(shù)與其對應的楔形窗.圖2a中，空間域內(nèi)的曲波函數(shù)由不同顏色的橢圓標記.其中，綠色代表第2尺度內(nèi)的曲波函數(shù)，藍色代表第3尺度的曲波函數(shù)，而紅色和黃色代表第4尺度內(nèi)不同方向的曲波函數(shù).圖2b中展示了波數(shù)域中對應的楔形窗.由圖可以看出，隨著尺度參數(shù)j增加，對應的曲波函數(shù)變小，其各向異性增強.值得注意的是，由于二維傅里葉變換存在90°的旋轉(zhuǎn)關(guān)系，空間域內(nèi)的曲波函數(shù)與波數(shù)域內(nèi)的窗函數(shù)在方向上互相正交.

圖3展示曲波函數(shù)在兩個互相垂直的方向上的分布特征.沿主軸方向，其波形從中心向兩側(cè)緩慢衰減，而在垂直于主軸方向，其波形從中心向兩側(cè)呈現(xiàn)快速震蕩衰減.此外，曲波函數(shù)波形的長度與寬度還滿足如下各向異性尺度關(guān)系：

(15)

曲波變換的各向異性特征使得其能更好地表示不同方向上的特征.這種特征非常適合于航空電磁數(shù)據(jù)中測線誤差調(diào)平.

1.1.3 基于Wrapping算法的離散曲波變換

Candès等(2006)提出了用于快速離散曲波變換的卷繞(Wrapping)和USFFT算法.本文中我們采用卷繞算法.(13)式中我們已經(jīng)定義了離散曲波函數(shù)，我們采取與(12)式類似的形式在波數(shù)域中計算曲波系數(shù)，即：

(16)

式(16)經(jīng)離散化可得：

(17)

必須指出，由于楔形區(qū)域無法直接進行傅里葉變換，利用式(17)進行離散曲波變換存在難度.為了解決此問題，我們在圖4給出“Wrapping”算法的實施步驟.圖4中高度為L1j、寬度為L2j黑色平行四邊形為波數(shù)域中楔形的支撐空間，該空間內(nèi)包含的信息由灰色橢圓形區(qū)域表示.為進行二維傅里葉變換，我們定義一個圍繞坐標原點、邊長分別為L1j和L2j的正四邊形，并將包含在黑色平行四邊形中的數(shù)據(jù)信息Wrapping到圍繞原點的正四邊形中，從而可進行二維傅里葉變換.

圖2 空間域與波數(shù)域中的曲波(a) 空間域中四個不同的曲波函數(shù)； (b) 波數(shù)域中對應的楔形支撐空間(參考Herrmann and Hennenfent, 2008).Fig.2 Curvelets in spatial and wavenumber domains(a) Four different curvelets in spatial domain; (b) Wedge-shaped support in wavenumber domain(Refer to Herrmann and Hennenfent, 2008).

圖3 不同方向曲波函數(shù)波形(a) 沿主軸方向； (b) 垂直于主軸方向.Fig.3 Curvelet waveforms in varied directions(a) Along the major axis; (b) Perpendicular to the minor axis.

Wrapping算法可歸納為如下4個步驟：

圖4 曲波變換中的Wrapping算法(參考Candès et al.，2006)Fig.4 Wrapping algorithm in curvelet transform(Refer to Candès et al.，2006)

1.2 基于閾值法的電磁數(shù)據(jù)調(diào)平

本節(jié)我們首先解釋曲波變換可用于航空電磁數(shù)據(jù)調(diào)平的原因，然后討論具體調(diào)平方法.我們這里不考慮測區(qū)調(diào)平誤差，而將注意力集中于測線調(diào)平誤差和短波長調(diào)平誤差.這是由于前者可采用較為簡單的加/減背景場實現(xiàn)，而測線調(diào)平相對來說要困難的多.然而，曲波變換特別適合于對測線誤差進行調(diào)平.事實上，如前所述，測線調(diào)平誤差主要是由飛行方向改變導致的系統(tǒng)狀態(tài)變化引起的，它們在垂直于測線的方向上具有相對穩(wěn)定的空間尺度.另外，調(diào)平誤差沿測線方向比較穩(wěn)定(或者伴隨緩慢的非線性漂移).所有這些導致該誤差具有沿測線極強的方向性.因此，我們可以對信號進行多尺度和多方向分析，進而在曲波域特定尺度和角度方向上找出與調(diào)平誤差相關(guān)的信息，進行有用信號和調(diào)平誤差的有效分離，實現(xiàn)調(diào)平誤差去除.

圖5展示了理論模型的重構(gòu)結(jié)果.圖5a給出理論數(shù)據(jù)，我們分別將振幅為最大異常的12.5%及振幅為0的調(diào)平誤差交替加入到圖5a相鄰測線的數(shù)據(jù)中，得到了圖5b中包含調(diào)平誤差的數(shù)據(jù).考慮到理論模型共有60×60單元，我們在曲波域設置了4個尺度，分別在粗尺度上設置了小波系數(shù)，在第2尺度上設置了16個角度方向，在第3和第4尺度上設置了32個角度方向.進而，我們提取不同尺度和不同角度方向中的曲波系數(shù)，將信號重構(gòu)回空間域，得到特定尺度和角度方向中的信號成分.從圖5給出的結(jié)果可以看出，隨著尺度參數(shù)j的增大，空間尺度逐漸變得精細.從圖5f可以看出，調(diào)平誤差集中在第4尺度，而從圖5g—k可以看出，調(diào)平誤差具有明顯的方向性.對比圖5h和圖5k可以看出，第4與第20個角度具有相同的方向，這是由于它們在波數(shù)域中的支撐楔處于同一方向上的兩個頂端.然而，由于這兩個支撐楔存在位置、形態(tài)和對稱性差異，它們包含的信息完全不同.最后，考慮到調(diào)平誤差僅存在于第4尺度的第4、5、20、21角度方向，我們僅對這些尺度和方向進行數(shù)據(jù)調(diào)平，得到的結(jié)果由圖5l給出.對比圖5l和圖5a可以看出，利用本文算法數(shù)據(jù)中的調(diào)平誤差被去除，信號得到有效恢復.

圖5 理論模型及不同尺度、不同角度方向的重構(gòu)結(jié)果(a) 原始數(shù)據(jù)； (b) 加入調(diào)平誤差的數(shù)據(jù)； (c) 粗尺度重構(gòu)結(jié)果； (d) 第2尺度重構(gòu)結(jié)果； (e) 第3尺度重構(gòu)結(jié)果； (f) 第4尺度重構(gòu)結(jié)果； (g) 第4尺度第1角度重構(gòu)結(jié)果； (h) 第4尺度第4角度重構(gòu)結(jié)果； (i) 第4尺度第5角度重構(gòu)結(jié)果； (j) 第4尺度第12角度重構(gòu)結(jié)果； (k) 第4尺度第20角度重構(gòu)結(jié)果； (l) 曲波變換調(diào)平結(jié)果.Fig.5 Theoretical model and reconstruction results from curvelet transform at different scales and angles(a) Original data；(b) Data with leveling errors added; (c) Coarse-scale reconstruction; (d) 2nd scale reconstruction result; (e) 3rd scale reconstruction; (f) 4th scale reconstruction; (g) Reconstruction of 1st angle from 4th scale; (h) Reconstruction of 4th angle from 4th scale; (i) Reconstruction of 5th angle from 4th scale; (j) Reconstruction of 12th angle from 4th scale; (k) Reconstruction of 20th angle from 4th scale; (l) Leveling by curvelet transform.

在確定了調(diào)平誤差聚集的尺度和角度后，我們可以利用閾值法消除這些誤差.本文中我們將硬閾值法、軟閾值法兩種方法應用于航空電磁數(shù)據(jù)調(diào)平.

硬閾值定義為

(18)

其中，x為曲波系數(shù)，T為閾值.硬閾值的主要思想是通過對低于閾值的系數(shù)進行置零來消除干擾.硬閾值法雖然易于實現(xiàn)，但其缺點是可能受局部畸變的干擾.為了使處理后的曲波系數(shù)更具連續(xù)性，我們定義如下形式的軟閾值：

(19)

由(19)式可知，軟閾值對所有系數(shù)都進行了修改.與硬閾值法相比，軟閾值法可以得到更加平滑的結(jié)果，但可能給數(shù)據(jù)造成較大的失真.

本文將結(jié)合兩種閾值方法的優(yōu)點.為此，我們采取了2個步驟來消除調(diào)平誤差：

(1)我們通過硬閾值法x=hard(x,Ta),x∈Xja,la消除典型的條帶狀誤差；

(2)我們通過軟閾值法x=soft(x,Tb)，x∈X去除上一步驟的殘余誤差及其他短波長誤差.其中，集合X由所有的曲波系數(shù)組成，Ta、Tb是兩個步驟中相應的閾值.

上述處理程序可以給出更進一步的說明.事實上，第一步中硬閾值法及相應閾值Ta主要用于調(diào)平誤差集中的尺度和角度方向ja,la.由于調(diào)平誤差的幅值未知，我們在這一步驟需要設置較大的閾值Ta，并采用硬閾值法以避免信號失真.當ja,la中調(diào)平誤差幅值相對有用信號很大時，我們可以將所有ja,la內(nèi)的曲波系數(shù)置零.然而，調(diào)平誤差是非常復雜的，可能存在方向性不明顯的局部擾動，或者存在第一步調(diào)平后的非線性殘余.在這種情況下，我們可以對所有的尺度和角度應用軟閾值法并設置閾值Tb，進一步去除前一步的調(diào)平殘差和短波長誤差.

在實際的數(shù)據(jù)處理中，首先通過曲波變換將數(shù)據(jù)轉(zhuǎn)換至曲波域，然后根據(jù)調(diào)平誤差的特征，選擇合適的方向閾值或全局閾值進行處理.由于調(diào)平誤差的復雜性，我們通過特定尺度和特定角度數(shù)據(jù)的重構(gòu)判斷所選閾值的合理性.若重構(gòu)結(jié)果中方向性明顯的誤差消失，說明所選閾值合理可行.最后，我們將數(shù)據(jù)變換到空間域，如此便可以消除數(shù)據(jù)中的調(diào)平誤差.

1.3 航空電磁數(shù)據(jù)調(diào)平步驟

航空電磁信號對飛行高度很敏感，因此在計算視電阻率或進行反演時，飛行高度往往作為變量.在對航空電磁數(shù)據(jù)進行調(diào)平時，我們需要首先校正飛行高度的影響.為此，我們采用Huang(2008)和Beiki等(2010)提出的方法，即利用收發(fā)距與飛行高度之比的立方對測量信號的實虛分量進行歸一化處理：

(20)

(21)

其中，I和Q分別為電磁信號實虛分量，m和n為受飛行高度影響較小的歸一化響應.Huang(2008)指出，當?shù)叵陆橘|(zhì)相對良導時，(20)和(21)式很準確；然而，對于高阻介質(zhì)，為獲得準確結(jié)果，式(21)中的(h/s)的階數(shù)應該稍微降低.Beiki等(2010)利用這一近似將航空電磁信號校正到某一參考高度進行調(diào)平.本文中我們首先應用(20)和(21)式將數(shù)據(jù)進行歸一化，進而利用上述基于曲波變換的調(diào)平方法對歸一化的響應m和n進行調(diào)平，然后我們將數(shù)據(jù)轉(zhuǎn)換回到實虛分量.最后,通過已調(diào)平的數(shù)據(jù)計算視電阻率.如果視電阻率中沒有明顯的條帶狀誤差，則調(diào)平完成.由于曲波變換的多尺度和多方向特點，基于曲波變換的調(diào)平方法只對數(shù)據(jù)的某些細節(jié)尺度和角度進行局部改變，故對數(shù)據(jù)的原始相位關(guān)系影響很小，在去除調(diào)平誤差時不容易引起數(shù)據(jù)畸變.

2 合成數(shù)據(jù)調(diào)平

下面我們首先測試本方法在合成數(shù)據(jù)調(diào)平上的有效性.考慮到調(diào)平誤差可存在于各種航空數(shù)據(jù)中，本節(jié)僅給出不具有物理意義的合成模型.對于所有模型，我們在曲波域設置了5個尺度，角度數(shù)分別為第2尺度16個、第3、4尺度32個、第5尺度64個.圖6a—c給出了原始數(shù)據(jù)、加入調(diào)平誤差的數(shù)據(jù)以及調(diào)平誤差，圖6d給出了采用我們的曲波變換方法進行調(diào)平的結(jié)果，為對比起見，圖6e給出了Huang(2008)提出的基于測線相關(guān)性調(diào)平方法的結(jié)果.對比圖6d和圖6e，我們發(fā)現(xiàn)基于測線相關(guān)性的調(diào)平方法更容易造成數(shù)據(jù)失真，這主要是由于該方法難以利用多項式擬合復雜的非線性調(diào)平誤差.相比之下，曲波變換由于出色的二維信號表達能力，很好地還原了原始數(shù)據(jù).從圖6f和圖6g展示的原始數(shù)據(jù)與調(diào)平后數(shù)據(jù)差中，我們也可以清楚地看到，本文方法的殘差比其他方法的殘差小.在點位(4000 m, 3000 m)附近，有一個延伸方向與測線相同的異常體.從理論上講，這類異常體對于調(diào)平而言是較為困難的.然而，由于曲波變換的多尺度和多方向特性，我們成功地將調(diào)平誤差從有效信號中分離出來，很好地恢復了該異常特征.

圖6 (a) 原始數(shù)據(jù)； (b) 加入調(diào)平誤差的數(shù)據(jù)； (c) 調(diào)平誤差； (d) 曲波變換調(diào)平結(jié)果； (e) 基于測線相關(guān)性的調(diào)平結(jié)果(Huang，2008)； (f) 圖(a)與(d)的差值； (g) 圖(a)與(e)的差值Fig.6 (a) Original data； (b) Data with leveling errors； (c) Leveling errors； (d) Leveling of curvelet transform; (e) leveling of line-to-line correlation by Huang (2008)； (f) Difference between (a) and (d)； (g) Difference between (a) and (e)

航空電磁中經(jīng)常會出現(xiàn)測線方向與大地坐標系不一致的情況.此時，傳統(tǒng)的調(diào)平方法(例如基于二維濾波的方法)需要對坐標系進行旋轉(zhuǎn).由于曲線變換的方向性，我們可以省去這些旋轉(zhuǎn)，僅通過選擇不同的角度參數(shù)很容易地解決問題.這意味著基于曲波變換的調(diào)平方法具有更大的靈活性.圖7給出一個測線與大地坐標軸斜交的模型和對應的曲波變換調(diào)平結(jié)果.由圖可以看出，調(diào)平誤差被很好地消除.

圖7 (a) 原始數(shù)據(jù)； (b) 加入調(diào)平誤差的數(shù)據(jù)； (c) 曲波變換調(diào)平結(jié)果Fig.7 (a) Original data； (b) Data with leveling errors； (c) Leveling of curvelet transform

航空電磁法中另一個常見的問題是測區(qū)形狀不規(guī)則.在這種情況下，數(shù)據(jù)調(diào)平不僅費時，而且調(diào)平結(jié)果也會受到影響.例如，基于測線相關(guān)性的調(diào)平方法必須對擬合多項式進行外插，而基于二維濾波的調(diào)平方法在不規(guī)則的邊緣處會產(chǎn)生一定的邊緣效應.然而，基于曲波變換的調(diào)平方法由于具有多尺度性，可以有效地處理不規(guī)則測區(qū)獲取的數(shù)據(jù).事實上，在調(diào)平過程中，我們用一個常數(shù)來填充缺失的區(qū)域，然后進行調(diào)平.由于缺失區(qū)域具有一定的尺度，該區(qū)域的信息大多出現(xiàn)在曲波域的粗尺度上，對調(diào)平結(jié)果影響較小，因此在完成數(shù)據(jù)調(diào)平后我們僅需將該區(qū)域的數(shù)據(jù)剔除即可.圖8給出形狀不規(guī)則的合成模型的調(diào)平結(jié)果.由圖可以看出，盡管測區(qū)由于存在三角形缺失變得不規(guī)則，調(diào)平誤差仍然被很好地消除，數(shù)據(jù)得到有效地恢復.

圖8 (a) 原始數(shù)據(jù)； (b) 加入調(diào)平誤差的數(shù)據(jù)； (c) 曲波變換調(diào)平結(jié)果Fig.8 (a) Original data； (b) Data with leveling errors； (c) Leveling of curvelet transform

3 實測數(shù)據(jù)調(diào)平

為了進一步驗證本文方法的有效性，我們將其應用于愛爾蘭Waterford地區(qū)利用頻率域航空系統(tǒng)采集的電磁數(shù)據(jù).該項目由愛爾蘭地質(zhì)調(diào)查局(GSI)牽頭，由加拿大Sander地球物理公司與Unicorn礦業(yè)有限公司合作組織實施.圖9給出Waterford測區(qū)的位置及測線布設情況.其中，測線方向為N15°W，線距200 m，切割線方向為E15°N，線距2000 m.我們選擇圖9中紅色線框標注的面積約為11.5 km×13 km的區(qū)域進行測試.我們在曲波域中設置5個尺度，其中第2尺度包含32個角度,第3、4尺度包含64個角度，第5尺度包含128個角度.我們以頻率912 Hz為例對實測數(shù)據(jù)實虛部進行調(diào)平.

圖9 Waterford項目測區(qū)位置及測線布設(修改自Bates和Cho, 2016)Fig.9 Survey area and survey lines for the Waterford project(Modified from Bates and Cho, 2016)

圖10a和圖10d為實、虛部原始數(shù)據(jù).我們首先利用(20)、(21)式將數(shù)據(jù)轉(zhuǎn)換為歸一化響應m和n，然后利用本文方法對m和n進行調(diào)平，再將數(shù)據(jù)轉(zhuǎn)換到調(diào)平后的實虛分量.圖10b和圖10e為調(diào)平后的數(shù)據(jù).為了便于比較，我們在圖10c和圖10f中還展示了GSI采用的基于微分多項式擬合的調(diào)平結(jié)果(Beiki et al.，2010).由對比結(jié)果可見，對于虛部數(shù)據(jù)，基于曲波變換和基于GSI的多項式擬合方法均獲得很好的結(jié)果；對于實部數(shù)據(jù)，基于曲波變換的調(diào)平方法有效去除了調(diào)平誤差，未對數(shù)據(jù)產(chǎn)生明顯的畸變，然而，GSI利用的多項式擬合方法在消除調(diào)平誤差的同時，也去除了同方向上的異常(圖10c).

圖10 愛爾蘭Waterford測區(qū)實部(左側(cè))與虛部(右側(cè))數(shù)據(jù)調(diào)平結(jié)果(a)與(d)為原始數(shù)據(jù)； (b)與(e)為曲波變換調(diào)平結(jié)果； (c)與(f)為GSI的調(diào)平結(jié)果.Fig.10 Leveling results of in-phase (left column) and quadrature component (right column) for the data from Waterford project, Ireland(a) and (d) Original data； (b) and (e) Leveling by curvelet transform； (c) and (f) Leveling by GSI.

由于本文方法對實虛分量的處理是獨立進行的，因此對于調(diào)平結(jié)果來說，檢驗其視電阻率和相位變化至關(guān)重要.圖11給出圖10中調(diào)平前后數(shù)據(jù)計算的視電阻率和相位.由圖可以看出，由于實虛部數(shù)據(jù)中含有調(diào)平誤差，故利用原始數(shù)據(jù)計算視電阻率及相位時，誤差被引入其中，導致計算結(jié)果中存在明顯的條帶狀誤差.經(jīng)曲波變換調(diào)平后，視電阻率數(shù)據(jù)中的條帶狀誤差明顯減少，相鄰測線間電性趨于相容、合理.作為對比，由GSI給出的調(diào)平結(jié)果計算得到的視電阻率在中間相對高阻區(qū)域仍存在明顯的條帶狀誤差，且造成了一定程度的數(shù)據(jù)畸變，特別是相位數(shù)據(jù)畸變更為嚴重.這主要是由于該方法在分別對實虛部進行調(diào)平時對原始數(shù)據(jù)畸變過多造成的.從圖11d、e、i、j可以得出相似的結(jié)論，即本文的調(diào)平方法對原始數(shù)據(jù)的視電阻率和相位畸變較小.

圖11 愛爾蘭Waterford測區(qū)調(diào)平前后數(shù)據(jù)計算得到的視電阻率(左側(cè))與相位(右側(cè))(a)與(f)由原始數(shù)據(jù)得到； (b)與(g)由曲波變換調(diào)平結(jié)果得到； (c)與(h)由GSI提供的調(diào)平結(jié)果得到； (d)與(i)分別為(a)與(b)、(f)與(g)的差； (e)與(j)分別為(a)與(c)、(f)與(h)的差.Fig.11 Apparent resistivity (left column) and phase (right column) before and after leveling for the data from Waterford project, Ireland(a) and (f) From original data; (b) and (g) From results by curvelet transform； (c) and (h) From results of GSI; (d) and (i) Differences between (a) and (b), (f) and (g)；(e) and (j) Differences between (a) and (c), (f) and (h).

4 結(jié)論

本文我們基于曲波變換成功地開發(fā)了一種航空電磁數(shù)據(jù)調(diào)平方法.該方法較之于傳統(tǒng)數(shù)據(jù)調(diào)平方法具有如下優(yōu)勢：

(1)利用曲波變換多尺度和多方向的優(yōu)點，可以方便地提取與飛行方向相關(guān)的調(diào)平誤差，同時有效保留了原始信號.

(2)對于實虛部數(shù)據(jù)的調(diào)平僅在曲波域的部分尺度和方向上進行，保留了實虛部數(shù)據(jù)中的主要成分，因此由其計算的視電阻率在調(diào)平誤差被極大削弱的同時，原始相位關(guān)系得到很好保證.

(3)對于傾斜測線及不規(guī)則測區(qū)，本文調(diào)平方法無需坐標旋轉(zhuǎn)和測區(qū)分割，同時避免了傳統(tǒng)方法產(chǎn)生的邊緣效應，具有很好的實用性.

理論和實測數(shù)據(jù)的實驗結(jié)果均表明，本文基于曲波變換的調(diào)平方法比傳統(tǒng)調(diào)平方法取得了更好的效果.由于調(diào)平誤差在航空地球物理數(shù)據(jù)中廣泛存在，而本文提出的方法實質(zhì)上與數(shù)據(jù)類型無關(guān)，我們可以期待將該方法進一步推廣應用于其他航空地球物理數(shù)據(jù)處理中.

致謝感謝愛爾蘭地質(zhì)調(diào)查局(GSI)與加拿大Sander地球物理有限公司允許本文使用他們的Waterford數(shù)據(jù)并提供技術(shù)支持.