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

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

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

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

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

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

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

1 等效源原理

1.1 基本理論

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

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

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

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

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

1.2 正則化改造

1.2.1 正則化方法

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

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

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

式中，U是A的左奇異向量矩陣，為m×n的正交矩陣；V是A的右奇異向量矩陣，為n×n的正交矩陣；Λ是n×n的對(duì)角矩陣，其對(duì)角線上的元素為A的遞減的奇異值，即。截?cái)嗥娈愔稻褪侵槐Ａ羟懊婀瞜個(gè)較大的奇異值，截掉后面較小的奇異值，截掉的奇異值直接取零。則式(5)未知參數(shù)的TSVD 正則解為：

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

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

正則化解為：

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

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

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

對(duì)于TSVD 法，如何選擇截?cái)鄥?shù)k是整個(gè)奇異值截?cái)喾ǖ年P(guān)鍵，而Tikhonov 正則化法需要確定合理的正則化參數(shù)α。本文將k和α統(tǒng)稱為正則化參數(shù)，對(duì)于正則化參數(shù)的優(yōu)選，有L 曲線法和廣義交叉驗(yàn)證法GCV。

采用GCV 選擇最佳正則化參數(shù)，就是選擇一個(gè)參數(shù)值，使得GCV 函數(shù)最?。?/p>

2 數(shù)值模擬試驗(yàn)

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

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

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

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

方案2、基于GCV 方法確定截?cái)鄥?shù)的TSVD 方法，確定的截?cái)鄥?shù)k=193，如圖2(b)；

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

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

為了檢驗(yàn)等效源模型對(duì)全張量重力梯度的降噪效果，將圖1 第2 列的所有梯度分量作為觀測(cè)值，采用方案1 的正則化等效源降噪，見(jiàn)圖1 第4 列，精度統(tǒng)計(jì)見(jiàn)表1。從表1 看出，采用全張量重力梯度降噪，精度有大幅的提升，尤其是3 個(gè)直線梯度分量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í)測(cè)數(shù)據(jù)實(shí)驗(yàn)

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

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

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

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

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

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

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

表2 等效源轉(zhuǎn)換結(jié)果同USGS 轉(zhuǎn)換結(jié)果差異統(tǒng)計(jì)表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)建了重力場(chǎng)梯度與等效源密度參數(shù)的線性方程組。構(gòu)造的系數(shù)矩陣只與采樣點(diǎn)坐標(biāo)和等效源節(jié)點(diǎn)坐標(biāo)有關(guān)，因此觀測(cè)數(shù)據(jù)不依賴于平面網(wǎng)格化數(shù)據(jù)，可避免曲化平以及網(wǎng)格化引入的誤差。針對(duì)虛擬等效源解算過(guò)程可能存在欠適定性的問(wèn)題，引入截?cái)嗥娈愔礣SVD 和Tikhonov 正則化方法。通過(guò)模擬和實(shí)測(cè)數(shù)據(jù)實(shí)驗(yàn)得出結(jié)論：

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

（2）等效源法具有不同梯度分量之間轉(zhuǎn)換的功能。通過(guò)梯度分量觀測(cè)值反演等效源參數(shù)后，再利用梯度張量與等效源的解析模型，可以計(jì)算其它梯度分量，實(shí)現(xiàn)不同梯度分量之間的轉(zhuǎn)換。該功能對(duì)于類似FALCON 這種只能獲取水平重力梯度的系統(tǒng)來(lái)說(shuō)，尤為重要，因?yàn)榭梢詫⑺教荻绒D(zhuǎn)換成垂直梯度和重力異常；

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