趙德軍，孫中苗，趙東明

（1.信息工程大學(xué)地理空間信息學(xué)院，鄭州 450001；2.西安測(cè)繪總站，西安 710054；3.地理信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，西安 710054；4.西安測(cè)繪研究所，西安 710054）

相對(duì)于航空標(biāo)量重力測(cè)量和矢量重力測(cè)量而言，航空重力梯度測(cè)量有兩大顯著優(yōu)勢(shì)[1-3]：1）最大的優(yōu)勢(shì)在于其不受載體運(yùn)動(dòng)和厄特弗斯效應(yīng)的影響；2）重力梯度以重力場(chǎng)的曲率描述重力場(chǎng)區(qū)域結(jié)構(gòu)，包含地球物理學(xué)和大地測(cè)量中急需的局部重力場(chǎng)幾何信息（水準(zhǔn)面的平均曲率、鉛垂線的曲率等），因而重力梯度測(cè)量更能反映重力場(chǎng)的精細(xì)結(jié)構(gòu)，敏感重力場(chǎng)的短波變化。

航空重力梯度測(cè)量能獲取高分辨率和高頻譜的重力梯度信號(hào)，能更好地描述小的異常特征。因此，在航空重力梯度數(shù)據(jù)處理時(shí)，既要去除高頻噪聲，又要保留高頻的有效重力梯度信號(hào)，是一個(gè)巨大的挑戰(zhàn)。在生產(chǎn)作業(yè)中，像巴特沃斯這樣的經(jīng)典低通數(shù)字濾波器仍被廣泛采用，濾波器的截止波長(zhǎng)一般根據(jù)數(shù)據(jù)的采樣率和已知地質(zhì)構(gòu)造的特征信息來(lái)估計(jì)，并需要經(jīng)過(guò)多次試驗(yàn)得到。低通濾波器存在選擇最優(yōu)參數(shù)的難題，以及方法本身無(wú)法針對(duì)具體數(shù)據(jù)進(jìn)行參數(shù)自適應(yīng)選擇的問(wèn)題[1-3]。

克里金插值是地球物理統(tǒng)計(jì)學(xué)中常用的空間數(shù)據(jù)插值估計(jì)方法，可以客觀地估計(jì)數(shù)據(jù)的噪聲水平和平滑度。在地球物理中，觀測(cè)到的信號(hào)既可以看成是區(qū)域分量與局部分量的疊加，又可看成是深部地質(zhì)體與淺層地質(zhì)體所引起的異常的綜合反映，同時(shí)還含有不可忽略的系統(tǒng)誤差。因子克里金法允許將區(qū)域化變量分解成多個(gè)不同尺度的分量（即因子），每個(gè)因子都有其對(duì)應(yīng)的變異函數(shù)，每個(gè)因子的變程代表地質(zhì)體的平均直徑大小。

因子克里金法具有多尺度分析等諸多功能和優(yōu)點(diǎn)，因而廣泛應(yīng)用于地球物理與地球化學(xué)勘探[4]，圖像處理與分析[5,6]，水土環(huán)境監(jiān)測(cè)[7,8]等領(lǐng)域。但在重力梯度降噪方面未見(jiàn)相關(guān)報(bào)道，因此本文首次將其應(yīng)用于航空重力梯度數(shù)據(jù)降噪處理，取得了顯著的效果。模擬試驗(yàn)表明，相對(duì)于傳統(tǒng)的低通濾波器，重力梯度各分量精度平均提升了36%，垂直梯度分量尤為顯著，精度提升了42%。

1 因子克里金降噪法

1.1 變異函數(shù)

克里金理論，就是用變異函數(shù)來(lái)描述區(qū)域化變量以解決有關(guān)地學(xué)問(wèn)題的理論，其核心是變異函數(shù)。區(qū)域化變量的結(jié)構(gòu)分析，就是通過(guò)區(qū)域化變量有限的空間觀測(cè)值來(lái)構(gòu)建相應(yīng)的理論變異函數(shù)模型，以表征該變量的主要結(jié)構(gòu)特征。

設(shè)()Z x是區(qū)域化變量，有了二階平穩(wěn)假設(shè)或本征假設(shè)，就可以給出變異函數(shù)的計(jì)算公式[5]：

1.2 因子克里金

因子克里金將區(qū)域化變量分解成多個(gè)尺度的因子，而每個(gè)因子都對(duì)應(yīng)一個(gè)變異函數(shù)。因子克里金濾波的思路：從實(shí)驗(yàn)變異函數(shù)中，將各因子對(duì)應(yīng)的變異函數(shù)提取出來(lái)，認(rèn)為小尺度變異函數(shù)對(duì)應(yīng)的因子是需要濾掉的噪聲，保留大尺度的變異函數(shù)對(duì)應(yīng)的因子作為有效信號(hào)，從而實(shí)現(xiàn)濾波的過(guò)程。區(qū)域化變量可以認(rèn)為是不同地質(zhì)體相互作用的疊加，因此將區(qū)域化變量分解為[5]：

根據(jù)估計(jì)方差最小原則，對(duì)每個(gè)因子()k Y x導(dǎo)出克里金方程組[5-7]：

從因子克里金方程組可看出，因子克里金除了需要知道原始數(shù)據(jù)的變異函數(shù)()hγ，還要知道每個(gè)因子的變異函數(shù)。利用分解后各因子的變異函數(shù)按式(5)解出權(quán)系數(shù)，然后用式(4)就可得到對(duì)原空間數(shù)據(jù)各個(gè)因子的估計(jì)值。

理論變異函數(shù)通常是非線性模型，因此在采用加權(quán)最小二乘擬合之前，需要將非線性模型轉(zhuǎn)換為線性模型。一般采用泰勒級(jí)數(shù)將非線性模型轉(zhuǎn)換為線性模型，也可以采用多項(xiàng)式逼近的方法將其轉(zhuǎn)換成線性模型。另外，除了從數(shù)學(xué)上將非線性模型線性化外，還可以利用一些數(shù)學(xué)軟件里的遺傳算法、模擬退火法等最優(yōu)方法來(lái)解決非線性回歸問(wèn)題[9]。

2 模擬數(shù)據(jù)分析

對(duì)于模擬實(shí)驗(yàn)，可以采用噪聲衰減因子β來(lái)估計(jì)噪聲衰減量。用X表示模型的真實(shí)梯度值，Xn表示增加噪聲后的梯度值，Xr為經(jīng)過(guò)濾波處理后的梯度值:

式中var 表示計(jì)算的方差。當(dāng)β為正時(shí)，說(shuō)明噪聲被濾除，β的大小表示去除噪聲的能力，當(dāng)β趨于1，表示Xr中幾乎不含噪聲成分。

利用下面的均方根誤差來(lái)表示濾波后數(shù)據(jù)的精度：

模擬一個(gè)地質(zhì)特征顯著的重力梯度場(chǎng)，假定在平面坐標(biāo)系中，分布8 個(gè)已知梯度的樣本點(diǎn)，其坐標(biāo)和值如表1所示。采用套合變異函數(shù)模型，用普通克里金插值，內(nèi)插坐標(biāo)網(wǎng)格的梯度值。網(wǎng)格的x和y坐標(biāo)范圍均從0 m 到100 m，坐標(biāo)間隔1 m。為了使模擬的模型更有說(shuō)服力，變異函數(shù)采用1 個(gè)指數(shù)模型和1個(gè)球狀模型的套合結(jié)構(gòu)。指數(shù)模型的基臺(tái)值為20 E2，變程10 m，球狀模型的基臺(tái)值為4 E2，變程為60 m。模擬的模型梯度如圖1(a)，然后在模型梯度值的基礎(chǔ)上，添加方差為10 E2的高斯白噪聲，如圖1(b)。

根據(jù)圖1(b)統(tǒng)計(jì)出實(shí)驗(yàn)變異函數(shù)，然后采用1 個(gè)塊金模型，2 個(gè)球狀模型的套合變異函數(shù)來(lái)逼近實(shí)驗(yàn)變異函數(shù)，并提取出各因子的理論變異函數(shù)，其參數(shù)見(jiàn)表2。是塊金值的變異函數(shù)，通常認(rèn)為塊金值表示隨機(jī)噪聲，這里塊金值為8.24 E2，與假設(shè)的噪聲方差10 E2比較吻合。的變程較小，認(rèn)為是局部分量的變異函數(shù)，代表小尺度的地質(zhì)體；的變程較大，一般認(rèn)為是區(qū)域分量的變異函數(shù)，代表大尺度的地質(zhì)體。這樣就將圖1(b)分解成3 個(gè)不同尺度的分量。圖1(c)、圖1(d)和圖1(e)分別是對(duì)應(yīng)的大尺度區(qū)域分量，對(duì)應(yīng)的小尺度局部分量，以及對(duì)應(yīng)的隨機(jī)噪聲分量。

將大尺度分量和小尺度分量疊加，認(rèn)為是有效信號(hào)，如圖1(f)。對(duì)比圖1(f)和圖1(a)，可看出采用因子克里金濾波后的圖像仍然有部分噪聲殘留。同時(shí)，采用空域高斯濾波低通對(duì)圖1(b)降噪，高斯濾波器的方差設(shè)為10 E2，濾波后的梯度為圖1(g)。

表3 是不同方法濾波后的衰減因子和誤差表。對(duì)比圖1 和表3，可看出，高斯濾波雖然圖像比較平滑，但是其衰減因子比較低，只有0.76，誤差為1.56 E。因子克里金濾波后的圖像不太平滑，但衰減因子達(dá)到了0.93，誤差為0.91 E。所以可以說(shuō)，空域高斯濾波過(guò)度平滑，濾除噪聲的同時(shí)，還濾除了部分有用信號(hào)，而因子克里金濾除高頻噪聲的同時(shí)，還保留了高頻信號(hào)。

同時(shí)還采用普通克里金濾波對(duì)圖1(b)進(jìn)行降噪處理，濾波后的圖像接近圖1(f)。從表3 看出普通克里金濾波的衰減因子為0.90，誤差為0.97 E，濾波效果略低于因子克里金。

表1 模擬樣本點(diǎn)的坐標(biāo)和梯度值Tab.1 Coordinates and gradient values of simulated points

表2 擬合的套合變異函數(shù)的參數(shù)Tab.2 Parameters of fitted nested variogram

表3 噪聲衰減因子和誤差比較Tab.3 Noise attenuation factor and error comparison

圖1 因子克里金分析圖Fig.1 Factor Kriging Analysis diagram

本文還統(tǒng)計(jì)了濾波前后信號(hào)的徑向平均功率譜密度（圖2），從中看出因子克里金降低了高頻段的噪聲能量，但還殘留了部分噪聲，保留了梯度的真實(shí)成分。而高斯濾波降低高頻段的噪聲的同時(shí)，又濾掉了部分有用信號(hào)，過(guò)度平滑了信號(hào)。

圖2 濾波前后功率譜曲線對(duì)比圖Fig.2 Comparison of power spectrum curves before and after filtering

為了充分說(shuō)明因子克里金的降噪效果，本文還對(duì)重力梯度的其它5 個(gè)分量進(jìn)行了模擬實(shí)驗(yàn)。圖3 中從上到下分別是Txx、Txy、Txz、Tyy 和Tyz 共5 個(gè)分量的降噪效果圖。其中，圖3 第1 列為模擬的無(wú)噪聲梯度分量值，其參數(shù)與圖1(a)中垂直梯度分量Tzz 相同。第2 列為在第1 列模型梯度值的基礎(chǔ)上，添加方差為10 E2高斯白噪聲后的梯度分量。第3 列為采用因子克里金降噪后的梯度分量圖，從中看出對(duì)于不同的梯度分量，因子克里金都能有效地降低噪聲的干擾。

圖3 重力梯度分量降噪圖Fig.3 Noise reduction of gravity gradient component

表4 為根據(jù)式(9)統(tǒng)計(jì)的重力梯度分量降噪誤差表，可以看出普通克里金方法和因子克里金方法都能有效地降低重力梯度分量的噪聲，且因子克里金法略優(yōu)于普通克里金降噪法。因子克里金法相對(duì)于高斯空間濾波，各梯度分量精度均提升30%以上，平均精度提升了36%，尤其是垂直梯度分量精度提升了42%。

表4 重力梯度分量降噪誤差統(tǒng)計(jì)(單位：E)Tab.4 Error statistics for noise reduction of gravity gradient component(unit:E)

3 實(shí)測(cè)數(shù)據(jù)分析與降噪處理

數(shù)據(jù)來(lái)源于貝爾地理空間公司在加拿大圣喬治灣測(cè)量的重力垂直梯度Tzz。沿飛行測(cè)線數(shù)據(jù)的采樣間隔約為80 m，主測(cè)線間距為500 m，副測(cè)線距為5000 m，貼地起伏飛行高度為100 m。選取垂直梯度做實(shí)驗(yàn)，區(qū)域選取地質(zhì)結(jié)構(gòu)明顯的區(qū)域，范圍為緯度48 °11 ′～48 °25 ′，經(jīng)度-59 °14 ′～-58 °54 ′，區(qū)域內(nèi)共24121 個(gè)測(cè)點(diǎn)。

3.1 空間變異分析

將測(cè)線數(shù)據(jù)按照通用橫軸墨卡托投影，轉(zhuǎn)換成平面直角坐標(biāo) ?設(shè)置滯后距分段為500 m，最大搜索半徑設(shè)置為區(qū)域?qū)蔷€長(zhǎng)度的一半，約15000 m，角度容差為22.5 °，分別計(jì)算了重力梯度在方位角為0 °，45 °，90 °，135 °四個(gè)方向上的實(shí)驗(yàn)變異函數(shù)，并繪制出擬合后的球狀模型變異函數(shù)（圖4），其參數(shù)見(jiàn)表5。從圖4 和表5 中看出，45 °和90 °方向的實(shí)驗(yàn)變異函數(shù)較為理想，變異值在滯后距約10000 m 處趨于穩(wěn)定。而0 °和135 °方向的變異結(jié)構(gòu)較差，當(dāng)滯后距大于變程后，實(shí)驗(yàn)變異值并未像理論值一樣趨于穩(wěn)定，與理論變異函數(shù)偏差較大，擬合誤差分別達(dá)到5732 E2和3682 E2。

從圖4 的擬合曲線和表5 的擬合參數(shù)可知，此處的垂直梯度存在明顯的各向異性。各向異性是由于特定的地質(zhì)體構(gòu)造引起的，按其性質(zhì)可分為幾何各向異性和帶狀各向異性。幾何各向異性指的是變異函數(shù)在不同方向上具有相同的基臺(tái)值、不同的變程，即在相同距離上不同方向的兩點(diǎn)間平均變異程度不同。

從表5 中可看出，此處的垂直梯度明顯屬于幾何各向異性。二維各向異性的變程圖是以東西方向（方位角90 °）為主軸（變程9923 m），南北方向（方位角0 °）為次軸（變程4812 m）的橢圓。對(duì)于幾何各向異性的變異函數(shù)，可通過(guò)坐標(biāo)旋轉(zhuǎn)變換和構(gòu)建伸縮矩陣將各向異性結(jié)構(gòu)變換成各向同性結(jié)構(gòu)，也即將橢圓形狀的變程圖轉(zhuǎn)換成圓形的變程圖。對(duì)于本區(qū)域，旋轉(zhuǎn)角度為0 °，各向異性比例系數(shù)K=9923/4812，因此只需要將原始數(shù)據(jù)的y坐標(biāo)（指向北向）乘以比例系數(shù)K，得到新的坐標(biāo)系。在新坐標(biāo)系下，垂直梯度具有各向同性的空間結(jié)構(gòu)。

表5 不同方向擬合的球狀變異函數(shù)參數(shù)Tab.5 Spherical variogram parameters fitted in different directions

對(duì)經(jīng)過(guò)坐標(biāo)旋轉(zhuǎn)和伸縮變換后的垂直梯度，重新統(tǒng)計(jì)各向同性的實(shí)驗(yàn)變異函數(shù)，如圖5。從圖5 中看出，各向同性的空間實(shí)驗(yàn)變異結(jié)構(gòu)明顯好于各向異性的空間結(jié)構(gòu)，滯后距大于變程后，實(shí)驗(yàn)變異值逐漸趨于平穩(wěn)，與理論變異函數(shù)相符。擬合球狀模型塊金值為720 E2，基臺(tái)值438 E2，變程10190 m，擬合誤差的均方差為18.3 E2。

圖5 各向同性變異圖Fig.5 Isotropic Variogram

圖6 套合變異函數(shù)圖Fig.6 Variogram of nested fitting

3.2 降噪處理

表6 理論變異函數(shù)模型參數(shù)Tab.6 Theoretical variogram model parameters

用分解出來(lái)的變異函數(shù)，計(jì)算不同尺度的因子，插值間隔為250 m，然后將y坐標(biāo)（指向北向）除以各向異性比例系數(shù)K，恢復(fù)成舊坐標(biāo)系，見(jiàn)圖7，坐標(biāo)單位是km。

圖7(a)是原始垂向重力梯度值，由于數(shù)據(jù)中含有大量噪聲，因此掩蓋了有效信息。通過(guò)因子克里金將原始梯度信息分解成3 個(gè)尺度的因子：圖7(b)是對(duì)應(yīng)的高頻率因子，認(rèn)為它是隨機(jī)噪聲；圖7(c)是對(duì)應(yīng)的小尺度因子；圖7(d)是對(duì)應(yīng)的大尺度因子，從圖中看出，隨機(jī)噪聲能夠有效地被壓制，呈現(xiàn)出明顯的地質(zhì)特征，且地質(zhì)特征的走向信息也有很清晰的反映。

圖7 垂直重力梯度的因子克里金濾波Fig.7 Factor Kriging filtering of vertical gravity gradient

4 結(jié)論

因子克里金分析的關(guān)鍵是分解各因子的變異函數(shù)，獲取了各因子的變異函數(shù)后，就可以通過(guò)因子克里金方程組，將區(qū)域變量分解成不同尺度的分量。繼而將重力梯度分解成隨機(jī)噪聲、小尺度分量和大尺度分量，以此實(shí)現(xiàn)降噪的目的。從模擬和實(shí)測(cè)數(shù)據(jù)的實(shí)驗(yàn)可看出：

（1）因子克里金濾波可以自動(dòng)獲得各項(xiàng)濾波參數(shù)。因子克里金分析法通過(guò)對(duì)采樣數(shù)據(jù)的統(tǒng)計(jì)分析來(lái)構(gòu)建空間變異函數(shù)，從而獲取塊金值，基臺(tái)值、變程等濾波參數(shù)。該過(guò)程是完全基于數(shù)據(jù)驅(qū)動(dòng)的，不需要任何先驗(yàn)信息，而傳統(tǒng)的空間濾波器的參數(shù)需要根據(jù)經(jīng)驗(yàn)確定；

（2）因子克里金濾波直接采用離散點(diǎn)數(shù)據(jù)。傳統(tǒng)的空間濾波方法往往需要規(guī)則網(wǎng)格，而因子克里金濾波是基于離散點(diǎn)的濾波，不需要事先將離散點(diǎn)數(shù)據(jù)網(wǎng)格化；

（3）模擬實(shí)驗(yàn)表明，因子克里金降噪法相對(duì)于傳統(tǒng)的空間濾波精度有顯著的提高，重力梯度各分量精度均提升30%以上，尤其是垂直梯度分量精度提升了42%，說(shuō)明因子克里金降噪法適用于地質(zhì)特征明顯區(qū)域的梯度場(chǎng)降噪；

（4）通過(guò)實(shí)測(cè)航空重力梯度的空間變異分析發(fā)現(xiàn)，該地區(qū)的變異函數(shù)呈現(xiàn)出較強(qiáng)的各向異性，因此需要將各向異性的重力梯度轉(zhuǎn)換為各向同性的重力梯度才能采用因子克里金降噪，經(jīng)降噪后該地區(qū)的重力梯度呈現(xiàn)出明顯的地質(zhì)特征。

致謝：感謝加拿大自然資源部（Natural Resources Canada）提供的圣喬治灣實(shí)測(cè)航空重力梯度數(shù)據(jù)，允許使用此數(shù)據(jù)進(jìn)行研究并發(fā)表研究成果?