戴亦軍，童孝忠，張連偉，曹創(chuàng)華

（1.中南大學(xué) 信息物理工程學(xué)院，長沙 410083；2.鐵道第三勘察設(shè)計(jì)院，天津 300251）

利用一維正則化反演進(jìn)行大地電磁測深數(shù)據(jù)擬二維反演解釋

戴亦軍1，童孝忠1，張連偉2，曹創(chuàng)華1

（1.中南大學(xué) 信息物理工程學(xué)院，長沙 410083；2.鐵道第三勘察設(shè)計(jì)院，天津 300251）

大地電磁測深的反演問題是不適定的，其反演結(jié)果不穩(wěn)定，且具有非唯一性。通過在目標(biāo)函數(shù)中采用正則化方法，可以使得不適定反演問題具有穩(wěn)定的反演結(jié)果，并改善解的穩(wěn)定性和非唯一性問題。為了提高野外大地電磁測深數(shù)據(jù)的處理效率和初步解釋的精度，提出了大地電磁測深數(shù)據(jù)的一維正則化反演進(jìn)行擬二維反演解釋方法。這里所述的大地電磁測深一維反演解釋，與以往的解釋方法不同，其思路首先用Bostick反演的深度來控制層參數(shù)，使反演計(jì)算的模型參數(shù)僅存在電阻率；最后采用阻尼高斯－牛頓算法進(jìn)行反演計(jì)算，并將Bostick反演結(jié)果作為反演計(jì)算的初始模型。通過模型試算，結(jié)果表明其處理速度快、解釋直觀，對野外大地電磁測深數(shù)據(jù)進(jìn)行初步反演解釋是可行的。

大地電磁測深；反演；正則化方法；阻尼高斯～牛頓算法

0 前言

大地電磁測深法（MT）是前蘇聯(lián)學(xué)者Tikhonov［1］和 Cagniard［2］于上世紀(jì)五十年代提出來的，是利用天然交變電磁場研究地球電性結(jié)構(gòu)的一種地球物理勘探方法。由于它不用人工供電，成本低，工作方便，不受高阻層屏蔽，對低阻層分辨率高，且勘探深度僅與電磁場的頻率有關(guān)，淺可達(dá)幾十米，深可達(dá)數(shù)百公里，因此在許多領(lǐng)域內(nèi)得到了廣泛的應(yīng)用。MT資料反演解釋就是根據(jù)地表測得的地球電磁場響應(yīng)，如視電阻率、阻抗相位、表面阻抗、傾子等，通過一定的優(yōu)化處理求得一個合理的地電模型［3］。目前二維數(shù)據(jù)處理技術(shù)日趨成熟，比較典型的二維反演解釋方法有：OCCAM法［4］、RRI法［5］、REBOCC 法［6］、NLCG 法［7］和 ABIC法［8］。同時，一維反演方法也很多，應(yīng)用廣泛的有：OCCAM反演［9］、一維連續(xù)介質(zhì)反演的曲線對比法［10］、“正演修正法”一維反演［11］、二次函數(shù)逼近非線性反演［12］和自適應(yīng)正則化反演［13］。

另外，一維反演通常是二維反演計(jì)算的第一步驟，并且可以采用一維反演方法對地質(zhì)模型進(jìn)行擬二維反演解釋，往往能對MT數(shù)據(jù)解釋帶來有用的效果。擬二維反演方法的精度比不上二維反演［14］，但由于其反演速度快、內(nèi)存需求少，所以在實(shí)際工作中，可以對野外大地電磁測深數(shù)據(jù)進(jìn)行初步解釋。因此，使用一維反演方法進(jìn)行擬二維反演解釋還是非常有意義的。

作者在本文采用Bostick反演的深度來控制層參數(shù)，使反演計(jì)算的模型參數(shù)僅存在電阻率，并將Bostick反演結(jié)果作為反演計(jì)算的初始模型。

1 MT一維正演

假定地電剖面是水平分層均勻的，如果地電剖面共有N層，則共有2 N－1個參數(shù)，hi（i＝1，2，…，N－1）和ρi（i＝1，2，…，N）分別代表第i層的厚度和電阻率。

對于上述的一維層狀介質(zhì)模型，計(jì)算視電阻率ρa(bǔ)和相位φa的公式如式（1）。

其中 Z1表示第一層地表的波阻抗；μ為導(dǎo)磁率；為角頻率；Z可用下面的遞推公式計(jì)算。1

0ii層的特征阻抗；Zi是第i層頂面的波阻抗。

由此可見，對N層地電斷面而言，若固定每一層的厚度，則視電阻率ρa(bǔ)和相位φa可表示為關(guān)于信號周期及層電阻率參數(shù)之間的函數(shù)，即

2 MT擬二維反演

2.1 反演問題描述

大地電磁反演問題，可以抽象地描述為已知觀測數(shù)據(jù)求之于對應(yīng)模型的過程。設(shè)d是觀測數(shù)據(jù)向量，m是模型參數(shù)向量，F(xiàn)為把地球模型映射到理想數(shù)據(jù)的函數(shù)，則

式中 F為正演函數(shù)。

大地電磁反演問題是不適定的，其反演結(jié)果是不穩(wěn)定的，并且具有非唯一性，這也就意味著不同的地電模型對觀測數(shù)據(jù)的擬合具有同樣的精度。在地球物理反演中，為了改善解的穩(wěn)定性和非唯一性問題，通常是引入Tikhonov的正則化思想［15、16］，見式（6）。

其中 Pα（m）為總目標(biāo)函數(shù)；α為正則化因子；φ（m）為觀測數(shù)據(jù)與預(yù)測數(shù)據(jù)之差的平方和（即數(shù)據(jù)目標(biāo)函數(shù)）；s（m）為穩(wěn)定器（即模型約束目標(biāo)函數(shù)），這里采用基于先驗(yàn)?zāi)Ｐ偷淖钚∧Ｐ图s束。

因此，大地電磁反演問題的總目標(biāo)函數(shù)表示為式（7）。

式中 Wd為數(shù)據(jù)權(quán)系數(shù)矩陣；Wm為模型權(quán)系數(shù)矩陣；mref為先驗(yàn)?zāi)Ｐ汀?/p>

將F（m）用泰勒公式展開為式（8）。

其中 mk為模型的第k次迭代值。

于是可得：這里 dk＝F（mk）；△m＝mk＋1－mk；Jk是雅克比靈敏度矩陣。

于是有：

2.2 擬二維反演實(shí)現(xiàn)

對實(shí)測數(shù)據(jù)進(jìn)行Bostick反演，求取地電模型的層厚度，并且使其在整個反演過程中保持不變，使反演計(jì)算的模型參數(shù)僅存在電阻率。這樣層厚度參數(shù)不參與反演計(jì)算，即可實(shí)現(xiàn)一維正則化自動迭代反演計(jì)算。

線性方程組（11）求解采用阻尼～高斯牛頓算法，并將Bostick反演的電阻率作為反演計(jì)算的初始模型參數(shù)，將同一斷面的多條MT測深曲線進(jìn)行一維自動迭代反演，并把反演結(jié)果繪制成擬斷面圖，對其進(jìn)行二維地質(zhì)推斷解釋。

2.2.1 靈敏度矩陣的計(jì)算方法

在解決參數(shù)反演問題中的難點(diǎn)之一，就是計(jì)算觀測數(shù)據(jù)關(guān)于模型參數(shù)的偏導(dǎo)數(shù)，即靈敏度矩陣，它們反映了參數(shù)的相對變化對觀測數(shù)據(jù)相對變化的影響情況。同時它們還在模型參數(shù)與正演模擬響應(yīng)之間，建立起一種線性關(guān)系。這里，我們采用擾動法來計(jì)算靈敏度矩陣，即利用一階有限差分公式［17］（12）來近似計(jì)算它們關(guān)于第k個參數(shù) △mk的擾動，并通過正演問題來活動擾動后的大地電磁響應(yīng)。

2.2.2 正則化因子的選取

在總目標(biāo)函數(shù)中，正則化因子α為數(shù)據(jù)目標(biāo)函數(shù)與模型約束目標(biāo)函數(shù)的權(quán)重系數(shù)，它的大小決定了反演的擬合對象；α過大主要擬合先驗(yàn)?zāi)Ｐ?，α過小則主要擬合觀測數(shù)據(jù)。因此，正則化反演問題的關(guān)鍵，在于正則化因子的選取方法。在這里，我們采用廣義交叉驗(yàn)證方法（Generalized Cross Validation，GCV）［18］來選取正則化因子，其正則化因子α的廣義交叉驗(yàn)證函數(shù)可以寫為式（13）。

其中

利用廣義交叉驗(yàn)證方法，確定“最優(yōu)”的正則化因子α，也就是尋找使GCV（α）函數(shù)達(dá)到極小時的α值。但該方法有時并不像我們想象的那樣理想，在許多實(shí)際問題中，GCV函數(shù)在達(dá)到極小點(diǎn)時過于平坦，難于確定哪一點(diǎn)為最小值點(diǎn)。于是，我們采用類似離差原理的選取方法，將第k次迭代的正則化因子表示為式（16）。

αk＝ max（cαk－1，α＊） （16）其中 0.01≤c≤0.5；α＊為GCV函數(shù)的最小值。

3 算例

3.1 一維模型反演

選用K型地電模型正演理論數(shù)據(jù)進(jìn)行反演計(jì)算，其真實(shí)的模型參數(shù)為：

在反演計(jì)算中，將Bostick反演結(jié)果作為初始模型，先驗(yàn)?zāi)Ｐ腿ref＝0Ω·m，數(shù)據(jù)權(quán)系數(shù)矩陣取為單位矩陣，迭代次數(shù)為十次，圖1和圖2（見下頁）顯示了反演結(jié)果。

從圖1可以看出，反演結(jié)果與實(shí)際的地電模型吻合較好。相對于Bostick反演，一維反演的結(jié)果更加準(zhǔn)確地反映出中間高阻層的位置及大小。

從數(shù)據(jù)的擬合角度（見下頁圖2）看，一維反演的結(jié)果對理論數(shù)據(jù)的擬合非常好，也優(yōu)于Bostick反演的擬合曲線。

3.2 二維模型反演

構(gòu)造的地電模型如下頁圖3所示。在電阻率為50Ω·m的圍巖中，存在電阻率為5Ω·m的低阻異常體，異常體頂部埋深1 000m。采用雙二次插值的有限單元法進(jìn)行正演計(jì)算［19］，模擬測點(diǎn)數(shù)為25個，采用40個記錄頻點(diǎn)（0.01Hz～100Hz），并繪制出視電阻率和阻抗相位擬斷面圖（如下頁圖4所示）。

圖1 三層地電模型的真實(shí)結(jié)果與反演結(jié)果對比圖Fig.1 Comparison of the true and inversion results of three－layer geo－electrical model

圖2 大地電磁響應(yīng)的真實(shí)結(jié)果與反演結(jié)果對比圖Fig.2 Comparison of the true and inversion results of MT response

作者在對地電模型正演產(chǎn)生的TE模式數(shù)據(jù)和TM模式數(shù)據(jù)進(jìn)行了一維反演，并繪制出視電阻率擬二維反演等值線圖（見下頁圖5）。可以看出，雖然反演的電阻率參數(shù)值還存在較大誤差，但異常體的形態(tài)還是得到了較好的反映。

圖3 地電模型Fig.3 The geo－electrical model

4 結(jié)論

（1）大地電磁測深反演問題是不適定的，其反演結(jié)果是不穩(wěn)定的，并且具有非唯一性。這也就意味著不同的地電模型，對觀測數(shù)據(jù)的擬合具有同樣的精度。在目標(biāo)函數(shù)中采用正則化方法，可以使得不適定反演問題，具有穩(wěn)定的反演結(jié)果，改善了解的穩(wěn)定性和非唯一性問題。

（2）與傳統(tǒng)的一維大地電磁測深反演方法相比，由于將Bostick反演結(jié)果作為反演計(jì)算的初始模型，因此不需要人工輸入初始模型參數(shù)，從而提高了工作效率。

（3）雖然擬二維反演方法的精度比不上二維反演，但由于其反演速度快、內(nèi)存需求少，所以在實(shí)際工作中，可以對野外大地電磁測深數(shù)據(jù)進(jìn)行初步解釋。

［1］ TIKNONOV A N.On determining electrical charac－teristics of the deep layers of the earth’s crust［J］.Deki Akud Nuck，1950（73）：295.

［2］ CAGNIARD L.Basic theory of the magnetotelluric methods of geophysicalprospecting［J］.Geophysics，1953（18）：605.

［3］ 陳樂壽，劉任，王天生.大地電磁測深資料處理與解釋［M］.北京：石油工業(yè)出版社，1989.

［4］ DEGROOT－HEDLIN C，CONSTABLE S.Occam’s inversion to generate smooth，two－dimensional models from magnetotelluric data［J］.Geophysics，1990，55（12）：1613.

［5］ 高才坤，湯井田，王燁，等.基于RRI反演的高頻大地電磁測深在深邊部礦產(chǎn)勘探中的試驗(yàn)研究［J］.地球物理學(xué)進(jìn)展，2009，24（1）：309.

［6］ SIRIPUNVARAPORN W，EGBERT G.An efficient data－subspace inversion method for 2－D magnetotelluric data［J］.Geophysics，2000，65（3）：791.

［7］ RODI W，MACKIE R L.Nonlinear conjugate gradient algorithm for 2－D magnetotelluric inversion［J］.Geophysics，2001，66（1）：174.

［8］ UCHIDA T.Smooth 2－D inversion for magnetotelluric for magnetotelluric data based on statistical criterion ABIC［J］.Journal of Geomagnetism and Geoelectricity，1993，45（9）：841.

［9］ 吳小平，徐果明.大地電磁數(shù)據(jù)的Occam反演改進(jìn)［J］.地球物理學(xué)報(bào)，1998，41（4）：547.

［10］徐世浙，劉斌.大地電磁一維連續(xù)介質(zhì)反演的曲線對比法［J］.地球物理學(xué)報(bào)，1995，38（5）：676.

［11］蘇朱劉，羅延鐘，胡文寶.大地電磁測深“正演修正法”一維反演［J］.石油地球物理勘探，2002，37（2）：138.

［12］嚴(yán)良俊，胡文寶.大地電磁測深資料的二次函數(shù)逼近非線性反演［J］.地球物理學(xué)報(bào)，2004，47（5）：935.

［13］陳小斌，趙國澤，湯吉，等.大地電磁自適應(yīng)正則化反演算法［J］.地球物理學(xué)報(bào)，2005，48（4）：937.

［14］TONG X Z，LIU J X，XU L H.Damped Gauss－Newton optimization algorithm for two－dimensional magnetotelluric regularization inversion［A］.International Conference on Information Engineering and Computer Science，2009.

［15］TIKHONOV A N，ARSENIN V Y.Solution of illposed problems［M］.New York：Wiley，1977.

［16］FARQUHARSON C G，OLDENBURG D W.A comparison of automatic technique for estimating the regularization parameter in non－linear inverse problems［J］.Geophysical Journal International，2004，156（3）：411.

［17］FARQUHARSON，C G.Approximate sensitivities for the multi－dimensional electromagnetic inverse problem［D］.University of Edinburgh，1990.

［18］HABER E，OLDENBURG D.A GCV based method for nonlinear ill－posed problems［J］.Computational Geosciences，2000，4（1）：41.

［19］柳建新，蔣鵬飛，童孝忠，等.不完全LU分解預(yù)處理的BICGSTAB算法在大地電磁二維正演模擬中的應(yīng)用［J］.中南大學(xué)學(xué)報(bào)：自然科學(xué)版，2009，40（2）：484.

P 631.3＋25

A

1001—1749（2012）01—0033—06

湖南省自然科學(xué)基金項(xiàng)目（10JJ6059）；湖南省科研條件創(chuàng)新專項(xiàng)基金項(xiàng)目（2010TT2056）

2011－08－31

戴亦軍（1975－），男，博士，主要從事電磁法正演模擬與反演成像。