基于因子克里金分析方法的航空重力梯度降噪

趙德軍，孫中苗，趙東明

（1.信息工程大學地理空間信息學院，鄭州 450001；2.西安測繪總站，西安 710054；3.地理信息工程國家重點實驗室，西安 710054；4.西安測繪研究所，西安 710054）

相對于航空標量重力測量和矢量重力測量而言，航空重力梯度測量有兩大顯著優勢[1-3]：1）最大的優勢在于其不受載體運動和厄特弗斯效應的影響；2）重力梯度以重力場的曲率描述重力場區域結構，包含地球物理學和大地測量中急需的局部重力場幾何信息（水準面的平均曲率、鉛垂線的曲率等），因而重力梯度測量更能反映重力場的精細結構，敏感重力場的短波變化。

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

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

因子克里金法具有多尺度分析等諸多功能和優點，因而廣泛應用于地球物理與地球化學勘探[4]，圖像處理與分析[5,6]，水土環境監測[7,8]等領域。但在重力梯度降噪方面未見相關報道，因此本文首次將其應用于航空重力梯度數據降噪處理，取得了顯著的效果。模擬試驗表明，相對于傳統的低通濾波器，重力梯度各分量精度平均提升了36%，垂直梯度分量尤為顯著，精度提升了42%。

1 因子克里金降噪法

1.1 變異函數

克里金理論，就是用變異函數來描述區域化變量以解決有關地學問題的理論，其核心是變異函數。區域化變量的結構分析，就是通過區域化變量有限的空間觀測值來構建相應的理論變異函數模型，以表征該變量的主要結構特征。

設()Z x是區域化變量，有了二階平穩假設或本征假設，就可以給出變異函數的計算公式[5]：

1.2 因子克里金

因子克里金將區域化變量分解成多個尺度的因子，而每個因子都對應一個變異函數。因子克里金濾波的思路：從實驗變異函數中，將各因子對應的變異函數提取出來，認為小尺度變異函數對應的因子是需要濾掉的噪聲，保留大尺度的變異函數對應的因子作為有效信號，從而實現濾波的過程。區域化變量可以認為是不同地質體相互作用的疊加，因此將區域化變量分解為[5]：

根據估計方差最小原則，對每個因子()k Y x導出克里金方程組[5-7]：

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

理論變異函數通常是非線性模型，因此在采用加權最小二乘擬合之前，需要將非線性模型轉換為線性模型。一般采用泰勒級數將非線性模型轉換為線性模型，也可以采用多項式逼近的方法將其轉換成線性模型。另外，除了從數學上將非線性模型線性化外，還可以利用一些數學軟件里的遺傳算法、模擬退火法等最優方法來解決非線性回歸問題[9]。

2 模擬數據分析

對于模擬實驗，可以采用噪聲衰減因子β來估計噪聲衰減量。用X表示模型的真實梯度值，Xn表示增加噪聲后的梯度值，Xr為經過濾波處理后的梯度值:

式中var 表示計算的方差。當β為正時，說明噪聲被濾除，β的大小表示去除噪聲的能力，當β趨于1，表示Xr中幾乎不含噪聲成分。

利用下面的均方根誤差來表示濾波后數據的精度：

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

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

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

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

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

表1 模擬樣本點的坐標和梯度值Tab.1 Coordinates and gradient values of simulated points

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

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

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

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

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

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

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

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

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

3 實測數據分析與降噪處理

數據來源于貝爾地理空間公司在加拿大圣喬治灣測量的重力垂直梯度Tzz。沿飛行測線數據的采樣間隔約為80 m，主測線間距為500 m，副測線距為5000 m，貼地起伏飛行高度為100 m。選取垂直梯度做實驗，區域選取地質結構明顯的區域，范圍為緯度48 °11 ′～48 °25 ′，經度-59 °14 ′～-58 °54 ′，區域內共24121 個測點。

3.1 空間變異分析

將測線數據按照通用橫軸墨卡托投影，轉換成平面直角坐標 ?設置滯后距分段為500 m，最大搜索半徑設置為區域對角線長度的一半，約15000 m，角度容差為22.5 °，分別計算了重力梯度在方位角為0 °，45 °，90 °，135 °四個方向上的實驗變異函數，并繪制出擬合后的球狀模型變異函數（圖4），其參數見表5。從圖4 和表5 中看出，45 °和90 °方向的實驗變異函數較為理想，變異值在滯后距約10000 m 處趨于穩定。而0 °和135 °方向的變異結構較差，當滯后距大于變程后，實驗變異值并未像理論值一樣趨于穩定，與理論變異函數偏差較大，擬合誤差分別達到5732 E2和3682 E2。

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

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

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

對經過坐標旋轉和伸縮變換后的垂直梯度，重新統計各向同性的實驗變異函數，如圖5。從圖5 中看出，各向同性的空間實驗變異結構明顯好于各向異性的空間結構，滯后距大于變程后，實驗變異值逐漸趨于平穩，與理論變異函數相符。擬合球狀模型塊金值為720 E2，基臺值438 E2，變程10190 m，擬合誤差的均方差為18.3 E2。

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

圖6 套合變異函數圖Fig.6 Variogram of nested fitting

3.2 降噪處理

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

用分解出來的變異函數，計算不同尺度的因子，插值間隔為250 m，然后將y坐標（指向北向）除以各向異性比例系數K，恢復成舊坐標系，見圖7，坐標單位是km。

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

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

4 結論

因子克里金分析的關鍵是分解各因子的變異函數，獲取了各因子的變異函數后，就可以通過因子克里金方程組，將區域變量分解成不同尺度的分量。繼而將重力梯度分解成隨機噪聲、小尺度分量和大尺度分量，以此實現降噪的目的。從模擬和實測數據的實驗可看出：

（1）因子克里金濾波可以自動獲得各項濾波參數。因子克里金分析法通過對采樣數據的統計分析來構建空間變異函數，從而獲取塊金值，基臺值、變程等濾波參數。該過程是完全基于數據驅動的，不需要任何先驗信息，而傳統的空間濾波器的參數需要根據經驗確定；

（2）因子克里金濾波直接采用離散點數據。傳統的空間濾波方法往往需要規則網格，而因子克里金濾波是基于離散點的濾波，不需要事先將離散點數據網格化；

（3）模擬實驗表明，因子克里金降噪法相對于傳統的空間濾波精度有顯著的提高，重力梯度各分量精度均提升30%以上，尤其是垂直梯度分量精度提升了42%，說明因子克里金降噪法適用于地質特征明顯區域的梯度場降噪；

（4）通過實測航空重力梯度的空間變異分析發現，該地區的變異函數呈現出較強的各向異性，因此需要將各向異性的重力梯度轉換為各向同性的重力梯度才能采用因子克里金降噪，經降噪后該地區的重力梯度呈現出明顯的地質特征。

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