孫 勇，周 帥

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

0 引言

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

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

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

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

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

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

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

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

式（1）中，x、y代表測點(diǎn)坐標(biāo)，每一個(gè)測點(diǎn)都被單獨(dú)賦予了一個(gè)權(quán)值b。

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

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

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

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

圖2的核心步驟是選擇合適的滑動(dòng)窗口對(duì)相鄰的測線進(jìn)行平均處理，窗口大小由窗口內(nèi)測線的條數(shù)來定義。相鄰測線包含的噪聲相類似，因此，可以取經(jīng)過平滑處理的相鄰測線與原始測線差值的平均值作為噪聲的近似值，再用原始數(shù)據(jù)減去這個(gè)近似值，以作為輸出。在圖2中，統(tǒng)計(jì)調(diào)平模塊的編制參考了胥值禮等[12]對(duì)航磁數(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為切割線的編號(hào)，j為測線的編號(hào)，Tcei，k表示第i條切割線上與測線相交的第k個(gè)點(diǎn)的重力梯度值，p1為第i條切割線上與測線的交點(diǎn)的個(gè)數(shù)，Tclj，k為k點(diǎn)處測線上的重力梯度值。因此，式（3）中ΔTei的物理意義為第i條切割線在與測線相交的各點(diǎn)處差值的算數(shù)平均值。測線的調(diào)整方法為

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

2 實(shí)測數(shù)據(jù)驗(yàn)證

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

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

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

在圖4中，I1～I(xiàn)3、C1～C3的表達(dá)式分別為

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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