孫 勇，周 帥

（吉林大學地球探測科學與技術(shù)學院,長春130026）

0 引言

航空地球物理測量的特點之一是測點間距遠遠小于測線間距，出于避免插值失真的考慮，網(wǎng)格化間距通常被定為測線間距的1/4左右。因此，測線與測線之間彼此的約束小，沿測線的低頻噪聲就會以條帶狀的形式顯示在網(wǎng)格化后的圖上。調(diào)平是用來去除數(shù)據(jù)中的低頻干擾的常用方法，如移除線性或非線性的漂移、降低機械噪聲的干擾等，最終使測量數(shù)據(jù)接近真實的場值。調(diào)平已經(jīng)成為航空地球物理測量數(shù)據(jù)預(yù)處理工藝中的標準步驟之一，可以為后續(xù)的解釋工作提供高質(zhì)量的地球物理數(shù)據(jù)。

在具體的操作方法上，已有大量學者做出過研究。已有的調(diào)平方法可分為統(tǒng)計調(diào)平（切割線調(diào)平）與經(jīng)驗調(diào)平（微調(diào)平）。在統(tǒng)計調(diào)平這一范疇內(nèi)，F(xiàn)oster等[1]用一組正交基來擬合測線與切割線上系統(tǒng)誤差，假設(shè)系統(tǒng)誤差緩慢變化，用最小二乘擬合原始數(shù)據(jù)與用正交基擬合的數(shù)據(jù)，以此來分離出這種意義下的 “真實”解。Harold等[2]在處理航空磁測數(shù)據(jù)時，通過最小化交叉點上的值來去除日變場的干擾，并且事先不需要地面基站的測量值。Bandy等[3]在Foster的研究基礎(chǔ)之上推導了當測線與切割線在水平面上的投影為直線時的調(diào)平方法，簡化了計算步驟，其研究結(jié)果可被視為Forster提出的最小二乘統(tǒng)計調(diào)平的一個特例。Huang等[4-5]在研究航空電磁場數(shù)據(jù)與航空磁測數(shù)據(jù)的調(diào)平時提出了一種基于每條測線逐步調(diào)平的方法。以誤差干擾小的切割線為基準，將切割線與鄰近的測線結(jié)合起來，用參數(shù)方程擬合了測線上的需要調(diào)平的誤差，再利用最小二乘準則求取了切割線與測線在二范數(shù)意義上的最小，求得了方程中的參數(shù)。Gleave等[6]在研究航空重力數(shù)據(jù)的調(diào)平時發(fā)現(xiàn)，在用濾波切除噪聲干擾時，容易將與噪聲等波長的信號也切除掉，因此其建議將區(qū)域場先進行分離再做調(diào)平處理。統(tǒng)計調(diào)平的優(yōu)點是可以對大范圍的測線進行某種規(guī)則下的統(tǒng)一處理，缺點是在統(tǒng)計的意義下難以準確區(qū)分真實信號與噪聲信號。Davydenko等[7]在處理航空電磁數(shù)據(jù)的調(diào)平時，引入了一種主成分分析法，來設(shè)計方向濾波器。這種方法假設(shè)切割線中的噪聲差異大于測線，從而對方向濾波器的參數(shù)做出了設(shè)定。Zhang等[8]在Davydenko的主成分分析法上做出了改進，進一步降低了對切割線的依賴，使方向濾波器在沒有切割線的地區(qū)同樣可以應(yīng)用。微調(diào)平專門用來處理統(tǒng)計調(diào)平所不能妥善處理的測線，微調(diào)平是一種經(jīng)驗調(diào)平，針對個別測線可用特定參數(shù)的濾波器進行處理，濾波器的參數(shù)則根據(jù)處理人員對該地區(qū)的了解而定。Nuady等[9]設(shè)計了一種非線性濾波器，該濾波器后被我國國土資源航遙中心進行了引進及改良。Minty[10]根據(jù)航磁數(shù)據(jù)調(diào)平誤差與飛行方向的相關(guān)性設(shè)計出了一種方向濾波器。對統(tǒng)計調(diào)平與微調(diào)平的詳細技術(shù)描述，也可以參見商用軟件Oasis Montaji Geosoft中關(guān)于調(diào)平的幫助文檔[11]。

本文在上述研究的基礎(chǔ)上，針對航空重力梯度數(shù)據(jù)，提出了一種組合調(diào)平方法。重力梯度數(shù)據(jù)與航磁數(shù)據(jù)或航空電磁數(shù)據(jù)相比，受到的交叉測線位置不同、高度不同，飛行方向與載體質(zhì)量干擾更大，相同之處在于同樣包括儀器的漂移干擾。用商用軟件Geosoft中的調(diào)平模塊直接對數(shù)據(jù)進行處理，往往得不到良好的效果。因此，本文在統(tǒng)計調(diào)平的基礎(chǔ)之上增加了3個處理步驟，分別是原始數(shù)據(jù)的去趨勢、融合過密的測線與增加方向濾波器，3個處理步驟均在統(tǒng)計調(diào)平之前進行。經(jīng)過實測數(shù)據(jù)驗證，這種組合調(diào)平的方法相比單純進行低通濾波或只進行統(tǒng)計調(diào)平可以取得更佳的效果。

1 組合調(diào)平的原理

傳統(tǒng)的濾波手段會將噪聲及與噪聲波長相等的信號一起濾除，而組合調(diào)平利用沿測線的濾波器保留高頻信息，再用網(wǎng)格化的數(shù)據(jù)進行低通濾波，以保留低頻信息。最后，重新采樣到原測線，與高頻數(shù)據(jù)合并以得到最終結(jié)果。數(shù)據(jù)流程圖如圖1所示。

圖1 組合調(diào)平的數(shù)據(jù)流程圖Fig.1 Data flow diagram of combined leveling

圖1左側(cè)分支主要的步驟為去趨勢，去趨勢可以在空間域與頻率域進行。在空間域進行計算時，除了早期的圖像法與多項式擬合法之外，多數(shù)要進行褶積運算，計算復(fù)雜度高，時間成本高。因此，出于方便計算的考慮，大部分算法均在頻率域中進行。在頻率域中，常用的背景場去除方法包括了向上延拓法、維納濾波、小波變換等。針對不同測區(qū)數(shù)據(jù)的特點，可以采用不同的方法。在本節(jié)的調(diào)平處理中，由于不需要準確地分離出區(qū)域場值，分離處理是為了避免在處理高頻的噪聲干擾時損失一部分信號。因此，可選用較為方便的空間域多項式擬合方法來分離趨勢場。

以測線l為基本單位，用n階多項式來擬合的區(qū)域場為Test。由于越高階的項在趨勢場中的占比越小，因此權(quán)因子b可以構(gòu)造為[0，1]區(qū)間內(nèi)的形式

式（1）中，x、y代表測點坐標，每一個測點都被單獨賦予了一個權(quán)值b。

單條測線的重力梯度趨勢場擬合值Test可表達為

式（2）中，B為用坐標信息求出的系數(shù)矩陣，l為未知數(shù)，n為多項式的階數(shù)，m為單條測線上測點的個數(shù)。設(shè)趨勢場的真值為Treal，在二范數(shù)最小的前提下求取最小二乘解，可以得到

將L帶回到式（2）中，可以得到趨勢場的最小二乘計算結(jié)果。

圖1所示的右側(cè)分支包括網(wǎng)格化、方向濾波與統(tǒng)計調(diào)平，其算法流程圖如圖2所示。

圖2的核心步驟是選擇合適的滑動窗口對相鄰的測線進行平均處理，窗口大小由窗口內(nèi)測線的條數(shù)來定義。相鄰測線包含的噪聲相類似，因此，可以取經(jīng)過平滑處理的相鄰測線與原始測線差值的平均值作為噪聲的近似值，再用原始數(shù)據(jù)減去這個近似值，以作為輸出。在圖2中，統(tǒng)計調(diào)平模塊的編制參考了胥值禮等[12]對航磁數(shù)據(jù)調(diào)平的研究成果，其核心思路是使切割線與測線的值相接近。其步驟是首先調(diào)切割線，其次用調(diào)平后的切割線調(diào)平測線。設(shè)Te為修正前切割線的重力梯度值，Tl為修正前測線的重力梯度值，Te′為修正后切割線的重力梯度值，Tl′為修正后測線的重力梯度值。切割線的調(diào)整方法為

圖2 圖1右側(cè)分支的算法流程圖Fig.2 Algorithm flow chart of the right branch in Fig.1

在式（3）中，i為切割線的編號，j為測線的編號，Tcei，k表示第i條切割線上與測線相交的第k個點的重力梯度值，p1為第i條切割線上與測線的交點的個數(shù)，Tclj，k為k點處測線上的重力梯度值。因此，式（3）中ΔTei的物理意義為第i條切割線在與測線相交的各點處差值的算數(shù)平均值。測線的調(diào)整方法為

式（4）中，ΔTlj為第j條測線與切割線相交處差值的算數(shù)平均值，Tce′i，k為在第k個交點上調(diào)整后切割線上的重力梯度值，Tclj，k為第k個交點上測線的重力梯度原值。經(jīng)過式（3）與式（4）的處理，可以拉近切割線與測線的差異，最終完成統(tǒng)計調(diào)平處理。

2 實測數(shù)據(jù)驗證

以2003年在加拿大西北部的Diavik鉆石礦區(qū)進行的航空全張量重力梯度測量為例，來驗證調(diào)平方法的有效性。圖3（a）展示了以金伯利巖為主的礦區(qū)，圖3（b）展示了測線軌跡，切割線共有13條。

圖3 礦區(qū)位置圖Fig.3 Map of mine location

圖 4（a）～圖 4（f）展示了調(diào)平之前I1～I3、C1～C3數(shù)據(jù)的形態(tài)。可以看出，在未經(jīng)調(diào)平處理時，數(shù)據(jù)包含了很大的與測線相關(guān)的方向噪聲。為了方便說明調(diào)平的效果，在圖4中，用紅色線框圈出了相應(yīng)的計算區(qū)域。

在圖4中，I1～I3、C1～C3的表達式分別為

圖4 各分量在調(diào)平之前的形態(tài)Fig.4 Shape of each component before leveling

在式（5）中，u、v、w為傘形坐標系下三個坐標軸的方向。為與本節(jié)提出的調(diào)平方法所得到的結(jié)果作出對比，首先需對數(shù)據(jù)進行常規(guī)的低通濾波處理。以圖 5（a）所示的分量C1為例，經(jīng)過圖5（b）所示的截至波數(shù)為 1.84（單位 1/km）的空間域低通濾波，可得到如圖5（c）所示的濾波結(jié)果，其與原數(shù)據(jù)的差如圖5（d）所示。可以看到，大尺度的低通濾波可以去除相當多的干擾，濾波前后信號的最大值從165.8E降低到了146.9E。但對過濾之后的圖像進行分析，可以看到濾波過程也濾掉了一部分有效的信號，即濾掉的部分不完全是條帶狀的干擾。

下面用本文提出的組合調(diào)平方法對C1進行處理，趨勢場的分離結(jié)果如圖6所示。

從圖6可以看出，經(jīng)過去趨勢的處理，大部分沿測線方向的噪聲都已經(jīng)被去除，但在測區(qū)的北部及南北向的切割線上，仍能看到沿測線方向的強干擾的噪聲。

圖5 對C1進行低通濾波的結(jié)果與相應(yīng)的濾波器參數(shù)設(shè)置Fig.5 Results of low-pass filtering for C1and the corresponding filter parameter setting

圖6 C1去趨勢前后的對比Fig.6 Comparison of C1before and after detrending

用此方法進行處理的關(guān)鍵之處在于相鄰測線的框定范圍，即選定合適的窗口大小。在垂直于測線方向的窗口分別為15條測線、10條測線、5條測線這3種不同的條件下，對C1的處理結(jié)果如圖7所示。

圖7 不同窗口大小對C1進行滑動的平均結(jié)果對比Fig.7 Comparison of sliding average results for C1with different window sizes

由圖7（a）～圖 7（c）可知，當窗口為 5 時，可以取得最佳的效果。此時，仍能觀測到沿測線方向的噪聲，因此需要做進一步的統(tǒng)計調(diào)平處理。對I1～I3、C1～C3做出相同大小的測線數(shù)據(jù)集與切割線數(shù)據(jù)集，之后統(tǒng)計交叉點上的信息，用式（3）與式（4）進行計算，得到的處理結(jié)果如圖8所示。

圖8 經(jīng)過統(tǒng)計調(diào)平的各道數(shù)據(jù)集Fig.8 After statistical adjustment of each channel data set

圖9 濾波與組合調(diào)平的結(jié)果對比Fig.9 Comparison of the results of filtering and combined leveling

濾波與組合調(diào)平的結(jié)果對比如圖9所示，組合調(diào)平的結(jié)果上異常形態(tài)更加完整，沿測線方向的噪聲也得到了更好的壓制。組合調(diào)平的結(jié)果異常幅值最大值為63.7E，這說明簡單進行低通濾波處理的結(jié)果不僅沒有保留完整的異常形態(tài)，而且噪聲也沒有被完全去除。

3 結(jié)論

航空地球物理測量的高動態(tài)性使其無法避免被與航向有關(guān)的噪聲干擾，如何在去除干擾的同時保留盡可能多的異常細節(jié)是處理過程中的難點，特別是對于高分辨率的航空重力梯度測量數(shù)據(jù)而言。與低通濾波處理或傳統(tǒng)的切割線調(diào)平處理不同，本文提出了一種組合調(diào)平方法。經(jīng)過實測數(shù)據(jù)驗證，該方式可以明顯提高原始數(shù)據(jù)的質(zhì)量。從圖9可以看出，在經(jīng)過處理后，仍有少量非地質(zhì)體信號的成分存在。因此，建議在組合調(diào)平后進行逐條測線的微調(diào)平處理，以進一步提高數(shù)據(jù)的質(zhì)量。