葉益信， 鄧居智， 楊海燕， 李澤林， 何梅興

(1.東華理工大學(xué)放射性地質(zhì)與勘探技術(shù)國防重點(diǎn)學(xué)科實(shí)驗(yàn)室，江西撫州 344000;2.中國地質(zhì)科學(xué)院地球物理地球化學(xué)勘查研究所，河北廊坊 065000)

目前，MT數(shù)據(jù)的正反演方法取得了很大進(jìn)展(吳信民等，2013)，二維反演方法更是百花齊放，種類繁多，但是大多數(shù)二維反演都得到光滑模型結(jié)果，如 OCCAM、RRI、REBOCC、NLCG 等(Constable et al.，1987;deGroot-Hedlin et al.，1990;Smith et al.，1991;Siripunvaraporn et al.，2000;Rodi et al.，2001)，這些方法對邊界位置的擬合較差，尤其是石油勘探中常常遇到大塊均勻電導(dǎo)率的電性單元電阻率差異較大情況，如何快速穩(wěn)定地得到反演結(jié)果和更清晰的電性分界面仍是當(dāng)前MT反演研究的一個(gè)重點(diǎn)問題。

為了克服光滑模型反演結(jié)果邊界位置模糊不清的問題，Smith等(1999)提出了反演二維MT數(shù)據(jù)的SBI方法，deGroot-Hedlin(2004)研究了用于反演二維MT數(shù)據(jù)的SBI算法，并認(rèn)為這種算法的優(yōu)點(diǎn)是能得到不同層的下底邊界深度及其電阻率值，但這種算法的收斂性較差，除非給一個(gè)好的初始模型;Zhang等(2010)對該方法做了改進(jìn)研究，通過引入對角梯度支撐改善了傾斜電性分解面的反演效果;楊長福等(2005)對該方法作了闡述，并與國際流行的二維反演方法作了對比研究。Smith等(1991)開發(fā)的RRI算法是我國MT工作者最為熟悉的一種二維反演方法，快速是其最大的優(yōu)點(diǎn)，但是它對復(fù)雜結(jié)構(gòu)的反映往往不是很準(zhǔn)確。胡祖志等(2005)對該方法與國際上通用的二維反演方法作了對比研究。由于各種方法各有其優(yōu)缺點(diǎn)，研究兩種或兩種以上方法的結(jié)合反演是當(dāng)前反演的一個(gè)熱點(diǎn)問題，歐東新等(2005)將SBI方法的模型構(gòu)建方式和RRI方法結(jié)合起來，提出了一種新的二維快速反演方法;Zhang等(2008)就SBI法與OCCAM法結(jié)合反演做過有意義的工作。本文采用RRI法與SBI法結(jié)合反演的辦法，先通過RRI反演得到一個(gè)光滑模型，然后采用模糊聚類方法得到一個(gè)合適的反演初始模型，再進(jìn)行SBI反演，最后得到一個(gè)折中的反演結(jié)果。通過理論模型試算、對比分析及實(shí)測資料反演對比，表明了該方法的可行性和有效性。

1 方法原理

1.1 SBI反演基本原理

SBI反演從模型構(gòu)建上著手，把地質(zhì)結(jié)構(gòu)看成由有限數(shù)量具有均勻電導(dǎo)率的電性單元組成，而這些單元的深度在橫向上是可變的，模型參數(shù)用各單元的電阻率值及其下底邊界深度表示，對于一個(gè)包含i層每層含n個(gè)最小單元的半空間模型，模型參數(shù)m可描述為:

式中ρi為第i層的電阻率，din為第i層第n個(gè)單元下底界面的深度，由于各參數(shù)的變化范圍較大，為了反演的穩(wěn)定性，采用它們的對數(shù)值進(jìn)行反演。

反演目標(biāo)函數(shù)為

通過求導(dǎo)得到模型參數(shù)的迭代系列為

對于給定的初始模型，可以通過這個(gè)迭代系列尋找極小可能構(gòu)造模型(deGroot-Hedlin et al.，2004;歐東新等，2005;葉益信等，2009)。

1.2 RRI反演基本原理

RRI反演通過解與一維相近的反演問題，來計(jì)算在每個(gè)測點(diǎn)位置下面的電阻率擾動(dòng)，把二維反演問題轉(zhuǎn)化為一系列一維反演問題。

其中δdxy和δdyx分別為TE和TM模式下觀測數(shù)據(jù)與理論數(shù)據(jù)的差值，σ0(z)為模型改變前的電導(dǎo)率值，Hy，0(yi，0)，Ex，0(yi，0)，Hx，0(yi，0)和 Ey，0(yi，0)分別是模型改變前第i個(gè)測點(diǎn)下地表的磁場值和電場值，Ey，0(yi，z)和 Ex，0(yi，z)是初始模型或者本次迭代前模型在第i個(gè)測點(diǎn)下某深度的理論電場。

綜合考慮模型橫向和垂向的不均勻性，構(gòu)造如下目標(biāo)函數(shù):

這是一個(gè)在各測點(diǎn)上的標(biāo)度Laplace范數(shù)。式中，函數(shù)f(z)可以控制標(biāo)度尺的長度，是用來度量不同深度的構(gòu)造，取f(z)=ln(z+z0)，常數(shù)z0通常是選取模型表層電阻率值和最高頻率情況下的趨膚深度。m=ln(σ)= －ln(ρ)。g(z)是起控制水平方向構(gòu)造的懲罰因子(高才坤等，2009)。

1.3 RRI與 SBI結(jié)合反演

從上面原理可以看到，由于兩種反演的模型參數(shù)不同，所以不能把兩種方法以加權(quán)的形式放在一起加以計(jì)算。因此我們先對數(shù)據(jù)作RRI反演，從半空間模型開始，得到一個(gè)較光滑的模型結(jié)果，記該結(jié)果為m1。將m1作為SBI反演的初始模型進(jìn)行SBI反演，得到反演結(jié)果m2，將m2作為最終反演結(jié)果。由于m1是一個(gè)以不同測點(diǎn)不同深度電阻率表示的結(jié)果，而SBI反演的初始模型為塊狀結(jié)構(gòu)模型，因此需要把m1轉(zhuǎn)化為塊狀結(jié)構(gòu)模型。本文采用模糊聚類方法將m1轉(zhuǎn)化為塊狀結(jié)構(gòu)模型。模糊聚類方法的原理如下:

將數(shù)據(jù)集分成L組，則可把數(shù)據(jù)集表示成:Xij(i=1，…，N;j=1，…，M)，在構(gòu)建反演初始模型時(shí)，N表示電阻率的個(gè)數(shù)，M通常等于3(電阻率的值、y、z的坐標(biāo))，L表示要構(gòu)建初始模型的層數(shù);那么數(shù)據(jù)集第j個(gè)矢量在第i類的成員函數(shù)可表示為:U=［Uij］LXN(U也稱分割矩陣);將數(shù)據(jù)集模糊聚類成L類的函數(shù)滿足以下條件(George et al.，2004):

然后通過解不含約束條件的目標(biāo)函數(shù)(Salski，2007):

式中，V=［v1，…，vL］(1 ＜ L ＜ M)為類的中心位置矢量，P為控制錄屬度權(quán)重的因子，‖·‖為范數(shù)，解上述目標(biāo)函數(shù)使其達(dá)到最小，可得類的中心位置及各自成員函數(shù)的表達(dá)式。(8)式中的X矢量表示RRI反演結(jié)果的模型參數(shù)，V矢量表示SBI反演的初始模型參數(shù)。

2 模型試驗(yàn)

2.1 模型Ⅰ試驗(yàn)

模型Ⅰ為一孤立體模型(圖1a)，一半空間存在一個(gè)異常體，電阻率為10 Ω·m，尺寸為1.6 km×0.8 km，頂面埋深800 m，半空間背景電阻率為100 Ω·m。測點(diǎn)個(gè)數(shù)為11，采集36個(gè)頻點(diǎn)(1 441，1 024，721，512，360，256，180，128，90，64，45，32，22.5，16，11.3，8，5.6，4 2，1.4，1，0.7，0.5，0.32，0.25，0.176，0.125，0.088，0.06，0.044，0.03，0.022，0.015，0.011，0.007，0.005)Hz，有限元正演網(wǎng)格為58×52。先用RRI法進(jìn)行反演(圖1b);然后用SBI法進(jìn)行反演(圖1c);最后采用RRI與SBI結(jié)合的反演(圖1d)。從反演結(jié)果圖可看出，RRI反演結(jié)果對異常體的結(jié)構(gòu)反映模糊，沒有反映出異常體的確切電阻率和邊界位置，SBI反演結(jié)果能反映出異常體的確切電阻率和邊界位置，但是其邊界輪廓有很多毛刺現(xiàn)象，不夠光滑，而RRI與SBI結(jié)合的反演結(jié)果不但可反映出異常體的確切電阻率和邊界位置，且其邊界輪廓更加光滑，與模型真實(shí)電阻率和邊界位置更加接近。從它們的迭代擬合差曲線變化(圖2)也可看出，RRI反演收斂速度較慢且其擬合差較大，SBI反演次之，采用RRI與SBI結(jié)合的反演收斂速度更快，其擬合差與SBI反演相當(dāng)。

2.2 模型Ⅱ試驗(yàn)

圖1 孤立體模型反演結(jié)果比較Fig.1 Single blocked model and inversion results

模型Ⅱ?yàn)橐恍ㄐ误w模型(圖3a)，包括四個(gè)電性單元，分別為:第一電性單元100 Ω·m，位于地表層;第二電性單元1 000 Ω·m，位于中間層左側(cè);第三電性單元10 Ω·m，位于中間層右側(cè);第四電性單元100 Ω·m，位于底層。正演計(jì)算時(shí)的頻點(diǎn)和網(wǎng)格與模型Ⅰ正演計(jì)算時(shí)一致，接收點(diǎn)共11個(gè)，點(diǎn)距為400 m。分別對數(shù)據(jù)進(jìn)行 RRI、SBI及SBI與RRI結(jié)合的反演，從它們的反演結(jié)果(圖3b，c，d)可以看出，RRI反演結(jié)果只能反映各個(gè)電性單元的大致電阻率值及電性分界面的大致形態(tài)，且對深部電性單元結(jié)構(gòu)反映不準(zhǔn);SBI反演結(jié)果能基本反映各個(gè)單元的電阻率及其電性分界面，但是對模型邊界位置的擬合稍差;兩種方法的結(jié)合反演不僅可以反映各個(gè)單元的電阻率，而且對模型邊界位置的擬合更好。它們的擬合差變化曲線變化與模型Ⅰ試驗(yàn)結(jié)果具有相似的結(jié)論。

圖2 迭代擬合差變化曲線圖Fig.2 Misfit error curves with iterations of different inversion approaches

圖3 楔形體模型反演結(jié)果比較Fig.3 Wedge model and inversion results

3 實(shí)測數(shù)據(jù)反演應(yīng)用

最后，采用該方法對一條實(shí)測MT數(shù)據(jù)進(jìn)行了反演，實(shí)測數(shù)據(jù)為武漢周邊地區(qū)地?zé)峥辈斓拇蟮仉姶艛?shù)據(jù)，該數(shù)據(jù)包含測點(diǎn)19個(gè)，頻點(diǎn)個(gè)數(shù)為34個(gè)，范圍從320 Hz到0.035 Hz，點(diǎn)距保持在200 m左右，總剖面長3.62 km。在反演開始前，首先對獲得的數(shù)據(jù)進(jìn)行二維濾波處理。先用RRI法對該數(shù)據(jù)進(jìn)行了反演(圖4)，然后再進(jìn)行SBI與RRI結(jié)合的反演，兩種反演結(jié)果的擬合差分別為20 rms和14.7 rms，表明結(jié)合反演的擬合差較小。通過結(jié)合反演結(jié)果與RRI反演結(jié)果的對比可看出，結(jié)合反演結(jié)果電性分離更明顯，電性界面位置更清晰。從圖4a可看出，該地區(qū)可分為四個(gè)電性單元:第(1)電性單元電阻率為50左右，分布于地表，厚度約0.1～0.2 km;第(2)電性單元電阻率1 000以上，分布于深度0.1～1.3 km之間，中間厚兩邊薄;第(3)電性單元電阻率為77左右，分布于橫向0～1.5 km深度0.5 ～1.5 km 范圍和橫向3 ～3.6 km 深度0.7～1.7 km范圍;第(4)電性單元電阻率為790左右，分布于深度1.5 km以下范圍。根據(jù)區(qū)域地質(zhì)資料，第(1)電性單元對應(yīng)為第四系覆蓋層;第(2)電性單元上部有三疊系的地層，下部主要為二疊系的下統(tǒng)棲霞組、孤峰組及二疊系上統(tǒng)龍?zhí)督M、大隆組地層，巖性對應(yīng)為硅質(zhì)頁巖、粉砂巖、炭質(zhì)頁巖、炭質(zhì)灰?guī)r，電阻率相對呈高阻反應(yīng);第(3)電性單元對應(yīng)為三疊系的大冶組和嘉陵江組，巖性對應(yīng)為白云質(zhì)灰?guī)r、泥質(zhì)灰?guī)r、砂質(zhì)頁巖，電阻率相對呈中低阻反應(yīng);第(4)電性單元主要為石炭系黃龍組、泥盆系上統(tǒng)五通組、志留系中統(tǒng)墳頭組，巖性對應(yīng)為生物碎屑灰?guī)r、白云質(zhì)灰?guī)r、泥質(zhì)灰?guī)r、石英砂巖夾粘土巖、石英質(zhì)礫巖、粉砂質(zhì)泥巖、細(xì)砂巖、石英砂巖，電阻率相對呈中高阻反應(yīng)。通過將反演結(jié)果與地質(zhì)資料比對，表明反演結(jié)果與地質(zhì)資料吻合的較好。

圖4 實(shí)測數(shù)據(jù)反演結(jié)果對比Fig.4 Real data inversion results

4 結(jié)論

(1)通過引入模糊聚類算法，將SBI和RRI反演結(jié)合起來反演，彌補(bǔ)了各自的局限性，同時(shí)發(fā)揮了各自的優(yōu)點(diǎn)。

(2)通過理論模型試算表明，與各自單獨(dú)反演相比，SBI和RRI的結(jié)合反演收斂速度更快，反演精度更高，對電性分界面特征反映效果更好。

(3)對武漢某測區(qū)一條實(shí)測MT資料的處理對比分析表明，這種方法的具有一定的實(shí)用性。

胡祖志，胡祥云，吳文鵬，等.2005.大地電磁二維反演方法對比研究［J］.煤田地質(zhì)與勘探，33(1):64-68.

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

歐東新，王家林.2005.二維塊狀結(jié)構(gòu)大地電磁快速反演［J］.石油物探，44(5):525-528.

吳信民，楊海燕，楊亞新，等.2013.論電法勘探的理論探測深度［J］.東華理工大學(xué)學(xué)報(bào):自然科學(xué)版，36(1):60-64.

楊長福，徐世浙.2005.國外大地電磁研究現(xiàn)狀［J］.物探與化探，29(3):243-247.

葉益信，胡祥云，金鋼燮，等.2009.大地電磁二維陡邊界反演應(yīng)用效果分析［J］.地球物理學(xué)進(jìn)展，24(2):668-674.

Constable S C，Parker R L，Constable C G.1987.Occam’s inversion:A practical algorithm for generating smooth models from electromagnetic sounding data［J］.Geophysics，52(3):289-300.

deGroot-Hedlin C，Constable S C.1990.Occam’s inversion to generate smooth two-dimensional models from magnetotelluric data［J］.Geophysical，55(12):1613-1624.

deGroot-Hedlin C，Constable S C.2004.Inversion of magnetotelluric data for 2D structure with sharp resistivity contrasts［J］.Geophysics，69(1):78-86.

George E，Tsekourasa，Haralambos Sarimveisb.2004.A new approach for measuring the validity of the fuzzy c-means algorithm［J］.Advances in Engineering Software，35:567-575.

Rodi W L，Mackie R L.2001.Nonlinear conjugate gradients algorithm for 2D magnetotelluric inversion［J］.Geophysics，66(1):174-187.

Salski A.2007.Fuzzy clustering of fuzzy ecological data［J］.Ecological Informatics，(2):262-269.

Siripunvaraporn W，Egbert G.2000.An efficient data-subspace inversion method for 2-d magnetotelluric data［J］.Geophysics，65(3):791-803.

Smith J T，booker J R.1991.Rapid inversion of two-and three-dimensional magnetotelluric data［J］.Geophys，Res.，96:3905-3922.

Smith T，Hoversten M，gasperikova E，et al.1999.Sharp boundary inversion of 2D magnetotelluric data［J］.Geophysical Prospecting，47:469-486.

Zhang L L，Yu P，Wang J L，et al.2008.Two-dimensional magnetotelluric inversion in combination of smoothest model and sharp boundary［C］//19th IAGA WG 1.2 workshop on electromagnetic induction in the earth，Beijing.

Zhang L L，Yu P，Wang J L，et al.2010.A study on 2D magnetotelluric sharp boundary inversion［J］.Geophys，(in chinese)，52(3):631-637.