基于正則化等效源模型的航空重力梯度測量信息降噪方法

趙德軍，孫中苗，趙東明，謝心和

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

航空重力梯度測量以其快速高效、機動靈活以及超高空間分辨率的優勢，在油氣勘探、固體礦產勘查、地理信息構圖和科學研究中發揮著越來越重要的作用。根據測量重力梯度的分量情況，可以分為以FALCON 為代表的部分張量梯度測量系統和以Air-FTG 為代表的全張量梯度測量系統兩類。航空重力梯度測量能獲取多個梯度分量，然而在不破壞其內部一致性的情況下處理多分量航空重力梯度數據是一個巨大的挑戰。

等效源法利用一系列簡單的虛擬地質體代替真實場源進行重磁場建模，通過虛擬地質體的正演對重磁場信息進行重構。等效源法重構重力場的目標是使虛擬場源產生的重力場無限接近真實的位場，其目的是恢復空間場值，而不是對真實場源的準確重建，不考慮地下空間的三維密度反演[1-6]。所以，通過構建合理的、簡單的等效場源模型，可以大幅減少計算難度。一般用長方體、圓柱體或球等規則幾何體來表示等效源模型，如果規則幾何體壓縮成沒有體積的點，則等效源模型就簡化為大地測量領域的點質量模型。

等效源模型可以直接處理原始離散點上的觀測數據，能夠聯合解算不同類型、高度及分辨率的重磁數據，因此被廣泛用于重磁位場數據處理。在重磁領域，等效源法被用于位場數據的轉換[1]、延拓[4]，以及多源重磁數據的融合[6]等。等效源法非常適合于重力梯度這樣的多分量數據處理，因為每個重力梯度分量都是由相同的地質體引起的。利用等效源模型重建重力梯度，能保持多個梯度分量之間固有的內在聯系。

但是，求解虛等效源模型所涉及的矩陣求逆常受到觀測誤差的影響，使數值解算失穩，可能出現病態性問題，因此需運用正則化方法處理其法方程的病態性。Tikhonov 正-則化方法在解決病態方程上得到了廣泛的應用，文獻[7-9]的討論均以Tikhonov 正則化為基礎。截斷奇異值TSVD 方法也是解決病態性問題的一種重要方法[10]。本文將TSVD 和Tikhonov 正則化方法引入等效源模型解算，以抑制設計矩陣的病態性，實現多分量航空重力梯度的降噪，以及不同重力梯度分量之間的轉換。

1 等效源原理

1.1 基本理論

矩形棱柱體作為最常用的模型構建單元被廣泛應用于重力梯度正演數值模擬之中。在局部直角坐標系——北東下NED，即x軸指向北，y軸指向東，z垂直向下，單個矩形棱柱體場源正演重力梯度張量和重力異常（大地測量鄰域稱為擾動重力）的無奇點解析公式為[11]：

從式(1)可看出，若已知計算點的坐標和長方體的坐標，則正演的重力梯度張量與密度ρ成線性關系：代表式(1)等式右端的線性系數項。假設埋藏的等效源共n個,則第i個觀測點的梯度分量為n個等效源的和：

觀測到的梯度數據位于三維空間中不規則分布的點上，若有m個觀測數據，則有：

將式(4)寫成矩陣形式：

式中L為m階的梯度觀測值向量，X為n階的等效源密度向量，A為聯系已知觀測值和未知密度的m×n階系數矩陣。根據最小二乘理論，可在均衡多個重力梯度信息的基礎上，得到自洽的密度參數估值：

1.2 正則化改造

1.2.1 正則化方法

由于虛擬等效源解算過程可能存在欠適定性，因此若不采用特殊的方法求解，將得不到合理的估值。對觀測矩陣A進行奇異值分解，就會發現觀測方程病態的根本原因：系數矩陣A的一系列奇異值中存在減趨于零的奇異值，正是由于存在這些特別小的奇異值，導致很小的觀測噪聲也會引起待估參數較大的偏差[7]。要想獲得穩定的等效源解，必須對不適定方程進行正則化改造，以抑制觀測誤差對待估參數的影響。為了獲得穩定的參數估值，需要根據病態觀測方程的特點，構建正則化解，常用的正則化方法有截斷奇異值TSVD 和Tikhonov 正則化。

截斷奇異值TSVD 法，顧名思義，首先將病態觀測方程進行奇異值分解，然后通過截掉系數矩陣的小奇異值和對應的特征向量，最后僅利用大的奇異值和對應的特征向量構建參數估值。通過這種拋棄系數矩陣小奇異值的方法，從而避免放大高頻觀測誤差對參數估值的影響。

對系數矩陣A作奇異值分解[10]：

式中，U是A的左奇異向量矩陣，為m×n的正交矩陣；V是A的右奇異向量矩陣，為n×n的正交矩陣；Λ是n×n的對角矩陣，其對角線上的元素為A的遞減的奇異值，即。截斷奇異值就是只保留前面共k個較大的奇異值，截掉后面較小的奇異值，截掉的奇異值直接取零。則式(5)未知參數的TSVD 正則解為：

式中采用了matlab 形式的子矩陣表示方法，

Tikhonov正則化方法實質是用相鄰的適定解去逼近原問題的解，通過構造穩定泛函準則并引入正則化參數來求解穩定的參數估值。對式(5)引入Tikhonov正則化算法，其正則化準則為：

正則化解為：

式中，α為正則化參數，根據觀測數據受干擾程度確定α取值，隨噪聲由弱變強，α取值相應由小增大；W為對角矩陣，其對角元素為等效源的體積權重，若設計的等效源單元體的體積相同，則可取為單位陣。由此可見，相對于最小二乘估計，Tikhonov正則化估計是通過適當地犧牲了待估參數的有偏性來換取方差的減小。

比較TSVD和Tikhonov這兩種正則化方法可以發現，二者在本質上是一致的，基本思路都是如何消除小奇異值的影響。二者的差異在于減少小奇異值對解的影響的程度上，TSVD 方法是直接拋棄了較小的奇異值對解的影響，相當于刪除引起病態矩陣中不可靠的部分，而Tikhonov 正則化則是對矩陣進行約束，降低奇異值對結果的影響。

1.2.2 正則化參數的選取

對于TSVD 法，如何選擇截斷參數k是整個奇異值截斷法的關鍵，而Tikhonov 正則化法需要確定合理的正則化參數α。本文將k和α統稱為正則化參數，對于正則化參數的優選，有L 曲線法和廣義交叉驗證法GCV。

采用GCV 選擇最佳正則化參數，就是選擇一個參數值，使得GCV 函數最小：

2 數值模擬試驗

全張量梯度模擬。為了檢驗等效源降噪的效果，采用文獻[9]的模擬實驗，虛擬場源模型參數、模擬數據參數，均與文獻[9]相同。采用文獻[12]的無奇點公式正演計算得到空中80m飛行高度處的理論重力梯度值，將其作為真值，見圖1 第1 列。圖1 中從上到下分別是Txx、Txy、Txz、Tyy、Tyz和Tzz共6 個分量，每個梯度分量東西方向和南北方向均模擬了31個測點。然后在每個梯度分量上加上均值為0 E、標準差為5 E 的高斯白噪聲，如圖1 第2 列，可以看出噪聲已經掩蓋了信號。

部分張量梯度模擬。為了檢驗等效源轉換重力梯度分量的效果，利用圖1 第2 列中添加噪聲后的Txy、Txx和Tyy這3 個分量模擬了FALCON 部分張量梯度測量系統的觀測值——水平梯度分量和分量：

等效層設置。采用100 m 邊長的立方體作為等效源，等效層的覆蓋范圍在每個方向上均超出數據區域5 個單位的等效源，因此等效層在東西和南北方向各有41 個等效源。等效層埋藏深度設置為貼近地面，即地下0～100 m。

方案1、基于L 曲線方法確定正則化參數的TSVD 方法，確定的截斷參數k=229，如圖2(a)；

方案2、基于GCV 方法確定截斷參數的TSVD 方法，確定的截斷參數k=193，如圖2(b)；

方案3、基于L 曲線方法確定正則化參數的Tikhonov 正則化方法，確定的正則化參數α=0.015959，如圖2(c)；

分別利用這4 種正則化方案確定出等效源模型的參數后，再重構重力梯度張量的6 個分量，并與無噪聲的理論梯度值作比較，統計均方差如表1。從表1看出，4 種方案都能有效地提高濾波精度。綜合來看，方案1 重構的重力梯度張量精度略高于其他3 種方案。盡管垂直梯度分量Tzz的精度有所提高，但效果不如其它分量，這也很好理解，重力梯度在垂直方向的量級較大，水平方向的量級較小，用較微弱的水平梯度信號來推算較強的垂直梯度信號，精度提升有限。圖1 第3 列是采用FALCON 兩個水平分量按照方案1 重構的重力梯度張量，從中看出，高頻噪聲被有效地濾掉了。

為了檢驗等效源模型對全張量重力梯度的降噪效果，將圖1 第2 列的所有梯度分量作為觀測值，采用方案1 的正則化等效源降噪，見圖1 第4 列，精度統計見表1。從表1 看出，采用全張量重力梯度降噪，精度有大幅的提升，尤其是3 個直線梯度分量Txx、Tyy、Tzz的精度提升了50%以上。

圖1 重力梯度分量轉換和全張量梯度降噪圖 /EFig.1 Gravity Gradient Component Conversion and Total Tensor Gradient Denoising Diagram /E

圖2 正則化參數的確定Fig.2 Determination of Regularization Parameters

表1 不同方案降噪后的均方差（單位：E）Tab.1 Root mean square error of different noise reduction schemes（unit:E）

3 實測數據實驗

FALCON 系統只能測量重力梯度的水平分量，但對實際應用而言，重力垂直梯度和重力異常具有更大的用處。因此，需要將水平梯度轉換為垂直梯度和重力異常。為了檢驗等效源法濾波和轉換的效果，采用HeliFALCON 實測數據進行轉換試驗。2014年，法國的通用地球物理公司 CGG（Compagnie Généralede Géophysique）公司承接了美國地質勘探局USGS 的合同，采用搭載于直升飛機的HeliFALCON系統在美國密蘇里州蘇利文北部開展了航空重力梯度和地磁測量。測區中心坐標為北緯38 °09 ′，西經91 °13 ′。直升飛機貼地飛行，平均飛行高度為90 m，測線為南北向，切割線為東西向，測線間距400 m，切割線間距4000 m，數據采樣間隔約5 m，共獲得了3537.7 公里的測線數據。

圖3 重力梯度分量轉換圖 /EFig.3 Gravity Gradient Component Conversion Diagram /E

等效源法計算量大，因此提取一塊10 km×10 km的范圍做實驗。將測線數據按通用橫軸墨卡托UTM投影到平面上，然后采用最小曲率內插法，將原始離散數據網格化為100 m 間距的網格。實際上，等效源法可直接對原始點位上的觀測數據進行計算，既不要求觀測數據分布在水平面上，也不需要規則網格。這里將離散數據網格化，是為了方便繪圖展示。圖3(a)是直升飛機激光測高獲取的試驗區數字地形圖，圖3(b)是獲取的TNE重力梯度分量，圖3(c)是獲取的TUV重力梯度分量，圖中坐標單位均是km。

采用100 m 邊長的立方體作為等效源，等效層設置在平均地形高處。基于L曲線法的TSVD 正則化方法，聯合TNE和TUV兩個水平梯度分量重構垂直梯度分量Tzz（圖3(d)）和重力異常Tz（圖3(e)）。Tzz展示的是地質體密度變化的俯視圖，所以圖3(d)和圖3(a)形態上相似。Tz含有相對低頻的波長信息，所以相對于Tzz圖像更光滑。

USGS 發布了其采用頻域位場轉換方法計算的垂直分量Tzz（圖3(f)）和重力異常Tz（圖3(g)），該方法利用重力梯度分量在頻率的內在轉換關系，聯合TNE和TUV兩個分量采用傅里葉變換得到重力梯度的頻譜[13]。

對比圖3(d)和圖3(f)，以及圖3(e)和圖3(g)，兩種方法計算的Tzz和Tz形態基本一致，但是細節上還存在一些差異。USGS 采用的頻域位場轉換方法會導致噪聲放大，因此該算法附加了一個低通濾波器以壓制噪聲。

用等效源法計算的Tzz和Tz分別減去USGS 公布的Tzz和Tz，統計二者的差異，如表2所示。

表2 等效源轉換結果同USGS 轉換結果差異統計表Tab.2 Statistical table of differences between equivalent source conversion results and USGS conversion results

Tzz平均差異為0.45 E，標準差為2.68 E，Tz平均差異為-0.14 mGal，標準差為0.36 mGal，表明等效源法濾波轉換結果與UGGS 的結果吻合。

4 結論

利用矩形棱柱體正演重力梯度全張量的解析模型，構建了重力場梯度與等效源密度參數的線性方程組。構造的系數矩陣只與采樣點坐標和等效源節點坐標有關，因此觀測數據不依賴于平面網格化數據，可避免曲化平以及網格化引入的誤差。針對虛擬等效源解算過程可能存在欠適定性的問題，引入截斷奇異值TSVD 和Tikhonov 正則化方法。通過模擬和實測數據實驗得出結論：

（1）等效源降噪能保證各梯度分量之間的內在聯系。聯合多個梯度分量能反演自洽的等效源參數，能確保不同梯度分量的同源性，能確保降噪后的梯度分量滿足拉普拉斯約束；

（2）等效源法具有不同梯度分量之間轉換的功能。通過梯度分量觀測值反演等效源參數后，再利用梯度張量與等效源的解析模型，可以計算其它梯度分量，實現不同梯度分量之間的轉換。該功能對于類似FALCON 這種只能獲取水平重力梯度的系統來說，尤為重要，因為可以將水平梯度轉換成垂直梯度和重力異常；

（3）正則化算法能有效改善等效源反演的穩定性和精度。模擬實驗充分驗證了正則化算法能降低病態方程系數矩陣的特征值對觀測噪聲的敏感度，有效抑制觀測噪聲的放大效應，提高了未知參數的穩定性和精度。綜合而言，基于L曲線方法確定正則化參數的TSVD 方法精度較高。