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

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

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

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

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

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

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

1 等效源原理

1.1 基本理論

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

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

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

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

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

1.2 正則化改造

1.2.1 正則化方法

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

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

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

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

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

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

正則化解為：

式中，α為正則化參數(shù)，根據(jù)觀測數(shù)據(jù)受干擾程度確定α取值，隨噪聲由弱變強，α取值相應(yīng)由小增大；W為對角矩陣，其對角元素為等效源的體積權(quán)重，若設(shè)計的等效源單元體的體積相同，則可取為單位陣。由此可見，相對于最小二乘估計，Tikhonov正則化估計是通過適當(dāng)?shù)貭奚舜绤?shù)的有偏性來換取方差的減小。

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

1.2.2 正則化參數(shù)的選取

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

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

2 數(shù)值模擬試驗

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

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

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

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

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

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

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

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

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

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

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

3 實測數(shù)據(jù)實驗

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

圖3 重力梯度分量轉(zhuǎn)換圖 /EFig.3 Gravity Gradient Component Conversion Diagram /E

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

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

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

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

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

表2 等效源轉(zhuǎn)換結(jié)果同USGS 轉(zhuǎn)換結(jié)果差異統(tǒng)計表Tab.2 Statistical table of differences between equivalent source conversion results and USGS conversion results

Tzz平均差異為0.45 E，標(biāo)準(zhǔn)差為2.68 E，Tz平均差異為-0.14 mGal，標(biāo)準(zhǔn)差為0.36 mGal，表明等效源法濾波轉(zhuǎn)換結(jié)果與UGGS 的結(jié)果吻合。

4 結(jié)論

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

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

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

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