不同幅度噪聲和缺失數(shù)據(jù)對大地電磁正則化反演的影響

李愛勇，曹創(chuàng)華，柳建新，郭榮文，童孝忠，柳 卓

(1. 中南大學 有色金屬成礦預測教育部重點實驗室，長沙 410083；2. 中南大學 地球科學與信息物理學院，長沙 410083；3. 有色金屬華東地質(zhì)勘查局，南京 210007)

自從物探反問題之父 BACKUS和應用數(shù)學家GILBERT提出了別具一格的 BG反演理論后，地球物理反問題成為地球物理學者研究的重要內(nèi)容。在實際應用方面，地球物理勘探中最關(guān)鍵的一個環(huán)節(jié)就是如何進行數(shù)據(jù)處理，以便能有效地識別異常體；然而數(shù)據(jù)處理時的反演技術(shù)仍是解決實際問題的關(guān)鍵。大地電磁反演問題也往往可以轉(zhuǎn)化為數(shù)學上不適定問題通過正則化方法求解，經(jīng)過幾十年來的發(fā)展，在地球物理的地震領(lǐng)域[1]、激電測井[2]、儀器儀表[3]、瞬變電磁測深[4]和垂直激電測深[6]領(lǐng)域基于正則化方法的反演已有所研究但實際應用效果還有待改進。正則化方法主要是利用改變正則化矩陣來改變反演的初始模型，然后進行最優(yōu)化反演[7]。正則化反演初始模型的選擇是否合理是反演結(jié)果的成敗基礎，在大地電磁領(lǐng)域，基于Wannamaker的正演模型成功實現(xiàn)了Occam的正則化反演[8]。

目前，國內(nèi)外對于正則化初始模型的選擇及其特點的分析較少，基本上是研究反演的過程，國內(nèi)學者柳建新等[7]、陳小斌等[9]、童孝忠[10]系統(tǒng)分析了正則化大地電磁反演的現(xiàn)狀，初步分析了基于混合遺傳算法的正則化反演，但對如何選擇最佳初始正則化模型以及在噪聲影響和缺失數(shù)據(jù)下初始模型對反演效果的影響基本上沒有討論。本文作者針對地球物理正則化反演里常用的最平坦模型進行的計算，利用改變后的正則化矩陣分別討論了最平滑模型和最小模型的反演效果，并進行了對比分析。然后就針對野外采集數(shù)據(jù)人文噪聲難以避免的情況利用高斯噪聲的形式加入到反演總體目標函數(shù)中進行了分析討論。

此外，針對最平坦、最平滑和最小化模型的反演效果，選擇這3種反演模型中較好的模型加入不同噪聲，比較對于噪聲的適應效果；然后針對某個頻段或者深度的數(shù)據(jù)缺失時進行反演；最后討論何種效果下的正則化反演是最適合在數(shù)據(jù)反演時利用。

1 大地電磁中的正則化反演技術(shù)

礦產(chǎn)地球物理勘探中，由于地形復雜，加之反演過程的復雜，結(jié)果往往不盡人意，會出現(xiàn)失真或者多解情況，這些情況在地球物理反演中被稱為不適定問題。隨著反演技術(shù)的發(fā)展，把非線性問題線性化，人們提出了廣義線性化的處理手段。但是地球物理反演中的問題大多數(shù)是非線性的，并不是每種線性化中的理論都能用于廣義的線性化反演，針對這種問題，很多學者提出了具有各自特色的方法來解決。總體上看，可以總結(jié)為求解不適定問題的一般方法，都利用一個數(shù)據(jù)集與原適定問題相近的解去逼近原始問題的解，統(tǒng)稱為正則化方法[7]。

1.1 正則化的定義

正則化的定義是[11]：假設X、Y為Hilbert空間向量，T：X→Y是緊算子，T+無界，若Ra：Y→X是一族線性有界算子，有如下關(guān)系：

如果式(1)中α為正則化因子，則稱Rα為算子T的正則化或正則化算子。

在實際問題中，方程 Tx=y的右端往往是觀測到的值y，它與y之間存在誤差，假設y ∈ D ( T+)，δδ且≤δ，Rαyδ作為方程Tx=y的精確解T+y的近似，其誤差為

差，反映了Rα與T+之間的逼近程度。當α→0時，正則誤差α→0時，

冗余誤差被放大。從Rα逼近T+的角度來說，正則因子α越小越好；但從解的穩(wěn)定性角度來說越小越好。這兩方面是矛盾的，因此，必須合適選取正則參數(shù)α以獲得較好的結(jié)果。

1.2 正則化因子的選取

目前，在地球物理反演當中，利用正則化方法來減少反演的多解性和增強過程的穩(wěn)定性是最為有效的手段。正則化因子的選取直接影響著反演的結(jié)果，其值越合適，先驗信息擬合的程度越高，就能更有效地利用野外所采集的數(shù)據(jù)。吳小平和徐果明[12]認為：在反演過程中可以選取一個折中化的因子，但這往往需要從外部輸入，不斷測算才能較好地符合我們想要的結(jié)果。CONSTABLE等[13]和 SMITH和 BOOKER[14]分別提出了Occam和RRI算法，這兩種反演方法都能進行二維自適應的反演計算，其中Occam是利用線性搜索自適應迭代法選取最佳的正則化因子[8]，實現(xiàn)了數(shù)據(jù)充分擬合和模型極度光滑；RRI是通過所建立的反演函數(shù)中的目標項和自適應因子之間的關(guān)系[8]，根據(jù)每次得到的目標項來推測自適應因子。這兩種反演算法已經(jīng)用在了實際生產(chǎn)中，效果較好，但是遇到數(shù)據(jù)量很大時，就有些力不從心，尤其是原始數(shù)據(jù)比較粗糙時，更難解決。

正則化通過人為地強加某種結(jié)構(gòu)信息(可以結(jié)合具體問題的結(jié)構(gòu)特點)，在數(shù)學上為地球物理反演問題提供唯一的、穩(wěn)定的結(jié)果。正則化因子的選擇是正則化反演的關(guān)鍵，其值過大反演結(jié)果只包含先驗正則化結(jié)構(gòu)信息；過小則只包含數(shù)據(jù)反映的結(jié)構(gòu)信息，因此，如何選取一個合適的正則化因子是非常重要的前提。目前，常用的正則化因子的求解方法有 GCV方法，L-Curve方法[10]和概率統(tǒng)計法，本文作者在加入噪聲時的討論中，采用的正則化因子的選擇方法是位對噪聲的Gaussian假設，其數(shù)據(jù)擬合服從Chi-square分布，屬于概率統(tǒng)計方法。

2 正則化矩陣的選取及其反演效果

則有：

可以得到：

反演總體目標函數(shù)可以用下式進行表示[12,15?18]：

如果先驗模型參數(shù)m?和不定的參數(shù)ξi，i=1，…，m都是可以加以利用，那么可以得到正則化矩陣如下：

設定數(shù)據(jù)和先驗信息都是對稱的矩陣，這意味著m?和標準差ξi伴隨著假設模型先驗參數(shù)的分解都是高斯分解法。如果已知的先驗信息是比較好的，一般情況下，正則化因子都被設定為α=1。這種情況下的參數(shù)設定不只是依賴于χ2(即數(shù)據(jù)擬合服從Chi-square分布)變化的。在這種情況下，即使估計的參數(shù)不確定，但通過奇異值分解的反演效果來看，與零估計相比可能仍然是比較好的選擇。

2.1 3種正則化矩陣初始模型

2.1.1 最平坦模型正則化

模型m代表著一個m(z)空間函數(shù)的一個分量，假設 0=m 和

最平坦模型是指使的函數(shù)

H ( m ? m ?)2的梯度利用拉格朗日法求取最小值時，仍然盡可能地擬合原始模型。如果Zi+1?Zi=Δ為常數(shù)，那么得到的矩陣如下：

式(8)即為最平坦模型正則化矩陣。2.1.2 最平滑模型正則化

正則化矩陣 H 能被離散成對角矩陣的形式(式(9))：

在其底部最后兩行為零。如果無法先驗知道基模型m?的結(jié)構(gòu)信息，可以假設m?=0，那么可以得到

這種加入了二階偏導數(shù)先驗信息稱之為最平滑模型正則化反演。

如果令Zi+1?Zi=Δ為常數(shù)，則有：

式(11)即為最平滑模型正則化矩陣。2.1.3 最小模型正則化

最小正則化模型的定義是假設H=I和m=0，其中I表示單位矩陣，然后再帶入反演目標函數(shù)進行正則化計算，也被稱為零階正則化阻尼最小二乘方法，其具體的實現(xiàn)過程如下：

給定生成矩陣的主對角線乘以正則化因子α，用以克服其奇異性或者改善病態(tài)條件。

令：

即得：

根據(jù) 2.1節(jié)提到的正則化初始模型，利用正則化算法和2D海洋電磁Occam反演原理[8]，通過對正則化矩陣初始的修正，針對上述加入3種不同幅度噪聲，反演結(jié)果分別如2.2.1～2.2.3小節(jié)所述，其中3種情況下的反演結(jié)果如圖1～3所示。

2.2.1 最平坦模型

如圖1所示，當深度為0～700 m時，加入3種不同噪聲后反演所得電阻率均為20 ?·m左右；當深度為700～2 000 m時，只有加入0.02%的噪聲時，反演所得電阻率為10 ?·m左右，而噪聲變大時，會掩蓋此信息；當深度為2 000～3 500 m時，隨著噪聲的增加，反演

求得其最大值λi，有

2.2 不同噪聲下的反演

對表 1中一維地層模型產(chǎn)生的合成數(shù)據(jù)(雖然誤差不切合實際，但可以更好地研究和分析地球物理反演問題)加入0.02%、2%、20%的噪聲，然后研究這3種情況下的正則化反演結(jié)果。結(jié)果越來越偏離真實結(jié)果；當深度趨于無限大時，反演結(jié)果都趨于地表表層電阻率，低于初始模型設置的電阻率。由此可見，噪聲為 0.02%時，由于噪聲干擾小，結(jié)合先驗信息可以很好地反演、重構(gòu)實際模型。隨著數(shù)據(jù)噪聲的增大(干擾信號的增加)，反演結(jié)果越來越差；當噪聲達到20%時，反演結(jié)果完全不能反映真實情況。這時反演結(jié)果大部分信息來自于最平整結(jié)構(gòu)信息(即所有相鄰參數(shù)間的變化率為最小)。因此，數(shù)據(jù)誤差較大時，即便對于簡單模型也很難重構(gòu)實際的模型結(jié)構(gòu)。從這些方面來看，在進行野外實測數(shù)據(jù)時，應盡量選用信噪比高的物探儀器，采用必要的措施減少數(shù)據(jù)誤差，這樣才有利于更好地反演野外數(shù)據(jù)。2.2.2 最平滑模型

表1 層狀合成模型參數(shù)Table 1 Parameters of layer model

圖1 加入不同噪聲最平坦模型反演結(jié)果擬合圖Fig. 1 Inversion fitting diagram of different noise with flattest model

圖2 加入不同噪聲最平滑模型反演結(jié)果擬合圖Fig. 2 Inversion fitting diagram of different noise with the smoothest model

圖3 加入不同噪聲最小模型反演結(jié)果擬合圖Fig. 3 Inversion fitting diagram of different noise with the smallest model

如圖2所示，當深度為0～700 m時，加入3種不同噪聲后反演所得電阻率均為20 ?·m左右，最平坦模型和最平滑模型具有相似的特點。當深度為 700～2 000 m時，加入0.02%的噪聲后，反演所得電阻率出現(xiàn)接近0 ?·m的情況；加入2%的噪聲后，反演所得電阻率幾乎為模型電阻率的2倍；而加入20%的噪聲后，反演所得電阻率已達模型電阻率的2倍以上；當深度為2 000～3 500 m時，加入0.02%的噪聲后，反演結(jié)果基本上符合真實值，而加入其他兩種噪聲后，反演結(jié)果都低于50 ?·m；當深度趨于無限大時，反演結(jié)果都趨于地表表層電阻率；加入20%的噪聲后，反演所得電阻率趨于0 ?·m。同樣，可以得到與最平坦模型相似的結(jié)果，隨著噪音的增加，曲線越來越不符合初始模型。當噪聲足夠小時，可以有效地反演出模型結(jié)構(gòu)信息；隨著數(shù)據(jù)噪聲的增大(干擾信號的增加)，反演結(jié)果越來越差。當噪聲達到20%時，反演結(jié)果完全不能反映真實電性結(jié)構(gòu)。這時反演結(jié)果大部分信息來自于最平滑結(jié)構(gòu)先驗信息(即所有相鄰參數(shù)間的二階偏導數(shù)最小)。與最平坦模型不一樣的是，最平滑反演后的結(jié)果更傾向于二次光滑結(jié)構(gòu)。

2.2.3 最小模型

圖3所示，當深度為0～700 m時，加入3種不同噪聲后反演所得電阻率變化較大，曲線直上直下，都在20 ?·m左右振蕩；當深度為700～2 000 m時，加入3種噪聲后，曲線都是遞減函數(shù)，3種情況相似，都從40 ?·m 降到 10 ?·m左右；當深度為 2 000～3 500 m時，加入 0.02%的噪聲后，反演結(jié)果基本上符合真實值，而加入其他兩種噪聲后，反演結(jié)果都低于50 ?·m，曲線的斜率都很?。划斏疃融呌跓o限大時，加入 0.02%的噪聲后，反演結(jié)果趨于地表表層電阻率，而加入2%和20%的噪聲后，反演結(jié)果相似，趨于30 ?·m左右?？傮w上看：隨著噪聲的增加，最小化模型的適應性越來越差，但是這種模型不像最平坦模型和最平滑模型那樣明顯受數(shù)據(jù)誤差的影響，它仍然可以反映真實模型的一些信息。因此，在以上3種反演方法中，最小化模型反演受數(shù)據(jù)誤差的影響程度最小。

2.3 “信息不完全”下的反演

當模型結(jié)構(gòu)信息“不完全”(即缺失數(shù)據(jù)時)或靈敏度矩陣病態(tài)時，為了研究以上3種正則化反演法的反演特點以及處理這些情況的方式，還是采用表1的模型實例，并在產(chǎn)生的反演數(shù)據(jù)中加入0.01%的噪聲。為了方便研究，通過以下等效方式模擬數(shù)據(jù)信息的不完全(即缺失數(shù)據(jù)時)或靈敏度矩陣的病態(tài)：在本例反演中，將靈敏度矩陣中部分列元素刪除掉用以表示電磁響應后所得數(shù)據(jù)缺失。反演結(jié)果如圖4所示。

如圖4所示，當深度為0～700 m時，3種正則化模型反演所得結(jié)果都很好，均為首層電阻率20 ?·m左右；當深度為700～2 000 m時，最平滑模型和最平坦模型反演結(jié)果為10 ?·m左右，而最小模型反演結(jié)果基本上與首層分不清，為20 ?·m左右，分辨率較低；當深度為2 000～3 500 m時，最平滑模型和最平坦模型反演結(jié)果都是50 ?·m左右，曲線圓滑，而最小模型在缺失數(shù)據(jù)時，相對高阻變?yōu)? ?·m左右；當深度趨于無限大時，3種正則化反演初始模型反演結(jié)果形態(tài)相似，都為40 ?·m左右，較為真實。由以上分析可知，對于最平坦模型反演，可采用線性插值相鄰反演結(jié)構(gòu)來獲得缺失信息部分的信息；對于最平滑模型反演，只是按照相鄰節(jié)點信息圓滑了缺失的數(shù)據(jù)；對于最小化模型反演，只采用了基模型m=0的信息。從以上反演結(jié)果可知，當某部分結(jié)構(gòu)信息在數(shù)據(jù)中反映不明顯時，這3種模型在反演中對這些情況的處理方式不一樣。但是具體哪種正則化反演效果好，需要具體問題具體分析，比如對于連續(xù)的大地地電參數(shù)反演，最平坦模型和最平滑模型的反演效果可能會好些，而最小化模型的反演效果可能較差。

3 結(jié)論

1) 選擇合適的正則化反演初始模型，對反演結(jié)果的影響至關(guān)重要，最平坦模型和最平滑模型反演效果較好，最小化模型反演效果一般。

2) 噪聲信息在反演中的影響起著決定作用，對于最平坦模型，噪聲較小時反演效果最好，隨著噪聲的增加，反演效果逐漸變差；對于最平滑模型，噪聲較小時反演效果很好，基本跟模型十分吻合，但隨著噪聲的增加，跳變比最平坦模型變化復雜，不能正確的反映初始模型；對于最小化模型，不管噪聲多大，反演效果都不理想。

圖4 信息“不完全”加入0.01%噪聲時3種正則化模型反演結(jié)果圖Fig. 4 Fitting diagram of regularization inversion with three missing data model by 0.01% noise

3) 缺失數(shù)據(jù)信息時，各種初始化模型在噪聲較小的情況下，均能反映真實的地電斷面情況，只是最小化模型反演出的數(shù)據(jù)跳變很大，有時難以辨認其真實特征。

致謝：

本文成文過程中得到了“有色資源與地質(zhì)災害探查”湖南省重點實驗室全體成員的協(xié)助，在此表示最衷心的謝意。

[1] 王振宇, 劉國華. 走時層析成像的迭代 Tikhonov正則化反演研究[J]. 浙江大學學報: 工學版, 2009, 37(12): 1685?1690.

WANG Zhen-yu, LIU Guo-hua. Study of travel time tomography by iterative Tikhonov regularization inversion [J].Journal of Zhe Jiang University: Engineering Science, 2009,37(12): 1685?1690.

[2] 汪宏年, 陶宏根, 其木蘇榮, 楊善德. 水平層狀介質(zhì)中雙側(cè)向資料的全參數(shù)正則化迭代反演與應用[J]. 地球物理學報,2002, 45(增刊): 387?399.

WANG Hong-nian, TAO Hong-gen, QI Mu-su-rong, YANG Shan-de. Regularized entire-parameter iterative inversion of dual-later log in horizontally stratify flied media and its application [J]. Chinese Journal of Geophysics, 2002, 45(Suppl):387?399.

[3] 郝燕玲, 成 怡, 孫 楓, 劉繁明. Tikhonov正則化向下延拓算法仿真實驗研究[J]. 儀器儀表學報, 2008(3): 225?228.

HAO Yan-ling, CHENG Yi, SUN Feng, LIU Fan-ming.Tikhonov regularization downward continuation algorithm simulation experimental study [J]. Chinese Journal of Scientific Instrument, 2008(3): 225?228.

[4] 毛立峰, 王緒本, 陳 斌. 直升機航空瞬變電磁自適應正則化一維反演方法研究[J]. 地球物理學進展, 2011, 26(1):300?305.MAO Li-feng, WANG Xu-ben, CHEN Bin. Study on an adaptive regularized 1D inversion method of helicopter TEM data [J]. Progress in Geophysics, 2011, 26(1): 300?305.

[5] 劉海飛, 柳建新, 阮百堯, 張賽民. 垂直激電測深二維自適應正則化反演[J]. 同濟大學學報: 自然科學版, 2009, 37(12):1685?1690.

LIU Hai-fei, LIU Jian-xin, RUAN Bai-yao, ZHANG Sai-min.2D self adaptive regularization inversion with vertical induced polarization sounding data [J]. Journal of Tongji University:Natural Science, 2009, 37(12): 1685?1690.

[6] WANNAMAKER P E, STODT J A, RIJO L. Two-dimensional topographic responses in magnetotelluric model using finite elements [J]. Geophysics, 1986, S1: 2131?2144.

[7] 柳建新, 孫 婭, 郭榮文, 童孝忠, 王 浩, 郭振威. 基于TSVD求解大地電磁不適定反演問題[J]. 中國科技論文在線,2009, 21(2): 2284?2290.

LIU Jian-xin, SUN Ya, GUO Rong-wen, TONG Xiao-zhong,WANG Hao, GUO Zhen-wei. Based on TSVD algorithm to solve the ill-posed inversion of magnetotelluric problems [J].China Science Paper Online, 2009, 21(2): 2284?2290.

[8] WANNAMAKER P E. Occam2DMT [EB/OL]. [2011?03?15].http://marineemlab.ucsd.edu/occam.

[9] 陳小斌, 趙國澤, 湯 吉, 詹 艷, 王繼軍. 大地電磁自適應正則化反演算法[J]. 地球物理學報, 2005, 48(4): 937?946.

CHEN Xiao-bin, ZHAO Guo-ze, TANG Ji, ZHAN Yan, WANG Ji-jun. Magnetotellurics adaptive regularized inversion algorithm[J]. Chinese Journal of Geophysics, 2005, 48 (4): 937?946.

[10] 童孝忠. 大地電磁測深有限單元法正演與混合遺傳算法正則化反演研究[D]. 長沙: 中南大學, 2008: 70?104.

TONG Xiao-zhong. Research of forward using finite element method and regularized inversion using hybrid genetic algorithm in magnetotelluric sounding [D]. Changsha: Central South University, 2008: 70?104.

[11] 朱華平. 不適定問題的正則化理論及其應用[D]. 武漢: 武漢理工大學, 2007: 10?32.

ZHU Hua-ping. Regularization theory and its application [D].Wuhan: Wuhan University of Technology, 2007: 10?32.

[12] 吳小平, 徐果明. 大地電磁數(shù)據(jù)的 Occam 反演改進[J]. 地球物理學報, 1998, 41(4): 547?554.

WU Xiao-ping, XU Guo-ming. The Occam inversion of magnetotelluric data [J]. Chinese Journal of Geophysics, 1998,41(4): 547?554.

[13] CONSTABLE SC, PARKERRL, CONSTABLE C G. Occam’s inversion: A practical algorithm for generating smooth models from electromagnetic sounding data [J]. Geophysics, 1987, 52(3):289?300.

[14] SMITH J T, BOOKER J R. Rapid inversion of two- and three-dimensional magnetotelluric data [J]. Geophysics, 1991, 96:3905?3992.

[15] 傅淑芳, 朱仁益. 地球物理反演問題[M]. 北京: 地震出版社,1998: 1?67.

FU Shu-fang, ZHU Ren-yi. The geophysical inversion problem[M]. Beijing: Earthquake Press, 1998: 1?67.

[16] 楊文采. 地球物理反演的理論與方法[M]. 北京: 地質(zhì)出版社,1996: 1?112.

YANG Wen-cai. Theory and method of geophysical inversion[M]. Beijing: Geological Publishing House, 1996: 1?112.

[17] 湯井田, 任政勇, 化希瑞. 地球物理學中的電磁場正演與反演[J]. 地球物理學進展, 2007, 22(4): 1181?1194.

TANG Jing-tian, REN Zheng-yong, HUA Xi-rui. Electromagnetic forward modeling and inversion in geophysical [J].Progress in Geophysics, 2007, 22(4): 1181?1194.

[18] 吉洪諾夫·A·N, 阿爾先寧·V·Y. 不適定問題的解法[M].王秉枕譯. 北京: 地質(zhì)出版社, 1979: 1?221.

TIKHONOV A N, ARSENIN V Y. Ill-posed problem solution[M]. WANG Bing-zhen, transl. Beijing: Geological Publishing House, 1979: 1?221.