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

李愛勇，曹創華，柳建新，郭榮文，童孝忠，柳 卓

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

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

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

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

1 大地電磁中的正則化反演技術

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

1.1 正則化的定義

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

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

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

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

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

1.2 正則化因子的選取

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

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

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

則有：

可以得到：

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

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

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

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

2.1.1 最平坦模型正則化

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

最平坦模型是指使的函數

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

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

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

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

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

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

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

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

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

令：

即得：

根據 2.1節提到的正則化初始模型，利用正則化算法和2D海洋電磁Occam反演原理[8]，通過對正則化矩陣初始的修正，針對上述加入3種不同幅度噪聲，反演結果分別如2.2.1～2.2.3小節所述，其中3種情況下的反演結果如圖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中一維地層模型產生的合成數據(雖然誤差不切合實際，但可以更好地研究和分析地球物理反演問題)加入0.02%、2%、20%的噪聲，然后研究這3種情況下的正則化反演結果。結果越來越偏離真實結果；當深度趨于無限大時，反演結果都趨于地表表層電阻率，低于初始模型設置的電阻率。由此可見，噪聲為 0.02%時，由于噪聲干擾小，結合先驗信息可以很好地反演、重構實際模型。隨著數據噪聲的增大(干擾信號的增加)，反演結果越來越差；當噪聲達到20%時，反演結果完全不能反映真實情況。這時反演結果大部分信息來自于最平整結構信息(即所有相鄰參數間的變化率為最小)。因此，數據誤差較大時，即便對于簡單模型也很難重構實際的模型結構。從這些方面來看，在進行野外實測數據時，應盡量選用信噪比高的物探儀器，采用必要的措施減少數據誤差，這樣才有利于更好地反演野外數據。2.2.2 最平滑模型

表1 層狀合成模型參數Table 1 Parameters of layer model

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

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

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

2.2.3 最小模型

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

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

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

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

3 結論

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

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

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

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

致謝：

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

[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] 汪宏年, 陶宏根, 其木蘇榮, 楊善德. 水平層狀介質中雙側向資料的全參數正則化迭代反演與應用[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] 吳小平, 徐果明. 大地電磁數據的 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]. 北京: 地質出版社,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].王秉枕譯. 北京: 地質出版社, 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.