利用陣列感應(yīng)測(cè)井進(jìn)行儲(chǔ)層滲透率評(píng)價(jià)

周峰，孟慶鑫，胡祥云，Evert Slob，潘和平，馬火林

1 中國(guó)地質(zhì)大學(xué)(武漢)機(jī)械與電子信息學(xué)院，武漢 4300742 中國(guó)地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室，武漢 4300743 河北地質(zhì)大學(xué)勘查技術(shù)與工程學(xué)院，石家莊 050031 4 代爾夫特理工大學(xué)土木與地球科學(xué)學(xué)院，荷蘭 2628 CN

高明亮，于生寶，鄭建波，徐暢，劉偉宇，欒卉*

吉林大學(xué)儀器科學(xué)與電氣工程學(xué)院，長(zhǎng)春 130026

?

利用陣列感應(yīng)測(cè)井進(jìn)行儲(chǔ)層滲透率評(píng)價(jià)

周峰1,2，孟慶鑫3，胡祥云2*，Evert Slob4，潘和平2，馬火林2

1 中國(guó)地質(zhì)大學(xué)(武漢)機(jī)械與電子信息學(xué)院，武漢 4300742 中國(guó)地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室，武漢 4300743 河北地質(zhì)大學(xué)勘查技術(shù)與工程學(xué)院，石家莊 050031 4 代爾夫特理工大學(xué)土木與地球科學(xué)學(xué)院，荷蘭 2628 CN

鉆井過(guò)程中儲(chǔ)層受到泥漿侵入影響的程度與儲(chǔ)層巖性有著密切關(guān)系，其中儲(chǔ)層滲透率對(duì)侵入深度有著較大影響，因此若可以獲知泥漿侵入深度，則有望對(duì)儲(chǔ)層滲透率進(jìn)行評(píng)估.本文首先建立含泥餅增長(zhǎng)的泥漿侵入數(shù)值模型，然后建立陣列感應(yīng)測(cè)井?dāng)?shù)值模型，兩者的聯(lián)合正演模擬顯示泥漿侵入對(duì)地層的影響可以反映在陣列感應(yīng)測(cè)井響應(yīng)上，利用阻尼最小二乘法對(duì)陣列感應(yīng)測(cè)井響應(yīng)進(jìn)行反演可以得到侵入深度.對(duì)侵入深度和儲(chǔ)層滲透率的關(guān)系進(jìn)行分析發(fā)現(xiàn)：在滲透率為1～100 mD(1 mD= 0.987×10-3μm2)數(shù)量級(jí)的儲(chǔ)層中，滲透率的變化可以在侵入深度上得到反映.以儲(chǔ)層和井?dāng)?shù)據(jù)進(jìn)行二維數(shù)值模擬發(fā)現(xiàn)：利用陣列感應(yīng)測(cè)井響應(yīng)反演出來(lái)的侵入深度曲線反映了滲透率在地層上的變化趨勢(shì)，采用解釋圖版的方法可以對(duì)儲(chǔ)層各層段的滲透率進(jìn)行粗略估算.

陣列感應(yīng)測(cè)井；滲透率評(píng)價(jià)；泥漿侵入

1 引言

在石油鉆井過(guò)程中，由于井底泥漿和儲(chǔ)層之間存在著壓力差，泥漿濾液滲入到儲(chǔ)層孔隙中，改變了井壁附近的流體成分，影響了測(cè)井工具對(duì)地層參數(shù)的準(zhǔn)確測(cè)量.因此，研究泥漿侵入對(duì)測(cè)井的影響并尋求有效的校正方法，對(duì)于儲(chǔ)層測(cè)井評(píng)價(jià)具有重要的意義(劉尊年等，2012).張建華等(1994)對(duì)泥漿侵入過(guò)程進(jìn)行了一維數(shù)值模擬，陳福煊和孫嘉戌(1996)建立了砂槽物理實(shí)驗(yàn)來(lái)研究泥漿侵入對(duì)地層電導(dǎo)率的影響，Dewan和Chenevert(2001)基于物理實(shí)驗(yàn)建立了泥餅動(dòng)態(tài)增長(zhǎng)的數(shù)學(xué)模型，李長(zhǎng)喜等(2006)分析了泥漿侵入中的低阻環(huán)帶的成因，常文會(huì)等(2010)在柱坐標(biāo)下對(duì)泥漿侵入進(jìn)行了二維數(shù)值模擬，Wu等(2004)和Deng等(2008)開(kāi)展了斜井泥漿侵入數(shù)值模擬，Salazar和Torres-Verdín(2009)比較了水基和油基泥漿侵入對(duì)儲(chǔ)層的影響.泥漿侵入對(duì)地層參數(shù)的影響在陣列感應(yīng)測(cè)井上有著明顯反映，仵杰等(2009)基于簡(jiǎn)化的活塞狀侵入模型進(jìn)行了感應(yīng)測(cè)井的數(shù)值模擬；鄧少貴等(2010)利用斜井泥漿侵入下的陣列側(cè)向測(cè)井來(lái)反演原狀地層電阻率；李虎等(2013)利用陣列感應(yīng)測(cè)井和自然電位測(cè)井進(jìn)行聯(lián)合反演來(lái)校正泥漿侵入影響，計(jì)算地層水電阻率.

前人所做的工作主要是將泥漿侵入視為不利因素，從而尋求校正泥漿侵入的方法，以減少其對(duì)測(cè)井產(chǎn)生的不利影響.然而受泥漿侵入影響的測(cè)井響應(yīng)也攜帶了有用的儲(chǔ)層信息，比如泥漿侵入深度與儲(chǔ)層滲透率或孔隙度及其他巖石物性參數(shù)有著密切關(guān)系.張海龍等(2005)建議利用泥漿對(duì)儲(chǔ)層流體的沖刷程度來(lái)預(yù)測(cè)油氣層產(chǎn)能，王建華等(2009)認(rèn)為儲(chǔ)層滲透率是影響鉆井濾液侵入深度的關(guān)鍵參數(shù)之一.本文認(rèn)為可以通過(guò)陣列感應(yīng)測(cè)井反演出泥漿侵入深度，再將泥漿侵入深度和儲(chǔ)層滲透率建立起關(guān)聯(lián)，從而實(shí)現(xiàn)儲(chǔ)層滲透率的半定量或者定量評(píng)價(jià).

本文建立含泥餅動(dòng)態(tài)增長(zhǎng)的泥漿侵入數(shù)值模型，基于Born幾何因子方法建立陣列感應(yīng)測(cè)井?dāng)?shù)值模型，并利用阻尼最小二乘法來(lái)反演泥漿侵入深度；分析泥漿侵入深度和儲(chǔ)層滲透率的關(guān)系，利用儲(chǔ)層數(shù)據(jù)進(jìn)行泥漿侵入和陣列感應(yīng)測(cè)井的正反演，以驗(yàn)證借助于泥漿侵入效應(yīng)的陣列感應(yīng)測(cè)井反演可以對(duì)儲(chǔ)層滲透率進(jìn)行有效評(píng)估.

2 泥漿侵入模型

除了儲(chǔ)層參數(shù)外，泥餅參數(shù)(如泥餅滲透率)對(duì)于侵入流速也有較大的影響(王建華等，2009)，因此要研究泥漿侵入和儲(chǔ)層巖性參數(shù)的關(guān)系，必須考慮泥餅在侵入中的作用.目前國(guó)內(nèi)的泥漿侵入模型一般都忽略了泥餅，或者將泥餅參數(shù)簡(jiǎn)單視為一個(gè)定值，這樣得到的侵入剖面特征雖然和真實(shí)情況相符，但是侵入液量以及侵入深度卻和真實(shí)情況偏差較大.為此將泥餅動(dòng)態(tài)增長(zhǎng)模型耦合到多相流多組分模型中，發(fā)展出一個(gè)更接近實(shí)際的泥漿侵入模型.

2.1 多相流多組分模型

以水基泥漿侵入為例，忽略氣相的參與，假設(shè)侵入過(guò)程符合等溫達(dá)西滲流，則近井區(qū)域中的孔隙壓力和含水飽和度可由油水兩相等溫達(dá)西滲流方程求解(Aziz and Settari,1979)：

(1)

(2)

其中ρw和ρo分別為水和油的密度(kg·m-3)，k為儲(chǔ)層滲透率(m2)，krw和kro分別為水和油的相對(duì)滲透率，g是重力加速度矢量(m·s-2)，D是儲(chǔ)層深度(m)，μw和μo分別為水和油的黏度(Pa·s)，Pw和Po分別是水和油的壓力(Pa)，φ為孔隙度，Sw和So分別為水和油的飽和度，t為侵入時(shí)間(s).

地層水與侵入的泥漿濾液水之間由于礦化度不同而發(fā)生鹽分混溶，采用對(duì)流擴(kuò)散方程求解水的礦化度(Navarro,2007)：

(3)

在柱坐標(biāo)下對(duì)上述公式進(jìn)行有限差分的離散，壓力和礦化度采用隱式求解，飽和度采用顯式求解，得到不同侵入時(shí)間的含水飽和度和地層水礦化度的徑向分布，再通過(guò)阿爾奇公式得到地層的電阻率(Archie,1942)：

(4)

(5)

其中Rf和Rw分別為地層綜合電阻率和地層水電阻率(Ωm)，m和n分別為膠結(jié)指數(shù)和飽和度指數(shù)，α為巖性導(dǎo)電附加系數(shù)，T是溫度(℃).

2.2 泥餅動(dòng)態(tài)模型

泥漿濾液中含有一定的固體顆粒成分，在侵入發(fā)生時(shí)泥漿中的泥質(zhì)成分逐漸沉積在井壁上形成泥餅.泥餅參數(shù)對(duì)于侵入流速至關(guān)重要，泥餅的滲透率、孔隙度隨時(shí)間的變化可以表示為(Wu et al.,2005)：

(6)

(7)

其中kmc和φmc分別為泥餅滲透率和孔隙度，Pmc為泥餅內(nèi)外壓差，kmc0和φmc0為泥餅參考滲透率和參考孔隙度，由1 psi壓差下的實(shí)驗(yàn)測(cè)量值所確定，v和δ分別為壓縮指數(shù)和相滲關(guān)系因子，均為實(shí)驗(yàn)室可測(cè)定值.該公式計(jì)算時(shí)，滲透率取mD為單位，壓力以psi為單位，其他參數(shù)無(wú)量綱.

泥漿濾液瞬時(shí)侵入速率表示為(Wu et al.,2005)：

(8)

其中qmc為瞬時(shí)泥漿侵入速率(m3·s-1)，h為目標(biāo)儲(chǔ)層厚度(m)，Pm為泥漿壓力(Pa)，μmc為泥漿黏度(Pa·s)，rw為井筒半徑(m)，rmc為泥餅內(nèi)徑(m)，下標(biāo)i為離散化后的網(wǎng)格序數(shù)，i=1表示泥餅網(wǎng)格，i=N表示模型外邊界.水基泥漿侵入時(shí)只考慮附在井壁內(nèi)側(cè)的外泥餅的影響(Salazar and Torres-Verdín,2009；王建華等，2009)，因此有rw>rmc.該式子中分母的第一部分和第二部分因式分別代表地層和泥餅的流阻.泥餅厚度隨侵入時(shí)間的增長(zhǎng)值可以用(9)式求解(Wu et al.,2005)：

(9)

其中fs為泥漿中泥質(zhì)成分的體積百分比含量.

求出每個(gè)時(shí)間步的泥餅滲透率、孔隙度和厚度，將其作為多相流模型中第一個(gè)網(wǎng)格變量，可以實(shí)現(xiàn)泥餅增長(zhǎng)模型和多相流多組分模型的耦合.整個(gè)模擬區(qū)域在徑向上采取非均勻網(wǎng)格，垂向采用均勻網(wǎng)格，內(nèi)外邊界均采取定壓邊界條件，利用MATLAB程序開(kāi)發(fā)，得到一個(gè)完整的泥漿侵入數(shù)值模擬器.

2.3 模型驗(yàn)證

Salazar和Torres-Verdín(2009)基于商用的多相流多組分模擬軟件CMG-STAR進(jìn)行二次開(kāi)發(fā)，實(shí)現(xiàn)含泥餅動(dòng)態(tài)增長(zhǎng)的泥漿侵入數(shù)值模擬.下面將本文自主開(kāi)發(fā)的泥漿侵入數(shù)值模擬程序與Salazar和 Torres-Verdín(2009)的模型在相同參數(shù)設(shè)置下進(jìn)行計(jì)算對(duì)比，以驗(yàn)證本模型的有效性.基本的儲(chǔ)層和流體參數(shù)如表1所示，分別對(duì)滲透率為10、30、100 mD的儲(chǔ)層進(jìn)行侵入100 h的模擬.對(duì)比發(fā)現(xiàn)：用本模擬器得到的侵入流速和泥餅厚度隨著時(shí)間的變化曲線(圖1a，1b)和文獻(xiàn)(Salazar and Torres-Verdín,2009)中對(duì)應(yīng)的曲線(圖2a，2b)基本吻合，證實(shí)了本文自主開(kāi)發(fā)的數(shù)值模擬程序的準(zhǔn)確性.

圖1 本文中的數(shù)值模擬程序得出的侵入流速(a)和泥餅厚度(b)隨時(shí)間的變化曲線Fig.1 Invasion rate (a) and mud cake thickness (b) versus time by the autonomous program

圖2 發(fā)表文獻(xiàn)中的侵入流速(a)和泥餅厚度(b)隨時(shí)間的變化曲線(Salazar and Torres-Verdín,2009)Fig.2 Invasion rate (a) and mud cake thickness (b) versus time in the publication (Salazar and Torres-Verdín,2009)

參數(shù)值單位井徑0.1454m儲(chǔ)層厚度1m儲(chǔ)層外徑610m徑向網(wǎng)格數(shù)51儲(chǔ)層初始?jí)毫?5166kPa井底壓力27580kPa束縛水飽和度8%殘余油飽和度10%地層水礦化度160×103mg/L泥漿礦化度3×103mg/L地層初始含水飽和度30%地層孔隙度25%泥餅參考滲透率0.03mD泥餅參考孔隙度25%泥餅最大厚度1.02×10-2m參數(shù)值單位泥餅壓縮指數(shù)0.4泥餅乘積因子0.1泥漿固態(tài)體積含量6%地層溫度100.6℃水相相對(duì)滲透率經(jīng)驗(yàn)指數(shù)3油相相對(duì)滲透率經(jīng)驗(yàn)指數(shù)3水相相對(duì)滲透率端點(diǎn)值30%油相相對(duì)滲透率端點(diǎn)值100%孔隙分布經(jīng)驗(yàn)指數(shù)15毛細(xì)壓力系數(shù)2.74×10-2Pa·m徑向彌散度0m膠結(jié)指數(shù)2飽和度指數(shù)2巖性導(dǎo)電附加系數(shù)1

3 陣列感應(yīng)測(cè)井正反演

以哈利伯頓高分辨率陣列感應(yīng)測(cè)井儀(HRAI)為參考，線圈系裝置提供探測(cè)深度分別為10 in、20 in、30 in、60 in、90 in(1 in=25.4 mm)的視電阻率曲線.考慮到小極距線圈系裝置反映的是淺層信號(hào)并主要用于校準(zhǔn)分辨率，因此在反演中舍棄了10 in探測(cè)深度的視電阻率曲線.

3.1 正演模型

本研究中陣列感應(yīng)測(cè)井正演基于幾何因子原理.該原理可直觀快速地模擬感應(yīng)測(cè)井響應(yīng)，并發(fā)展出適用于不同地電模型的修正幾何因子，其中Born幾何因子可描述非均勻地層中電導(dǎo)率在背景上的擾動(dòng)，適用于泥漿侵入地層，其視電導(dǎo)率算式可寫(xiě)為(王衛(wèi)等，2011)：

(10)

gBorn(r,z)=gDoll(1-ikrT)(1-ikrR)eikrR+ikrT，

(11)

其中σa為雙線圈系所得視電導(dǎo)率(S·m-1)，根據(jù)不同極距裝置組合可得各子陣列結(jié)果；gBorn(r,z)表示位于垂向z和徑向r處的單元環(huán)的Born幾何因子，gDoll為Doll幾何因子，k為復(fù)波數(shù)，rT、rR分別為發(fā)射線圈和接收線圈到單元環(huán)的距離(m)，i表示虛部；σ(r,z)為垂向z和徑向深度r處的單元環(huán)電導(dǎo)率(S·m-1).根據(jù)徑向分層模型可進(jìn)一步寫(xiě)出總式：

(12)

其中σt為原狀地層電導(dǎo)率(S·m-1)，σi為侵入帶電導(dǎo)率(S·m-1)，σm為井眼泥漿電導(dǎo)率(S·m-1).Gm、Gi、Gt為積分幾何因子，是由幾何因子理論得到的不同地層區(qū)域?qū)y(cè)井響應(yīng)的貢獻(xiàn)值.計(jì)算上述式子，可得各陣列感應(yīng)測(cè)井系的視電導(dǎo)率(取倒數(shù)即得視電阻率).

3.2 反演方法

基于上述正演算法，應(yīng)用N個(gè)測(cè)井?dāng)?shù)據(jù)觀測(cè)值求取具有w個(gè)參數(shù)的正演模型來(lái)完成擬合反演，該過(guò)程的數(shù)學(xué)表達(dá)式寫(xiě)為(李虎等，2013)：

dn=Fn(p1,p2,p3,…,pw),(n=1,2,…,N)

(13)

其中dn表示測(cè)井觀測(cè)數(shù)據(jù)，F(xiàn)n表示正演模型，pw表示待反演的w個(gè)模型參數(shù).采用阻尼最小二乘法求解，目標(biāo)函數(shù)可寫(xiě)為：

(14)

其中P=(p1,p2,p3,…,pw)T，表示待反演參數(shù)的轉(zhuǎn)置矩陣，P為極小解的必要條件需滿足O(P)在P處的梯度為零.將非線性函數(shù)F(P)在初始近似值P(0)處做泰勒展開(kāi)，并略去二次項(xiàng)及二次以上項(xiàng)進(jìn)行近似表示，目標(biāo)函數(shù)的矩陣形式可寫(xiě)為：

F[ΔP(0)]=[J0ΔP(0)-ε0]T[J0ΔP(0)-ε0]，

(15)

其中J0為初始剛度陣，ε0為初始?xì)埩?

若反演剛度陣J是滿秩，則由式JTJΔP=JTε可知，F(xiàn)(P)極小點(diǎn)的進(jìn)一步近似P(1)由(16)式確定：

(16)

(17)

式中P(k)為第k步迭代解，Jk為第k步反演剛度陣(Jacobi陣)，λk為阻尼因子，εk為每步殘量.

當(dāng)泥漿侵入是典型的低阻侵入時(shí)，沖洗帶和低阻環(huán)帶差別不明顯，可將其一起視為侵入帶作為一個(gè)整體來(lái)討論，故以三參數(shù)(侵入半徑、侵入層電阻率、原狀地層電阻率)進(jìn)行反演.

3.3 方法驗(yàn)證

將本研究中的陣列感應(yīng)正反演方法在一個(gè)典型的低阻侵入儲(chǔ)層中進(jìn)行驗(yàn)證，含水飽和度、地層水礦化度和地層電阻率的徑向剖面如圖3所示.可以看出，每隔24 h侵入前緣距離井壁的距離分別為0.50 m、0.75 m、0.95 m和1.1 m(圖3a)，由于侵入液量以環(huán)柱狀展開(kāi)，因此侵入深度的增加值是逐漸減緩的.地層水礦化度隨時(shí)間的變化趨勢(shì)與含水飽和度基本一致，但略滯后于含水飽和度的變化(圖3b)，分析認(rèn)為主要是由于鹽分的分子混溶效應(yīng)引起.由圖3c可以看出：地層電阻率在井壁附近存在一個(gè)10～15 Ωm的低阻區(qū)域，呈現(xiàn)出低侵的特征，對(duì)應(yīng)于泥漿沖洗帶；最外側(cè)47 Ωm為原狀地層電阻率；在沖洗帶與原狀地層之間存在5～10 Ωm的低阻環(huán)帶.低阻環(huán)帶主要是由于礦化度變化滯后于含水飽和度引起，低阻環(huán)帶和原狀地層之間的分界面對(duì)應(yīng)于含水飽和度前緣，而低阻環(huán)帶與沖洗帶之間的分界面則對(duì)應(yīng)于礦化度前緣，因此若在電測(cè)井中顯示出兩個(gè)電阻率突變的響應(yīng)特征，則第二個(gè)突變對(duì)應(yīng)于實(shí)際侵入前緣，這種影響在高阻侵入時(shí)應(yīng)尤為注意.在低阻侵入的時(shí)候，低阻環(huán)帶和沖洗帶之間的電性差別不十分明顯，因此在陣列感應(yīng)反演中可將低阻環(huán)帶和沖洗帶作為一個(gè)整體來(lái)考慮，只尋找侵入帶和原狀地層的分界面，進(jìn)而獲取侵入深度.

采用陣列感應(yīng)測(cè)井?dāng)?shù)值模型對(duì)上述地電模型進(jìn)行正演模擬，得到的陣列感應(yīng)響應(yīng)如圖4所示.可以看出：受低阻侵入的影響，侵入剖面上的視電阻率曲線偏小，分析原因是由于侵入帶電阻率較小，鄰近收發(fā)裝置增大了其信號(hào)貢獻(xiàn)率；隨著侵入半徑的增加，各子陣列視電阻率值逐漸降低，特別是深部徑向視電阻率值(R60和R90)的降低幅度較大，同時(shí)各子陣列感應(yīng)曲線受到逐步增大的低阻侵入帶的影響，分異性降低(視電阻率值趨向于低阻).陣列感應(yīng)測(cè)井響應(yīng)的變化能夠反映出泥漿侵入深度隨著時(shí)間的增加.

圖3 含水飽和度(a)、地層水礦化度(b)和綜合電阻率(c)的徑向分布Fig.3 Radial distributions of water saturation (a),water salinity (b),and comprehensive resistivity (c)

圖4 侵入24 h間隔對(duì)應(yīng)的陣列感應(yīng)測(cè)井視電阻率Fig.4 Apparent resistivities of array induction logging every 24 hours of invasion

圖5 反演得到的侵入帶電阻率、原狀地層電阻率和侵入深度Fig.5 Invasion zone resistivity,virgin zone resistivity,and invasion depth inversed from apparent resistivity

基于上述陣列感應(yīng)測(cè)井響應(yīng)進(jìn)行三參數(shù)反演，得到侵入帶電阻率、原狀地層電阻率和侵入深度的值如圖5所示，與泥漿侵入數(shù)值模擬進(jìn)行對(duì)比可以看出反演結(jié)果較為準(zhǔn)確地反映出受侵入影響儲(chǔ)層的電性結(jié)構(gòu)(圖3c)，對(duì)于本文所關(guān)心的侵入深度值而言，反演結(jié)果與數(shù)值模擬結(jié)果一致(圖3a)，證實(shí)了該反演方法的有效性.

4 滲透率評(píng)價(jià)應(yīng)用

前人進(jìn)行陣列感應(yīng)測(cè)井反演時(shí)，往往只關(guān)注如何校正出原狀地層的真電阻率，而未對(duì)侵入深度進(jìn)行利用.上面方法反演出泥漿侵入深度，其潛在的應(yīng)用是進(jìn)行儲(chǔ)層滲透率評(píng)價(jià).

4.1 滲透率與侵入深度關(guān)系

實(shí)際鉆井過(guò)程中，鉆頭以較高的剪切速度鉆開(kāi)地層，泥漿濾液以較大的速率滲入地層，此時(shí)雖然有部分固體顆粒沉積在井壁上，但由于泥漿隨鉆頭在井內(nèi)高速旋動(dòng)，使得泥餅厚度小、滲透率大，故泥餅對(duì)侵入的阻礙作用不大.起鉆后，井下液體呈現(xiàn)靜濾失狀態(tài)，大量泥漿固體顆粒沉積在井壁上，泥餅厚度增加，在內(nèi)外靜壓差作用下泥餅滲透率變小，泥餅流阻加大，泥漿濾液的濾失速率降低(Peng and Peden,1992).基于上述2.3節(jié)算例的參數(shù)進(jìn)行數(shù)值模擬，設(shè)鉆井時(shí)間為12 h，鉆井過(guò)程中泥餅厚度為0.1 mm，泥餅滲透率為0.1 mD，泥餅孔隙度為45% (Navarro,2007).起鉆后，儲(chǔ)層受靜態(tài)泥漿浸泡，泥餅孔滲參數(shù)變化按照公式(6)和(7)計(jì)算.通過(guò)對(duì)泥漿侵入進(jìn)行參數(shù)敏感性分析發(fā)現(xiàn)：隨著儲(chǔ)層滲透率增大，侵入流速也會(huì)增加，但是在相同的侵入液量下，孔隙度的增加會(huì)減少侵入深度(范宜仁等,2013).實(shí)際儲(chǔ)層的滲透率和孔隙度彼此具有相關(guān)性，因此在進(jìn)行敏感性分析時(shí)必須考慮兩者的同時(shí)變化對(duì)侵入深度的影響.以文獻(xiàn)(Salazar and Torres-Verdín,2009)給出的幾種典型的巖石物性參數(shù)為參考，設(shè)置儲(chǔ)層孔隙度分別是15%、20%、25%、28%、30%，對(duì)應(yīng)的儲(chǔ)層滲透率分別是3、10、30、100、300 mD，其他參數(shù)不變，進(jìn)行侵入100 h的數(shù)值模擬.

圖6a、6b分別記錄了侵入流速和累計(jì)侵入液量隨時(shí)間的變化曲線，可以看出：在鉆井的12 h時(shí)間，由于泥餅流阻較小，泥漿濾液以較大的速率侵入到儲(chǔ)層中，孔滲性較好的儲(chǔ)層對(duì)應(yīng)的侵入流速較大(圖6a)，侵入總液量也較大(圖6b)，此時(shí)儲(chǔ)層滲透率是決定侵入流速和液量的主要因素.隨著鉆井完成，井下泥漿呈現(xiàn)靜濾失的狀態(tài)，泥餅逐漸變厚且滲透率降低，侵入流速逐漸降至最小，且在不同孔滲性儲(chǔ)層中侵入流速差別不大(圖6a)，隨著侵入時(shí)間的增加侵入總液量不再呈現(xiàn)明顯增長(zhǎng)(圖6b)，分析認(rèn)為主要原因是靜濾失過(guò)程中泥餅滲透性小且厚度大，使得泥餅對(duì)侵入過(guò)程的影響顯著，儲(chǔ)層本身的滲透率不再是主要敏感性參數(shù)，這種趨勢(shì)在高滲透率儲(chǔ)層中特別明顯.

鑒于鉆井結(jié)束后侵入液量變化不大，因此提取侵入24 h時(shí)不同孔滲性儲(chǔ)層的侵入剖面進(jìn)行分析(圖7a)，得到儲(chǔ)層滲透率和侵入深度的關(guān)系如圖7b所示.對(duì)比可以看出：隨著儲(chǔ)層滲透率增加，侵入深度逐漸增加，但是兩者的關(guān)系呈現(xiàn)冪指數(shù)的遞增關(guān)系，即儲(chǔ)層滲透率增大，侵入深度相應(yīng)增加，但是其增加幅度越來(lái)越小.分析原因認(rèn)為：一方面儲(chǔ)層滲透率越大，泥餅上的壓降越大，儲(chǔ)層內(nèi)部壓力梯度越小(圖8)，因此儲(chǔ)層滲透率變化對(duì)侵入流速的影響越小；另一方面，儲(chǔ)層滲透率越大，儲(chǔ)層對(duì)應(yīng)的孔隙度越大，侵入液量雖然越大，但是由此帶來(lái)的侵入深度的增加值卻被孔隙度抵消了一部分.由此也可以推斷，在高滲透率儲(chǔ)層(103mD的量級(jí)以上)，隨著儲(chǔ)層滲透率的增加，侵入深度不再明顯增加.

以上分析可以看出：儲(chǔ)層滲透率在一定的范圍內(nèi)(1 mD至100 mD數(shù)量級(jí))可以由侵入深度所反映.在小于1 mD的儲(chǔ)層中，達(dá)西滲流理論受限，因此這里不做討論.

4.2 半定量預(yù)測(cè)

以華北某油田測(cè)井和取心數(shù)據(jù)為例進(jìn)行泥漿侵入的數(shù)值模擬.取深度從1036 m到1096 m的井段數(shù)據(jù)，井眼直徑為0.2 m，鉆井液礦化度為12000 mg/L，泥漿密度1130 kg·m-3，水和油的黏度分別為0.968 Pa·s和2.99 Pa·s，儲(chǔ)層和泥漿壓差為2 MPa，重力加速度為9.8 m·s-2，滲透率和孔隙度變化分別如圖9a、9b所示，這里的滲透率為水平滲透率，垂向滲透率為水平滲透率的1/10.

以該井段的數(shù)據(jù)進(jìn)行泥漿侵入二維數(shù)值模擬，得到侵入24 h后的含水飽和度剖面(圖10a)和電阻率剖面(圖10b).

將陣列感應(yīng)測(cè)井響應(yīng)進(jìn)行逐層反演，得到侵入深度曲線如圖11所示.比較侵入深度(圖11紅線)和儲(chǔ)層滲透率(圖11藍(lán)線)，可以看出反演出來(lái)的侵入深度和儲(chǔ)層滲透率在縱向上的變化趨勢(shì)基本一致，尤其是侵入深度曲線變化明顯反映出了滲透率曲線中兩個(gè)極低和一個(gè)極高的特征.

圖6 不同孔滲性儲(chǔ)層對(duì)應(yīng)的侵入流速(a)和侵入液量(b)隨時(shí)間的變化Fig.6 Invasion rate and volume vs time for different porosity and permeability

圖7 侵入24 h含水飽和度徑向分布(a)和儲(chǔ)層滲透率與侵入深度關(guān)系曲線(b)Fig.7 Radial distributions of water saturation (a),and invasion depths vs reservoir permeabilities (b) after 24 hours of invasion

圖8 侵入24 h時(shí)不同孔滲性儲(chǔ)層對(duì)應(yīng)的壓力徑向分布Fig.8 Radial distributions of formation pressure after 24 hours of invasion

圖9 華北某井段的滲透率(a)和孔隙度(b)曲線Fig.9 Permeability (a) and porosity (b) curves of an oil well in North China

圖10 侵入24 h的含水飽和度(a)和電阻率(b)剖面Fig.10 Profiles of water saturation (a) and resistivity (b) after 24 hours of invasion

圖11 反演的侵入深度(紅色)和儲(chǔ)層滲透率(藍(lán)色)對(duì)比Fig.11 Contrast between the curves of invasion depth (red) and permeability (blue)

圖12 二維軸對(duì)稱的分層地層模型Fig.12 Two-dimensional axisymmetric layered form ation model

圖13 侵入24 h的含水飽和度分布Fig.13 Water saturation distribution after 24 hours of invasion

圖14 反演的侵入深度曲線Fig.14 Inversed invasion depth curve

4.3 圖版定量評(píng)估

關(guān)于泥漿侵入的參數(shù)敏感性分析認(rèn)為：儲(chǔ)層的滲透率、孔隙度和含水飽和度是對(duì)侵入深度影響較大的地層參數(shù)，這幾個(gè)參數(shù)也是在同一口井中隨著深度變化較大的(Zhou et al.,2015).在實(shí)際生產(chǎn)中可以制作一系列解釋圖版集，建立起不同孔隙度、滲透率和含水飽和度對(duì)應(yīng)的侵入深度值曲線，這樣就可以利用獲取的侵入深度及其他可測(cè)定參數(shù)值來(lái)粗略評(píng)估未知的儲(chǔ)層滲透率.

在進(jìn)行評(píng)估之前，將測(cè)井、鉆井、取心中可以獲知的基本參數(shù)代入到泥漿侵入數(shù)值模型中，并設(shè)置不同的孔隙度、滲透率和含水飽和度值進(jìn)行一系列數(shù)值模擬，從而得到不同“孔-滲-飽”參數(shù)對(duì)應(yīng)的侵入深度值，制作成侵入深度圖版集.在現(xiàn)場(chǎng)測(cè)量時(shí)，孔隙度和含水飽和度通過(guò)常規(guī)測(cè)井獲取后找到其對(duì)應(yīng)的圖版曲線，再利用陣列感應(yīng)測(cè)井?dāng)?shù)據(jù)反演出侵入深度值，最后找到相應(yīng)曲線中該侵入深度值所對(duì)應(yīng)的儲(chǔ)層滲透率值.

以上述井?dāng)?shù)據(jù)為基礎(chǔ)，將地層“孔-滲-飽”隨地層深度的變化簡(jiǎn)化為一個(gè)概念化的層狀模型(圖12)，其中孔隙度和含水飽和度是常規(guī)測(cè)井方法可測(cè)的，忽略了重力和垂向滲透率影響，數(shù)值模擬泥漿侵入24 h后的徑向侵入剖面如圖13所示.對(duì)陣列感應(yīng)測(cè)井響應(yīng)進(jìn)行反演，計(jì)算得到各層的侵入深度分別為：0.37 m、0.78 m和1.07 m，如圖14所示.

基于這口井的數(shù)據(jù)制作了一組反映侵入深度和滲透率關(guān)系的圖版集，這里挑選出含水飽和度為20%、30%和40%的圖版，如圖15a—15c所示.針對(duì)第一層地層，找到在含水飽和度為20%的圖版中孔隙度為10%的曲線，由侵入深度為0.37 m的橫坐標(biāo)找到縱坐標(biāo)對(duì)應(yīng)的值，即為估算的儲(chǔ)層滲透率(如圖中虛線所示).同理可以得到其他兩個(gè)地層的滲透率估算值.表2的對(duì)比看出：估算的滲透率和預(yù)設(shè)的滲透率誤差在一個(gè)數(shù)量級(jí)內(nèi)，從儲(chǔ)層評(píng)估的尺度上來(lái)說(shuō)誤差在可以接受的范圍.誤差來(lái)源一方面是陣列感應(yīng)測(cè)井低的徑向空間分辨率，另一方面是漸變而非活塞狀的侵入前緣.可見(jiàn)使用解釋圖版可以粗略地估算出儲(chǔ)層各個(gè)層段的滲透率，然而這種方法的缺點(diǎn)是需要預(yù)先進(jìn)行大量的數(shù)值模擬計(jì)算來(lái)獲取足夠多的參數(shù)圖版，以對(duì)應(yīng)實(shí)際儲(chǔ)層參數(shù)較大的取值范圍.

表2 滲透率估算值和預(yù)設(shè)值(單位：mD)

5 結(jié)論

本文建立了一個(gè)含泥餅增長(zhǎng)的泥漿侵入數(shù)值模型，基于Born幾何因子方法建立了陣列感應(yīng)測(cè)井?dāng)?shù)值模型，基于阻尼最小二乘法對(duì)陣列感應(yīng)測(cè)井響應(yīng)進(jìn)行反演得到泥漿侵入深度值；分析了儲(chǔ)層滲透率對(duì)侵入深度的影響，并利用陣列感應(yīng)測(cè)井對(duì)儲(chǔ)層滲透率進(jìn)行評(píng)估.得出結(jié)論如下：在滲透率為1 ～100 mD數(shù)量級(jí)的儲(chǔ)層中，泥漿侵入深度與儲(chǔ)層滲透率有著較好的對(duì)應(yīng)關(guān)系，利用陣列感應(yīng)測(cè)井反演出的侵入深度值可以反映出儲(chǔ)層滲透率的變化趨勢(shì)，基于此方法建立的解釋圖版集可以在現(xiàn)場(chǎng)生產(chǎn)中用來(lái)粗略估算儲(chǔ)層滲透率.

圖15 含水飽和度為20%(a)、30%(b)和40%(c)時(shí)的滲透率解釋圖版Fig.15 Interpretation charts for water saturations of 20% (a),30% (b) and 40% (c)

盡管目前還沒(méi)有一個(gè)理想的數(shù)學(xué)模型將泥漿侵入深度和儲(chǔ)層滲透率之間的相關(guān)性進(jìn)行直接和準(zhǔn)確的量化，但是本文的研究為利用陣列感應(yīng)測(cè)井進(jìn)行滲透率定量評(píng)估提供可行性論證，也為將來(lái)借助陣列感應(yīng)測(cè)井?dāng)?shù)據(jù)對(duì)儲(chǔ)層滲透率進(jìn)行直接準(zhǔn)確的反演提供理論依據(jù).

Archie G E.1942.The electrical resistivity log as an aid in determining some reservoir characteristics.Transactions of the American Institute of Mining,Metallurgical and Petroleum Engineering,146(1):54-62,doi:10.2118/942054-G.

Aziz K,Settari A.1979.Petroleum Reservoir Simulation.London:Applied Science Publishers.

Chang H W,Pan H P,Zhou F.2010.Two-dimensional numerical simulation of mud Invasion.Earth Science-Journal of China University of Geosciences (in Chinese),35(4):674-680,doi:10.3799/dqkx.2010.082.

Chen F X,Sun J X.1996.Simulation test of radial electric conductivity during the mud filtrate invading porous formation.Chinese Journal of Geophysics (Acta Geophysica Sinica) (in Chinese),39(S1):371-378.

Deng S G,Fan Y R,Li G X,et al.2008.Numerical simulation of mud invasion in deviated wells in curvilinear coordinates.∥2nd IASME/WSEAS International Conference on Geology and Seismology.Cambridge,UK,140-148.

Deng S G,Li Z Q,Fan Y R,et al.2010.Numerical simulation of mud invasion and its array laterolog response in deviated wells.Chinese Journal of Geophysics (in Chinese),53(4):994-1000,doi:10.3969/j.issn.0001-5733.2010.04.024.

Dewan J T,Chenevert M E.2001.A model for filtration of water-base mud during drilling:determination of mudcake parameters.Petrophysics,42(3):237-250.

Fan Y R,Hu Y Y,Li H,et al.2013.Numerical simulation of mud-cake dynamic formation and reservoir mud filtrate invasion.Well Logging Technology (in Chinese),37(5):466-471,doi:10.3969/j.issn.1004-1338.2013.05.002.

Li C X,Ouyang J,Zhou C C,et al.2006.Forming mechanism and application of low resistivity annulus in oil reservoirs invaded by fresh drilling mud.Petroleum Exploration and Development (in Chinese),32(6):82-86,doi:10.3321/j.issn:1000-0747.2005.06.020.

Li H,Fan Y R,Hu Y Y,et al.2013.Joint inversion of HDIL and SP with a five-parameter model for estimation of connate water resistivity.Chinese Journal of Geophysics (in Chinese),56(2):688-695,doi:10.6038/cjg20130233.

Liu Z N,Sun J M,Chi X R,et al.2012.Analysis of research present situation of mud invasion.Progress in Geophysics (in Chinese),27(6):2594-2601,doi:10.6038/j.issn.1004-2903.2012.06.037.

Navarro D.2007.Effects of invasion transient on resistivity time-lapsed logging [Master′s thesis].Houston:Houston University.

Peng S J,Peden J M.1992.Prediction of filtration under dynamic conditions.∥SPE International Symposium on Formation Damage Control.Lafayette,Louisiana:SPE:503-511,doi:10.2118/23824-MS.

Salazar J M,Torres-Verdín C.2009.Quantitative comparison of processes of oil-and water-based mud-filtrate invasion and corresponding effects on borehole resistivity measurements.Geophysics,74(1):E57-E73,doi:10.1190/1.3033214.

Wang J H,Yan J N,Zheng M.2009.Prediction model for invasion radius of solids and filtrate in drilling fluids.Acta Petrolei Sinica (in Chinese),30(6):923-926,doi:10.7623/syxb200906023.

Wang W,Liao D L,Wu J.2011.Inversion of array induction based on Born geometric factor.Journal of Oil and Gas Technology (in Chinese),33(8):82-85，doi：10.3969/j.issn.1000-9752.2011.08.018.

Wu J H,Torres-Verdin C,Sepehrnoori K,et al.2004.Numerical simulation of mud-filtrate invasion in deviated wells.SPE Reservoir Evaluation &Engineering,7(2):143-154,doi:10.2118/87919-PA.

Wu J H,Torres-Verdín C,Sepehrnoori K,et al.2005.The influence of water-base mud properties and petrophysical parameters on mudcake growth,filtrate invasion,and formation pressure.Petrophysics,46(1):14-32.

Wu J,Feng J,Xie X C,et al.2009.Forward analysis of high definition induction logging response in mud invasion formation.Well Logging Technology (in Chinese),33(3):212-217,doi:10.3969/j.issn.1004-1338.2009.03.004.

Zhang H L,Liu G Q,Zhou C C,et al.2005.Reservoir productivity prediction by array induction logging data.Petroleum Exploration and Development (in Chinese),32(3):84-87,doi:10.3321/j.issn:1000-0747.2005.03.021.

Zhang J H,Hu Q,Liu Z H.1994.A theoretical model for mud-filtrate invasion in reservoir formations during drilling.Acta Petrolei Sinica (in Chinese),15(4):73-78,doi:10.7623/syxb199404010.

Zhou F,Hu X Y,Meng Q X,et al.2015.Model and method of permeability evaluation based on mud invasion effects.Applied Geophysics,12(4):482-492,doi:10.1007/s11770-015-0513-1.

附中文參考文獻(xiàn)

常文會(huì),潘和平,周峰.2010.泥漿侵入二維數(shù)值模擬.地球科學(xué):中國(guó)地質(zhì)大學(xué)學(xué)報(bào),35(4):674-680,doi:10.3799/dqkx.2010.082.

陳福煊,孫嘉戌.1996.泥漿濾液侵入孔隙地層徑向?qū)щ娞匦缘哪M實(shí)驗(yàn).地球物理學(xué)報(bào),39(S1):371-378.

鄧少貴,李智強(qiáng),范宜仁等.2010.斜井泥漿侵入仿真及其陣列側(cè)向測(cè)井響應(yīng)數(shù)值模擬.地球物理學(xué)報(bào),53(4):994-1000,doi:10.3969/j.issn.0001-5733.2010.04.024.

范宜仁,胡云云,李虎等.2013.泥餅動(dòng)態(tài)生長(zhǎng)與泥漿侵入模擬研究.測(cè)井技術(shù),37(5):466-471,doi:10.3969/j.issn.1004-1338.2013.05.002.

李長(zhǎng)喜,歐陽(yáng)健,周燦燦等.2006.淡水鉆井液侵入油層形成低電阻率環(huán)帶的綜合研究與應(yīng)用分析.石油勘探與開(kāi)發(fā),32(6):82-86,doi:10.3321/j.issn:1000-0747.2005.06.020.

李虎,范宜仁,胡云云等.2013.基于陣列感應(yīng)與自然電位聯(lián)合反演地層水電阻率.地球物理學(xué)報(bào),56(2):688-695,doi:10.6038/cjg20130233.

劉尊年,孫建孟,遲秀榮等.2012.泥漿侵入研究現(xiàn)狀分析.地球物理學(xué)進(jìn)展,27(6):2594-2601,doi:10.6038/j.issn.1004-2903.2012.06.037.王建華,鄢捷年,鄭曼.2009.鉆井液固相和濾液侵入儲(chǔ)層深度的預(yù)測(cè)模型.石油學(xué)報(bào),30(6):923-926,doi:10.7623/syxb200906023.王衛(wèi)，廖東良，仵杰.2011.基于Born幾何因子的陣列感應(yīng)反演研究.石油天然氣學(xué)報(bào)，33(8):82-85，doi：10.3969/j.issn.1000-9752.2011.08.018.

仵杰,馮娟,解茜草等.2009.泥漿侵入地層中高分辨率感應(yīng)測(cè)井響應(yīng)特征的正演分析.測(cè)井技術(shù),33(3):212-217,doi:10.3969/j.issn.1004-1338.2009.03.004.

張海龍,劉國(guó)強(qiáng),周燦燦等.2005.基于陣列感應(yīng)測(cè)井資料的油氣層產(chǎn)能預(yù)測(cè).石油勘探與開(kāi)發(fā),32(3):84-87,doi:10.3321/j.issn:1000-0747.2005.03.021.

張建華,胡啟,劉振華.1994.鉆井泥漿濾液侵入儲(chǔ)集層的理論計(jì)算模型.石油學(xué)報(bào),15(4):73-78,doi:10.7623/syxb199404010.

(本文編輯 何燕)

高明亮，于生寶，鄭建波等.2016.基于IGA算法的電阻率神經(jīng)網(wǎng)絡(luò)反演成像研究.地球物理學(xué)報(bào)，59(11):4372-4382,doi:10.6038/cjg20161136.

Gao M L，Yu S B，Zheng J B，et al.2016.Research of resistivity imaging using neural network based on immune genetic algorithm.Chinese J.Geophys.(in Chinese),59(11):4372-4382,doi:10.6038/cjg20161136.

高明亮，于生寶，鄭建波，徐暢，劉偉宇，欒卉*

吉林大學(xué)儀器科學(xué)與電氣工程學(xué)院，長(zhǎng)春 130026

摘要 為滿足地球物理資料反演解釋的高精度、快速、穩(wěn)定的要求，本文結(jié)合免疫遺傳算法尋優(yōu)速度快和BP神經(jīng)網(wǎng)絡(luò)反演不依賴初始模型等優(yōu)點(diǎn)，設(shè)計(jì)了一種將BP神經(jīng)網(wǎng)絡(luò)和免疫遺傳算法進(jìn)行有機(jī)結(jié)合的全局優(yōu)化反演策略，并將該策略成功地應(yīng)用于二維高密度電法數(shù)據(jù)反演.利用免疫遺傳算法(Immune Genetic Algorithm，簡(jiǎn)稱IGA)對(duì)神經(jīng)網(wǎng)絡(luò)的反演參數(shù)進(jìn)行同步優(yōu)化，提高了電阻率反演的精度.仿真和實(shí)驗(yàn)結(jié)果驗(yàn)證設(shè)計(jì)的全局優(yōu)化反演策略取得了較好的效果，通過(guò)與線性反演方法和BP法以及遺傳神經(jīng)網(wǎng)絡(luò)法等反演方法進(jìn)行比較，得出該方法具有反演精度更高，反演時(shí)間更短等顯著優(yōu)勢(shì)的結(jié)論.

關(guān)鍵詞 免疫遺傳算法；BP神經(jīng)網(wǎng)絡(luò)；高密度電阻率法；反演精度

基金項(xiàng)目 地面電磁探測(cè)(SEP)系統(tǒng)研制-野外試驗(yàn)研究(201311193-05(SinoProbe-09-02-05))資助

作者簡(jiǎn)介 高明亮，男，1987年生，博士，研究方向?yàn)橹绷麟姺ǚ蔷€性反演方法研究.E-mail：mlgao13@mails.jlu.edu.cn

*通訊作者 欒卉，女，1979年生，副教授，碩士生導(dǎo)師，研究方向?yàn)殡姶趴碧椒椒ㄑ芯?E-mail：luanhui@jlu.edu.cn

Abstract In order to meet the requirements of high precision,high speed and stability in geophysical inversion,we design a global optimization inversion strategy based on BP neural network and immune genetic algorithm,considering the fast searching speed of immune genetic algorithm and the independence of the initial model in BP neural network inversion.The strategy is successfully applied to the two-dimensional high density resistivity inversion.By using Immune Genetic Algorithm (IGA) for synchronous optimization of the neural network inversion parameters,the precision of the resistivity inversion can be improved .The results of experiment and simulation verify that the global optimization strategy achieves great results.Comparing with the linear inversion method,BP method,and genetic neural network method,our method has higher precision and shorter time of inversion.

Keywords Immune Genetic Algorithm;BP neural network;High density resistivity method;Inversion precision

1 引言

隨著計(jì)算機(jī)技術(shù)和最優(yōu)化方法的發(fā)展，非線性反演方法在地球物理領(lǐng)域逐步顯示出強(qiáng)大的生命力(徐海浪和吳小平，2006；劉斌，2010；張凌云，2011；戴前偉等，2014).近年來(lái)，國(guó)內(nèi)外關(guān)于電阻率非線性反演方法進(jìn)行了較深入的研究，取得了一些成果.目前非線性反演方法研究主要集中在一維、二維反演方面，其中神經(jīng)網(wǎng)絡(luò)算法以其較強(qiáng)的計(jì)算能力和學(xué)習(xí)能力在地球物理反演領(lǐng)域應(yīng)用最為廣泛(Poulton and El Fouly，1991；Calderon-Macis et al.，2000；Spichak et al.，2002；Fani et al.，2013).Jimmy(2004)、Singh(2005)等率先將神經(jīng)網(wǎng)絡(luò)算法引入電阻率測(cè)深的一維反演工作中.在2D、3D電阻率非線性反演中，采用神經(jīng)網(wǎng)絡(luò)法進(jìn)行了電阻率反演研究，不依賴于初始模型，不需要求解偏導(dǎo)數(shù)矩陣且反演速度較快，取得了較好的反演效果(El-Qady and Keisuke，2001；Mann，2006；Ahmad and Samsudin，2009；Singh et al.，2010).但是由于神經(jīng)網(wǎng)絡(luò)的固有特性，存在網(wǎng)絡(luò)結(jié)構(gòu)不穩(wěn)定，反演容易陷入局部極小、收斂慢等問(wèn)題(Poulton et al.，1992；Cranganu et al.，2007；Vladimir and Krasnopolsky，2007).國(guó)內(nèi)外學(xué)者對(duì)此進(jìn)行了較為深入的研究，對(duì)神經(jīng)網(wǎng)絡(luò)算法做了進(jìn)一步改進(jìn)，通過(guò)蟻群算法、遺傳算法等算法對(duì)網(wǎng)絡(luò)進(jìn)行優(yōu)化(Dai et al.，2014；Marina et al.，2014)，克服了神經(jīng)網(wǎng)絡(luò)的上述缺點(diǎn).Maiti等(2011)提出將蒙特卡羅法與神經(jīng)網(wǎng)絡(luò)相結(jié)合進(jìn)行混合二維電阻率反演，張凌云和劉鴻福(2011)利用蟻群算法對(duì)神經(jīng)網(wǎng)絡(luò)的結(jié)構(gòu)參數(shù)進(jìn)行了優(yōu)化，改善了神經(jīng)網(wǎng)絡(luò)容易陷入局部極小的問(wèn)題，但是網(wǎng)絡(luò)收斂速度慢，反演精度有待提高.通過(guò)對(duì)蟻群、模擬退火等多種非線性方法研究發(fā)現(xiàn)遺傳算法與神經(jīng)網(wǎng)絡(luò)相結(jié)合進(jìn)行混合反演，其反演精度最高(張凌云，2011).遺傳算法屬于隨機(jī)優(yōu)化算法，通過(guò)初始值的改變尋找最優(yōu)解，而在網(wǎng)絡(luò)維數(shù)較大的情況下，隨機(jī)優(yōu)化算法較組合排序算法具有更明顯的優(yōu)勢(shì).雖然遺傳算法與神經(jīng)網(wǎng)絡(luò)混合反演的反演精度最高，但是遺傳算法優(yōu)化神經(jīng)網(wǎng)絡(luò)參數(shù)時(shí)間較長(zhǎng)，遺傳算法的解空間容易被相似性高的適應(yīng)值個(gè)體充滿而缺乏多樣性，導(dǎo)致算法早熟，局部搜索能力變差(周明等，1999；單文桃等，2013).針對(duì)遺傳算法的不足，本文嘗試將免疫遺傳算法與神經(jīng)網(wǎng)絡(luò)進(jìn)行混合反演.通過(guò)在遺傳算法中引入免疫算子并將BP算法嵌入其中，優(yōu)化神經(jīng)網(wǎng)絡(luò)的反演參數(shù).這種混合反演方法解決了單一遺傳算法優(yōu)化神經(jīng)網(wǎng)絡(luò)時(shí)間過(guò)長(zhǎng)的問(wèn)題，加快了遺傳算法搜索速度，提高全局搜索能力，確?？焖偈諗坑谌肿顑?yōu)解.同時(shí)克服了BP神經(jīng)網(wǎng)絡(luò)容易陷入局部極小，容錯(cuò)能力差，網(wǎng)絡(luò)結(jié)構(gòu)不穩(wěn)定等問(wèn)題，減少了反演迭代次數(shù)，節(jié)約了運(yùn)算時(shí)間，并且可以獲得更好的電阻率反演結(jié)果，實(shí)現(xiàn)了快速、高精度反演.

2 電阻率二維免疫遺傳算法與BP神經(jīng)網(wǎng)絡(luò)混合反演

由于直流電阻率的建模與模型參數(shù)優(yōu)化直接影響反演成像質(zhì)量，因此對(duì)反演模型中的視電阻率進(jìn)行預(yù)處理建模顯得尤為重要.為滿足地球物理反演高精度、快速、穩(wěn)定的要求，并且適應(yīng)復(fù)雜地質(zhì)勘查情況，提高反演成像的快速性和準(zhǔn)確性，本文應(yīng)用BP神經(jīng)網(wǎng)絡(luò)算法對(duì)直流電阻率數(shù)據(jù)進(jìn)行建模，利用免疫遺傳算法對(duì)神經(jīng)網(wǎng)絡(luò)參數(shù)進(jìn)行同步優(yōu)化，模型輸入為X、Z位置和視電阻率值D矩陣，U為D歸一化處理的數(shù)據(jù)，yd和ya分別是參考模型和實(shí)際模型輸出的真實(shí)電阻率矩陣，參考模型空間可通過(guò)地質(zhì)信息或廣義線性反演結(jié)果給定模型參數(shù)的變化范圍.E為參考模型輸出與實(shí)際輸出電阻率值的均方差矩陣，BP神經(jīng)網(wǎng)絡(luò)反演參數(shù)通過(guò)免疫遺傳算法進(jìn)行在線和離線優(yōu)化獲得，從而使直流電阻率模型的實(shí)際電阻率能夠準(zhǔn)確跟蹤參考模型輸出，也就是均方差值E趨近于0.參考模型和BP神經(jīng)網(wǎng)絡(luò)模型都需要用BP神經(jīng)網(wǎng)絡(luò).其中，參考模型是用來(lái)訓(xùn)練BP神經(jīng)網(wǎng)絡(luò)的初始反演參數(shù)，它是以離線的方式通過(guò)有限差分法正演數(shù)值計(jì)算產(chǎn)生訓(xùn)練集進(jìn)而訓(xùn)練BP神經(jīng)網(wǎng)絡(luò)從而得到初始的反演參數(shù)(戴前偉等，2014).而作為反演模型的BP網(wǎng)絡(luò)，其作用則是用以控制被控對(duì)象的輸出能夠跟隨參考模型的輸出.

基于IGA算法的電阻率神經(jīng)網(wǎng)絡(luò)反演成像主要分為三部分：BP神經(jīng)網(wǎng)絡(luò)建模、IGA算法優(yōu)化神經(jīng)網(wǎng)絡(luò)反演參數(shù)、神經(jīng)網(wǎng)絡(luò)反演成像.圖1為IGA優(yōu)化的神經(jīng)網(wǎng)絡(luò)反演成像模型系統(tǒng)結(jié)構(gòu)圖.

2.1 BP神經(jīng)網(wǎng)絡(luò)建模

BP神經(jīng)網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)如圖2a所示.神經(jīng)網(wǎng)絡(luò)包括輸入層、隱含層和輸出層.其中，輸入層神經(jīng)元的數(shù)量等于輸入變量的個(gè)數(shù)，輸出層神經(jīng)元的數(shù)量等于輸出變量的個(gè)數(shù)，隱含層神經(jīng)元的數(shù)量可以依據(jù)輸入信息量的大小和復(fù)雜度作相應(yīng)的改變.本文將X、Z位置和視電阻率作為神經(jīng)網(wǎng)絡(luò)的輸入，以真實(shí)電阻率作為神經(jīng)網(wǎng)絡(luò)的輸出，這樣的數(shù)據(jù)類型最適合神經(jīng)網(wǎng)絡(luò)的訓(xùn)練和測(cè)試.具體討論參見(jiàn)第3節(jié).BP神經(jīng)網(wǎng)絡(luò)的主要思想就是通過(guò)樣本集訓(xùn)練網(wǎng)絡(luò)，使得神經(jīng)網(wǎng)絡(luò)的輸出與期望輸出的誤差平方和達(dá)到最小.它是通過(guò)反復(fù)迭代使相對(duì)誤差函數(shù)在斜率下降的方向上計(jì)算網(wǎng)絡(luò)權(quán)閾值和偏差的變化而逐漸達(dá)到目標(biāo).每一次權(quán)值和偏差的變化都與網(wǎng)絡(luò)誤差的響應(yīng)成正比，并以發(fā)射傳播的方式傳到每一層.

圖1 IGA優(yōu)化的反演成像模型系統(tǒng)結(jié)構(gòu)圖Fig.1 The system structure diagram of inversion imaging model

圖2 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)示意圖(a) 神經(jīng)網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu);(b) 神經(jīng)網(wǎng)絡(luò)算式.Fig.2 A neural network structure diagram(a) The neural network topology;(b) The expression of neural network.

BP神經(jīng)網(wǎng)絡(luò)具體的結(jié)構(gòu)算式如圖2b所示.其中神經(jīng)網(wǎng)絡(luò)的每一隱含層都需要用到激活函數(shù)，本文選擇logsig函數(shù)，其函數(shù)表達(dá)式如式(1)所示：

(1)

輸入層到隱含層單元J表達(dá)式如下：

(2)

其中xi是神經(jīng)網(wǎng)絡(luò)的輸入變量，netj是隱含層第j個(gè)神經(jīng)元的輸出，wij是神經(jīng)網(wǎng)絡(luò)的連接權(quán)值，bj是單元j第j個(gè)神經(jīng)元的閾值.

隱含層單元J到輸出層aj表達(dá)式為

(3)

BP神經(jīng)網(wǎng)絡(luò)訓(xùn)練算法一般采用梯度下降法來(lái)計(jì)算神經(jīng)網(wǎng)絡(luò)權(quán)值，權(quán)值更新規(guī)則如下：

(4)

式中，ε為學(xué)習(xí)動(dòng)量，η為學(xué)習(xí)率，t為訓(xùn)練次數(shù).

BP神經(jīng)網(wǎng)絡(luò)的學(xué)習(xí)誤差E為

(5)

最終BP神經(jīng)網(wǎng)絡(luò)通過(guò)不斷地迭代直到誤差函數(shù)E足夠小，網(wǎng)絡(luò)訓(xùn)練結(jié)束.式中P為樣本總數(shù)，N為輸出節(jié)點(diǎn)數(shù)，apj和dpj分別為第p個(gè)樣本中輸出層神經(jīng)網(wǎng)絡(luò)第j個(gè)神經(jīng)元實(shí)際輸出和期望輸出值.

2.2 免疫遺傳算法

免疫遺傳算法就是將免疫理論(Immune Algorithm，IA)和基本遺傳算法(Simple Genetic Algorithm，SGA)各自的優(yōu)點(diǎn)結(jié)合起來(lái)的一個(gè)多學(xué)科相互交叉、滲透的優(yōu)化算法.免疫遺傳算法可以看作一種新型融合算法，羅小平(2002)在遺傳算法中加入免疫算子，使遺傳算法變成具有免疫功能的新算法，是一種改進(jìn)的遺傳算法，該算法既保留了遺傳算法的搜索特性，又利用了免疫算法的多機(jī)制求解多目標(biāo)函數(shù)最優(yōu)解的自適應(yīng)特性，在很大程度上避免了“早熟”，收斂于局部極值.

2.3 IGA算法優(yōu)化神經(jīng)網(wǎng)絡(luò)反演參數(shù)

本文應(yīng)用免疫遺傳算法優(yōu)化神經(jīng)網(wǎng)絡(luò)的反演參數(shù).對(duì)神經(jīng)網(wǎng)絡(luò)反演參數(shù)進(jìn)行優(yōu)化設(shè)計(jì)時(shí)，應(yīng)在保證神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)穩(wěn)定的前提下，力求加快網(wǎng)絡(luò)收斂速度，以實(shí)現(xiàn)高精度反演.在優(yōu)化設(shè)計(jì)時(shí)，本文使用的目標(biāo)函數(shù)與約束函數(shù)均是高度非線性的，屬于多目標(biāo)優(yōu)化問(wèn)題，適合于免疫遺傳算法的求解.

2.3.1 建立目標(biāo)函數(shù)

對(duì)于神經(jīng)網(wǎng)絡(luò)的優(yōu)化設(shè)計(jì)可建立多種目標(biāo)函數(shù)，但通常多以目標(biāo)最小為目標(biāo)函數(shù)，本文設(shè)計(jì)的目標(biāo)函數(shù)為

(6)

將公式(1)、(2)、(3)代入(6)可以得到：

(7)

2.3.2 優(yōu)化變量設(shè)計(jì)

式(7)為IGA優(yōu)化的目標(biāo)函數(shù)，其中目標(biāo)函數(shù)內(nèi)變量wij、bpj為免疫遺傳算法優(yōu)化的變量X.根據(jù)相關(guān)文獻(xiàn)資料以及實(shí)驗(yàn)測(cè)試，本文確定選擇有一個(gè)隱含層的三層IGA-BP神經(jīng)網(wǎng)絡(luò)進(jìn)行訓(xùn)練.免疫遺傳算法優(yōu)化神經(jīng)網(wǎng)絡(luò)步驟如下：

(1) 編碼

將待尋優(yōu)變量X作為抗體.標(biāo)準(zhǔn)遺傳算法一般采用二進(jìn)制編碼，其搜索能力強(qiáng)，交叉變異計(jì)算方便，但需要頻繁的編碼與解碼，計(jì)算量較大而且精度不高.由于十進(jìn)制編碼精度更高，搜索空間更大，待優(yōu)化參量較多，因而本研究采用十進(jìn)制編碼，可以減小搜索空間，有利于加快收斂速度.

(2) 初始抗體的生成

隨機(jī)產(chǎn)生種群，為了避免無(wú)效抗體的產(chǎn)生，本文將變量的取值控制在規(guī)定范圍內(nèi).

(3) 疫苗的提取

在免疫遺傳算法中，疫苗是指從具體的待求解問(wèn)題的先驗(yàn)知識(shí)中提取的一種特征信息.免疫遺傳算法的運(yùn)行效率和性能很大程度上取決于疫苗的合理選擇.優(yōu)良的疫苗對(duì)免疫算法的收斂速度有著決定性的影響.本文將每一代的最優(yōu)解作為疫苗，動(dòng)態(tài)地建立疫苗庫(kù)，當(dāng)前的最優(yōu)解比疫苗庫(kù)中的最差疫苗的親和力高時(shí)，則取代該最差疫苗.

(4) 抗體產(chǎn)生(交叉與變異)

用設(shè)定的交叉概率Pc和變異概率Pm，按照標(biāo)準(zhǔn)遺傳算法進(jìn)行交叉和變異操作.

(5) 接種疫苗

從種群中選擇要進(jìn)行接種的個(gè)體，按照動(dòng)態(tài)計(jì)算的接種率將疫苗的基因片段依次接入，即用疫苗基因置換抗體相應(yīng)基因位上原有的基因，形成更好的免疫種群.

(6) 計(jì)算適應(yīng)度和濃度

本文選擇種群大小為100.運(yùn)用IGA來(lái)搜索最優(yōu)的函數(shù)參數(shù)X，計(jì)算抗體的適應(yīng)度.公式(7)為本文參數(shù)尋優(yōu)目標(biāo)函數(shù)，而免疫遺傳算法中的適應(yīng)度函數(shù)通常為目標(biāo)最大化問(wèn)題，因此，適應(yīng)度函數(shù)Fitness(x)為

(8)

(7) 免疫機(jī)制選擇

免疫機(jī)制選擇指每次遺傳操作后，隨機(jī)抽取一些個(gè)體注射疫苗，然后進(jìn)行免疫檢測(cè).若親和度提高，則繼續(xù)；反之，說(shuō)明在交叉、變異的過(guò)程中出現(xiàn)了嚴(yán)重的退化現(xiàn)象，這時(shí)該個(gè)體將被父代中對(duì)應(yīng)的個(gè)體所取代.

(8) 終止條件的判斷

免疫遺傳算法在神經(jīng)網(wǎng)絡(luò)的優(yōu)化過(guò)程中，本文定義終止條件為三種形式:

(a) 達(dá)到設(shè)定的循環(huán)迭代次數(shù)；

(b) 連續(xù)幾代優(yōu)化解沒(méi)有改善；

(c) 達(dá)到最大優(yōu)化時(shí)間.

滿足其一，則將與抗原親和度最高的抗體加入免疫記憶數(shù)據(jù)庫(kù)中，輸出最優(yōu)結(jié)果并進(jìn)行下一步.

(9) 獲取神經(jīng)網(wǎng)絡(luò)參數(shù)

通過(guò)免疫遺傳算法的深度優(yōu)化，提供了一組最佳的網(wǎng)絡(luò)參數(shù)，然后再進(jìn)行神經(jīng)網(wǎng)絡(luò)反演，若滿足反演效果，則進(jìn)行下一步，否則，繼續(xù)進(jìn)行訓(xùn)練，直到達(dá)到滿意結(jié)果.

(10) 仿真測(cè)試，結(jié)果輸出.免疫遺傳算法優(yōu)化神經(jīng)網(wǎng)絡(luò)框圖如圖3所示.

3 二維IGA-BP電阻率反演成像建模

電阻率的反演建模是一個(gè)極其復(fù)雜的問(wèn)題，Singh等(2010)使用采集的視電阻率和基單元作為網(wǎng)絡(luò)模型的輸入節(jié)點(diǎn)，真實(shí)的電阻率作為模型的輸出節(jié)點(diǎn)，通過(guò)改變異常體的大小和位置作為網(wǎng)絡(luò)訓(xùn)練樣本，對(duì)于復(fù)雜的地質(zhì)結(jié)構(gòu)，則需要建立大量訓(xùn)練樣本，耗費(fèi)運(yùn)算時(shí)間.徐海浪等(2006)考慮上述缺點(diǎn)，使用視電阻率作為網(wǎng)絡(luò)模型的輸入節(jié)點(diǎn)，真實(shí)電阻率作為網(wǎng)絡(luò)模型的輸出.但是由于輸入和輸出節(jié)點(diǎn)過(guò)大，導(dǎo)致網(wǎng)絡(luò)收斂速度降低，訓(xùn)練和測(cè)試需要大量的時(shí)間.針對(duì)這個(gè)問(wèn)題，戴前偉等(2014)提出了新的建模方式，將采集并預(yù)處理后的視電阻率數(shù)據(jù)按分布區(qū)域分組，以每個(gè)分組數(shù)據(jù)為神經(jīng)網(wǎng)絡(luò)的一個(gè)訓(xùn)練樣本.雖然通過(guò)分割樣本的大小來(lái)控制輸入輸出節(jié)點(diǎn)的數(shù)量，提高了網(wǎng)絡(luò)收斂的速度并且節(jié)約了運(yùn)算時(shí)間，但是依然存在輸入輸出節(jié)點(diǎn)過(guò)多的問(wèn)題.Ahmand等(2010)使用采集的視電阻率和相應(yīng)的水平、垂直位置作為網(wǎng)絡(luò)模型的輸入節(jié)點(diǎn)，相應(yīng)位置的真實(shí)電阻率作為輸出節(jié)點(diǎn)，其特點(diǎn)是輸入輸出節(jié)點(diǎn)少，網(wǎng)絡(luò)結(jié)構(gòu)簡(jiǎn)單，網(wǎng)絡(luò)收斂速度最佳，但是要求建立足夠多訓(xùn)練網(wǎng)絡(luò)的樣本.

圖3 免疫遺傳神經(jīng)網(wǎng)絡(luò)框圖Fig.3 The flowchart of Immune Genetic Algorithm

對(duì)于IGA-BP神經(jīng)網(wǎng)絡(luò)，以上四種建模方法均有不足，Singh等(2010)采用的建模方法網(wǎng)絡(luò)結(jié)構(gòu)復(fù)雜，訓(xùn)練樣本數(shù)目過(guò)多，網(wǎng)絡(luò)收斂速度慢；徐海浪等(2006)采用的建模方法輸入輸出節(jié)點(diǎn)過(guò)多，網(wǎng)絡(luò)結(jié)構(gòu)臃腫，參數(shù)過(guò)多，不便于IGA優(yōu)化；戴前偉等(2014)采用的樣本分割建模方法，雖然是輸入輸出節(jié)點(diǎn)有所減少，網(wǎng)絡(luò)收斂速度加快，但是依然存在節(jié)點(diǎn)過(guò)多的問(wèn)題.對(duì)于IGA優(yōu)化設(shè)計(jì)來(lái)講，參數(shù)過(guò)多容易導(dǎo)致優(yōu)化時(shí)間較長(zhǎng)；Ahmand等(2010a,2010b)提出了很好的建模方法，輸入輸出節(jié)點(diǎn)少，網(wǎng)絡(luò)結(jié)構(gòu)精煉，但是單個(gè)樣本集數(shù)量過(guò)大.

本文充分考慮以上多種建模的優(yōu)缺點(diǎn)，采用如下的建模方法:首先，確定模型的輸入輸出節(jié)點(diǎn)，對(duì)采集的視電阻率進(jìn)行預(yù)處理，然后將預(yù)處理后的視電阻率和相應(yīng)的水平、垂直位置作為網(wǎng)絡(luò)模型輸入節(jié)點(diǎn)，真實(shí)的電阻率作為網(wǎng)絡(luò)模型的輸出節(jié)點(diǎn)，這樣構(gòu)建的網(wǎng)絡(luò)結(jié)構(gòu)參數(shù)最小，模型結(jié)構(gòu)精簡(jiǎn)，收斂速度快.其次分割訓(xùn)練樣本集.由于單個(gè)樣本集整體輸入網(wǎng)絡(luò)訓(xùn)練，數(shù)量過(guò)大，我們將大的樣本集分隔開(kāi)成多個(gè)小樣本集，來(lái)訓(xùn)練網(wǎng)絡(luò)，這樣在保證網(wǎng)絡(luò)收斂快的同時(shí)節(jié)約了運(yùn)算時(shí)間，構(gòu)建的模型反演效果更好.

為了驗(yàn)證不同學(xué)習(xí)算法對(duì)神經(jīng)網(wǎng)絡(luò)反演的影響，本文采用溫納-斯倫貝格裝置進(jìn)行測(cè)量，共使用37個(gè)測(cè)量電極，測(cè)量極距為1 m.共采集210個(gè)數(shù)據(jù)點(diǎn).通過(guò)上述的討論，本文IGA-BP神經(jīng)網(wǎng)絡(luò)模型使用三層神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu).網(wǎng)絡(luò)的樣本集通過(guò)改變異常體的大小、位置、形態(tài)來(lái)獲得，其中測(cè)試樣本集不參與網(wǎng)絡(luò)的訓(xùn)練，圖4為構(gòu)建的部分訓(xùn)練樣本模型.圖5為用于測(cè)試的樣本模型.

其中在構(gòu)建訓(xùn)練樣本集時(shí)要盡可能考慮反演目標(biāo)中異常體的形態(tài)、大小以及反演目標(biāo)的周圍環(huán)境測(cè)試樣本集應(yīng)該與反演目標(biāo)較為接近.本文使用RES2DMOD軟件獲得的2D數(shù)據(jù)集作為訓(xùn)練與測(cè)試數(shù)據(jù)集.

圖4 用于訓(xùn)練的部分樣本模型Fig.4 The sample models for network training

圖5 用于測(cè)試的樣本模型Fig.5 The sample models for network test

在人工神經(jīng)網(wǎng)絡(luò)中，BP神經(jīng)網(wǎng)絡(luò)占據(jù)了重要的位置，許多研究對(duì)其內(nèi)部的訓(xùn)練算法做了進(jìn)一步的改進(jìn)，從而提高收斂速度.其中比較有代表的是附加動(dòng)量法、自適應(yīng)學(xué)習(xí)速率法、彈性BP算法(簡(jiǎn)稱RPROP).彈性BP算法只取偏導(dǎo)數(shù)符號(hào)，不考慮偏導(dǎo)數(shù)的幅值，偏導(dǎo)數(shù)的符號(hào)決定權(quán)值更新的方向.該訓(xùn)練算法收斂速度較前述兩種方法快，因而在大量實(shí)際應(yīng)用中，彈性BP算法非常有效.

Ahmand(2010a,2010b)、張凌云等(2011)均采用彈性BP算法作為網(wǎng)絡(luò)訓(xùn)練算法，對(duì)于大數(shù)據(jù)來(lái)說(shuō)，該算法收斂速度最快，因而在地球物理中得到廣泛應(yīng)用.

BP算法大多是局部?jī)?yōu)化算法，這些算法本身自然在訓(xùn)練過(guò)程容易陷入局部小的問(wèn)題，即在訓(xùn)練過(guò)程中可能找不到某個(gè)具體問(wèn)題的解.針對(duì)這個(gè)問(wèn)題，張凌云等(2010)提出混合算法進(jìn)行優(yōu)化設(shè)計(jì)，并指出遺傳算法優(yōu)化神經(jīng)網(wǎng)絡(luò)的反演精度最高，但是基本的遺傳算法優(yōu)化時(shí)間過(guò)長(zhǎng)，本文在基本遺傳算法的基礎(chǔ)上，又引進(jìn)了免疫算法，通過(guò)引入免疫機(jī)制提高遺傳算法的優(yōu)化速度，使其在優(yōu)化過(guò)程中，快速得到最佳解.

本文分別用三種算法——采用彈性BP訓(xùn)練算法的神經(jīng)網(wǎng)絡(luò)算法、遺傳神經(jīng)網(wǎng)絡(luò)算法(簡(jiǎn)稱GABP法)和免疫遺傳神經(jīng)網(wǎng)絡(luò)算法(簡(jiǎn)稱IGA-BP法)對(duì)所建立的訓(xùn)練樣本集進(jìn)行訓(xùn)練和測(cè)試.實(shí)驗(yàn)中算法的運(yùn)行參數(shù)分別設(shè)置如下:種群規(guī)模均為100，疫苗接種概率:a=0.7，交叉概率:Pc=0.67，變異概率:Pm= 0.05，最大代數(shù)為1000，動(dòng)態(tài)疫苗庫(kù)的疫苗數(shù)量M=8.目標(biāo)誤差均設(shè)置為0.01，具體的參數(shù)設(shè)置如表1所示.

三種不同學(xué)習(xí)算法性能如表2、3所示，從中我們可以清晰地看出非線性算法混合反演的優(yōu)越性，反演計(jì)算時(shí)間短、精度高.本文選取50個(gè)樣本集和12個(gè)測(cè)試集進(jìn)行試驗(yàn)驗(yàn)證，其中表2性能對(duì)比為50次訓(xùn)練的平均值，學(xué)習(xí)步數(shù)和時(shí)間均是在目標(biāo)收斂的情況下.表3、4為分別使用BP法、GABP法和IGA-BP法的進(jìn)行訓(xùn)練和測(cè)試的結(jié)果對(duì)比表，由此得到以下結(jié)論：

表1 不同算法的參數(shù)設(shè)置

(1) 針對(duì)神經(jīng)網(wǎng)絡(luò)容易陷入局部小，收斂性相對(duì)較差的問(wèn)題，應(yīng)用多種算法對(duì)神經(jīng)網(wǎng)絡(luò)進(jìn)行權(quán)重優(yōu)化，從表2中，IGA-BP神經(jīng)網(wǎng)絡(luò)算法明顯優(yōu)于GABP神經(jīng)網(wǎng)絡(luò)，優(yōu)化時(shí)間明顯少于GABP神經(jīng)網(wǎng)絡(luò)，而且迭代次數(shù)要明顯少于GABP神經(jīng)網(wǎng)絡(luò)，大大節(jié)約了運(yùn)算時(shí)間.

表2 反演算法權(quán)重優(yōu)化過(guò)程對(duì)比

表3 三種反演算法性能對(duì)比

表4 三種反演算法試驗(yàn)結(jié)果對(duì)比

(2) BP神經(jīng)網(wǎng)絡(luò)法進(jìn)行反演時(shí)需要用雙隱含層才可以獲得較好的網(wǎng)絡(luò)進(jìn)行預(yù)測(cè)，需要多次調(diào)整隱含層節(jié)點(diǎn)數(shù)，因此迭代次數(shù)較多、費(fèi)時(shí)較長(zhǎng)，而IGA-BP法反演和GABP反演時(shí)，由于先優(yōu)化了BP神經(jīng)網(wǎng)絡(luò)初始的權(quán)闕值，所以在訓(xùn)練時(shí)僅單隱含層網(wǎng)絡(luò)就可得到更優(yōu)的結(jié)果，反演速度較快.

(3) 收斂比是指在50次訓(xùn)練的過(guò)程中收斂的次數(shù)占總數(shù)的比例，表明算法的優(yōu)越性.綜合對(duì)比不難發(fā)現(xiàn)IGA-BP反演收斂效果最好.

(4) 反演時(shí)間是指訓(xùn)練好神經(jīng)網(wǎng)絡(luò)的權(quán)闕值參數(shù)之后，實(shí)際反演所需的時(shí)間.IGA-BP神經(jīng)網(wǎng)絡(luò)的反演時(shí)間在相同情況之下，所需時(shí)間最短，其次為GABP神經(jīng)網(wǎng)絡(luò)，而B(niǎo)P神經(jīng)網(wǎng)絡(luò)的時(shí)間較長(zhǎng)，表明神經(jīng)泛化能力不及前兩者.

(5) 預(yù)測(cè)誤差是指參考模型輸出與實(shí)際輸出電阻率值的均方差值(MSE)，判斷系數(shù)R2通常用于判定自變量和因變量之間的符合度，通常介于0～1之間，R2越大符合性越好(張凌云，2011；戴前偉，2014)，這兩者是表明數(shù)據(jù)穩(wěn)定性的參數(shù).從三種神經(jīng)網(wǎng)絡(luò)方法的預(yù)測(cè)誤差值看，IGA-BP要優(yōu)于GABP和BP算法，在測(cè)試的樣本模型中預(yù)測(cè)最大絕對(duì)值誤差，前者僅為12.3，后兩者高達(dá)32.19和82.28.因此，通過(guò)上述五點(diǎn)的對(duì)比，可見(jiàn)IGA-BP和GABP神經(jīng)網(wǎng)絡(luò)算法無(wú)論訓(xùn)練成功率還是數(shù)據(jù)穩(wěn)定性都要明顯好于BP算法.從收斂速度和運(yùn)算時(shí)間上看，IGA-BP算法反演效果是最優(yōu)的.

表5 反演算法試驗(yàn)結(jié)果對(duì)比

4 數(shù)據(jù)仿真與模型反演

為了進(jìn)一步驗(yàn)證上述反演算法的性能，本文選擇的算例模型為經(jīng)典的垂直雙阻模型和3個(gè)異常體的組合模型，并將其用于IGA-BP和最小二乘法RES2DINV、GABP、BP反演結(jié)果的對(duì)比.

模型1采用溫納-斯倫貝格裝置，電極極距為1 m，共使用37個(gè)電極.背景場(chǎng)的電阻率為100 Ωm，其中的兩處異常均為高阻500 Ωm，地下埋深分別為1 m和3 m，在水平16 m和19 m之間，異常體的尺寸大小均為1 m×3 m，如圖6a所示.在電阻率二維神經(jīng)網(wǎng)絡(luò)反演測(cè)試中，神經(jīng)網(wǎng)絡(luò)輸入節(jié)點(diǎn)為X、Z位置和視電阻率數(shù)據(jù)，網(wǎng)絡(luò)輸出為真實(shí)電阻率參數(shù)，網(wǎng)絡(luò)輸出的反演結(jié)果如圖6所示.

從圖6和表5的反演結(jié)果分析可知，從異常體形態(tài)上來(lái)看，IGA-BP神經(jīng)網(wǎng)絡(luò)反演的結(jié)果能較為準(zhǔn)確地反應(yīng)異常體的位置和電阻率，異常體形態(tài)和輪廓清晰，最小二乘法反演結(jié)果沒(méi)有區(qū)分開(kāi)異常體之間的間隙，IGA-BP反演效果明顯優(yōu)于RES2DINV軟件的反演結(jié)果.從反演時(shí)間上來(lái)看，IGA-BP算法首先利用IGA法對(duì)神經(jīng)網(wǎng)絡(luò)進(jìn)行初始權(quán)值、闕值的最優(yōu)尋找，為神經(jīng)網(wǎng)絡(luò)算法提供較好的訓(xùn)練初始權(quán)值、闕值，優(yōu)化好神經(jīng)網(wǎng)絡(luò)參數(shù)后，神經(jīng)網(wǎng)絡(luò)可以直接反演，而且不需要重復(fù)訓(xùn)練反演時(shí)間只需要1.6 s，而最小二乘法反演時(shí)間為2.7 s.僅單純從反演時(shí)間來(lái)看，IGA-BP反演時(shí)間最短，具有較高的反演效率.

模型2為3個(gè)異常體的組合模型，主要用于檢驗(yàn)IGA-BP和GABP、BP的反演性能的差異，模型2的基本參數(shù)設(shè)置如下：依然采用溫納-斯倫貝格裝置，接收電極37個(gè)，電極極距1 m，背景場(chǎng)的電阻率為150 Ωm，高阻異常體1的形態(tài)大小為1 m×3 m，電阻率值為800 Ωm，在水平13 m和16 m之間，低阻異常體2的形態(tài)大小均為0.5 m×1 m，頂部埋深2 m，在水平14.75 m和15.25 m之間，電阻率值為20 Ωm，低阻異常體3的形態(tài)大小均為12 m×2 m，頂部埋深1 m，在水平25 m和37 m之間，電阻率值為20 Ωm，頂部埋深1 m，在電阻率二維神經(jīng)網(wǎng)絡(luò)反演測(cè)試中，神經(jīng)網(wǎng)絡(luò)輸入節(jié)點(diǎn)為 X、Z位置和視電阻率數(shù)據(jù)，網(wǎng)絡(luò)輸出為真實(shí)電阻率參數(shù)，網(wǎng)絡(luò)輸出的反演結(jié)果如圖7所示.

從表5和圖7的反演結(jié)果分析可知，從預(yù)測(cè)誤差和判斷系數(shù)來(lái)看，IGA-BP算法的MSE和R2僅為3.07%和0.97，反演性能最優(yōu).從異常體形態(tài)上來(lái)看，三種神經(jīng)網(wǎng)絡(luò)反演的結(jié)果均能較為準(zhǔn)確地反映異常體的位置和電阻率，異常體形態(tài)和輪廓都較為清晰.BP算法的反演結(jié)果在高阻異常體1、低阻異常體2處出現(xiàn)了較大失真，多了一些冗余構(gòu)造，同時(shí)低阻異常體的電阻率值與實(shí)際差異較大.GABP反演算法較BP算法得到一些改善，但在低阻異常體2依然存在較大的失真.IGA-BP算法的反演結(jié)果與實(shí)際模型更為接近，異常體的形態(tài)結(jié)構(gòu)、位置、電阻率值均得到較好的反映，證實(shí)了算法具有較好的泛化能力和穩(wěn)定性.

圖6 樣本模型1的示意圖及不同方法的反演結(jié)果(a) 模型示意圖;(b) IGA-BP反演結(jié)果;(c) RES2DIV軟件反演結(jié)果.Fig.6 The model 1 and inversion results of different algorithms(a) The Schematic model;(b) IGA-BP inversion results;(c) RES2DIV software inversion results.

圖7 樣本模型2的示意圖及不同方法的反演結(jié)果(a) 模型示意圖;(b) IGA-BP反演結(jié)果;(c) GABP反演結(jié)果;(d) BP反演結(jié)果.Fig.7 The model 2 and inversion results of different algorithms(a) Schematic model;(b) IGA-BP inversion results;(c) GABP inversion results;(d) BP inversion results.

圖8 采空區(qū)示意圖及視電阻率斷面圖(a) 野外示意圖;(b) 視電阻率斷面圖.Fig.8 The schematic of Goaf and its apparent resistivity cross-sectional view(a) The Goaf schematic view;(b)The apparent resistivity cross-sectional view.

圖9 不同方法的反演結(jié)果(a) IGA-BP反演結(jié)果;(b) BP反演結(jié)果;(c) 線性反演結(jié)果.Fig.9 The inversion results of different algorithms(a) IGA-BP inversion results;(b) BP inversion results;(c) RES2DIV software inversion results.

5 野外實(shí)驗(yàn)驗(yàn)證

為了進(jìn)一步驗(yàn)證IGA-BP算法的優(yōu)越性，本文進(jìn)行了野外工程實(shí)驗(yàn)，實(shí)驗(yàn)區(qū)域?yàn)榧质￠L(zhǎng)春市吉林大學(xué)朝陽(yáng)校區(qū)原日偽時(shí)期留下的防空洞區(qū)域，圖8a為防空洞示意圖.圖中A處為下水管道口，下水管道口徑為1 m×2 m，深度為3 m，在第12個(gè)電極和第14個(gè)電極之間.B處為防空洞區(qū)域.防空洞走向?yàn)橛蓶|向北，頂部距離地面0.5 m，深度為3 m，在第25個(gè)電極和第37個(gè)電極之間，本次探測(cè)在測(cè)區(qū)內(nèi)由南向北布設(shè)1條測(cè)線，使用GenPenE60D分布式高密度電法儀器進(jìn)行實(shí)地測(cè)量，采用溫納-斯倫貝格裝置，電極極距為1 m，共使用37個(gè)電極，共完成測(cè)點(diǎn)210個(gè).

本文用三種算法分別對(duì)采集的視電阻率進(jìn)行反演成像，視電阻率斷面圖、反演結(jié)果分別如圖8b、圖9所示.由反演結(jié)果圖9可知，從異常體形態(tài)上來(lái)看，IGA-BP算法和BP算法均較正確地反映了異常體的位置和電阻率值，而在異常體輪廓上IGA-BP算法反演結(jié)果更為清晰.而線性反演方法反演結(jié)果與實(shí)際結(jié)果相比失真很大，并且異常體的阻值也與實(shí)際值相差很大，無(wú)法觀測(cè)數(shù)據(jù)末端的異常體.基于以上判斷，不難發(fā)現(xiàn)IGA-BP算法的優(yōu)越性，從而驗(yàn)證了IGA-BP算法具有較強(qiáng)的泛化能力和穩(wěn)定性.

6 結(jié)論

本文聯(lián)合多種算法實(shí)現(xiàn)了高密度電阻率法的非線性混合反演，克服了單獨(dú)使用BP神經(jīng)網(wǎng)絡(luò)易陷入局部極小、依賴初始權(quán)閾值、反演時(shí)間長(zhǎng)及收斂速度慢等缺點(diǎn).文中提出將免疫遺傳算法與神經(jīng)網(wǎng)絡(luò)混合進(jìn)行反演，提出了最小精煉的反演模型結(jié)構(gòu)，通過(guò)模型仿真與實(shí)驗(yàn)驗(yàn)證，并與BP等反演算法對(duì)比，該算法反演誤差較小，反演時(shí)間較短，較傳統(tǒng)的線性反演成像更加清晰準(zhǔn)確，具有很好的應(yīng)用前景.總之，IGA-BP反演方法是一種新型的直流電法非線性反演成像技術(shù)，對(duì)復(fù)雜地質(zhì)構(gòu)造和三維電阻率反演成像的問(wèn)題，IGA-BP反演方法有待進(jìn)一步進(jìn)行理論研究，需在以后的工作中不斷完善.

致謝 感謝吉林大學(xué)地球信息探測(cè)儀器教育部重點(diǎn)實(shí)驗(yàn)室提供了研究平臺(tái)，感謝吉林大學(xué)時(shí)間域電磁組的全體成員給予的幫助和指導(dǎo)，感謝審稿專家對(duì)該文的審定.

References

Ahmad N，Samsudin T.2009.Using artificial neural networks to invert 2D DC resistivity imaging data for high resistivity contrast regions:A MATLAB application.Computers and Geosciences，35(11)：2268-2274.

Ahmad N，W A T.Wan.A，Samsudin T.2010a.Inversion of quasi-3D DC resistivity imaging data using artificial neural networks.Journal of Earth System Science，119(1)：27-40.

Ahmad N，W A T.Wan.A，Samsudin H.T.2010b.3D inversion of DC data using artificial neural networks.Studia Geophysica et Geodaetica，54(3)：465-485.

Calderon-Macis C，Sen I K，Stoffa P L.2000.Artificial neural networks for parameter estimation in geophysics.Geophysical Prospecting，48(1)：21-47.

Cranganu C.2007.Using artificial neural networks to predict the presence of overpressured zones in the Anadarko Basin，Oklahoma.Pure and Applied Geophysics，164(10)：2067-2081.

Dai Q W，Jiang F B，Li D.2014.Nonlinear inversion for electrical resistivity tomography based on chaotic DE-BP algorithm.Journal of Central South University.，21(5)：2018-2025.

Dai Q W，Jiang F B，Dong L．2014．RBFNN inversion for electrical resistivity tomography based on Hannan-Quinn criterion．Chinese J.Geophys.(in Chinese)，57(4)：1335-1344．

El-Qady G，Keisuke U.2001.Inversion of DC resistivity data using neural networks.Geophysical Prospecting，49(1):417-430.

Fani E A，George J T，Ioannis F G，et al.2013.Estimation of seasonal variation of ground resistance using Artificial Neural Networks.Electric Power Systems Research，94：113-121.

Jimmy Stephen，Manoj C，Singh S B.2004.A direct inversion seheme for deep resistivity Sounding data using artificial neural networks.Journal of Earth System Science，113(l):49-66.

Liu B．2010．Study on the water-bearing structure advanced detection and water inrush hazards real-time monitoring in tunnel based on the electrical resistivity method and induced polarization method [Ph．D．thesis] (in Chinese)．Shandong：Shandong University．

Luo X P．2002．The research on artificial immune genetic learning algorithm and its application in engineering [Ph．D．thesis] (in Chinese)．Hangzhou：Zhejiang University．

Poulton M，El-Fouly A.1991.Preprocessing GPR signatures for cascading neural network classification.SEG Expanded Abstracts，10(1):507-509.

Poulton M，Sternberg K，Glass C.1992.Neural network pattern recognition of subsurface EM images.Journal of Applied Geophysic，29(1):21-36.

Spichak V V，F(xiàn)ukuoka K，Kobayashi T.et al.2002.Artificial neural network reconstruction of geoelectrical parameters of the Minou fault zone by scalar CSAMT data.Journal of Applied Geophysics，49(1):75-90.

Singh U K，Tiwari R K，Singh S B.2005.One-dimensional inversion of geo-electrical resistivity sounding data using artificial neural networks—a case study.Computers &Geosciences，31(1):99-108.

Singh U K，Tiwari R K，Singh S B.2010.Inversion of 2-D DC resistivity data using rapid optimizationand minimal complexity neural network.Nonlinear Processes in Geophysics，17(1):65-76.

Mann C J H .2006.Geophysical Applications of Artificial Neural Networks and Fuzzy Logic.Kybernetes，35(1):3-4.

Shan W T，Chen X A，He Y，et al． 2013．Application of immune genetic algorithm based fuzzy RBF neural Network in high-speed motorized spindles．Chinese Journal of Mechanical Engineering (in Chinese)，49(23)：167-173．

Maiti S，Gupta G，Erram V C.2011.Inversion of Schlumberger resistivity sounding data from the critically dynamic Koyna region using the Hybrid Monte Carlo-based neural network Approach.Nonlinear Processes in Geophysics.，18(2):179-192.

Marina R C，NiklasL，Thomas K，Jasper A.2014.Vrugt Two-dimensional probabilistic inversion of plane-wave electromagnetic data:methodology，model constraints and joint inversion with electrical resistivity data.Geophysical Journal International，196 (3)：1508-1524.

Vladimir M，Krasnopolsky.2007.Neural network emulations for complex multidimensional geophysical mappings:Applications of neural network techniques to atmospheric and oceanic satellite retrievals and numerical modeling.Reviews of Geophysics，45(1):RG3009.

Xu H L,Wu X P．2006．Inversion of 2D resistivity neural networks．Chinese J.Geophys.(in Chinese),49(2)：584-589．

Zhang L Y,Liu H F．2011． The application of ABP method in high-density resistivity method inversion． Chinese J.Geophys.(in Chinese)，54(1)：227-233．

Zhang L Y．2011．The study of nonlinear inversion method in high-density resistivity method inversion[Ph．D．thesis] (in Chinese)．Taiyuan：Taiyuan University of Technology．

Zhou M，Sun S D．1999．Genetic Algorithm Theory and Application (in Chinese)．Bengjing：National Defense Industry Press．

附中文參考文獻(xiàn)

戴前偉，江沸菠，董莉.2014.基于漢南-奎因信息準(zhǔn)則的電阻率層析成像徑向基神經(jīng)網(wǎng)絡(luò)反演.地球物理學(xué)報(bào),57(4):1335-1344.

單文桃，陳小安，合燁等.2013.基于免疫遺傳算法的模糊徑向基函數(shù)神經(jīng)網(wǎng)絡(luò)在高速電主軸中的應(yīng)用.機(jī)械工程學(xué)報(bào)，49(23):167-173.

劉斌.2010.基于電阻率法與激電法的隧道含水地質(zhì)構(gòu)造超前探測(cè)與突水災(zāi)害實(shí)時(shí)監(jiān)測(cè)研究[博士論文]．山東：山東大學(xué).

羅小平.2002.人工免疫遺傳學(xué)習(xí)算法及其工程應(yīng)用研究[博士論文].杭州:浙江大學(xué).

徐海浪，吳小平.2006.電阻率二維神經(jīng)網(wǎng)絡(luò)反演.地球物理學(xué)報(bào)，49(2):584-589.

張凌云，劉鴻福.2011.ABP法在高密度電阻率法反演中的應(yīng)用．地球物理學(xué)報(bào)，54(1)：227-233.

張凌云.2011.高密度電阻率勘探反演的非線性方法研究[博士論文].太原:太原理工大學(xué).

周明，孫樹(shù)棟.1999.遺傳算法原理與應(yīng)用.北京:國(guó)防工業(yè)出版社.

(本文編輯 汪海英)

Evaluation of reservoir permeability using array induction logging

ZHOU Feng1,2,MENG Qing-Xin3,HU Xiang-Yun2*,SLOB Evert4,PAN He-Ping2,MA Huo-Lin2

1 School of Mechanics and Electronic Information,China University of Geosciences,Wuhan 430074,China2 Hubei Subsurface Multi-scale Imaging Key Laboratory,Institute of Geophysics and Geomatics, China University of Geosciences (Wuhan),Wuhan 430074,China3 Exploration Technology and Engineering College,Hebei GEO University,Shijiazhuang 050031,China 4 Faculty of Civil Engineering and Geosciences,Delft University of Technology,Delft 2628 CN,the Netherlands

During drilling,the mud column sustains a slightly higher pressure than the formation to maintain the stability of the well wall,which causes the mud filtrate to penetrate into formation pores and displace in-situ fluids.The invasion depth is affected by reservoir properties,especially the reservoir permeability.Therefore,it is possible to estimate the reservoir permeability if the invasion depth can be measured.

Array induction logging;Permeability evaluation;Mud invasion

基于IGA算法的電阻率神經(jīng)網(wǎng)絡(luò)反演成像研究

Research of resistivity imaging using neural network based on immune genetic algorithm

GAO Ming-Liang，YU Sheng-Bao，ZHENG Jian-Bo，XU Chang，LIU Wei-Yu，LUAN Hui*

College of Instrumentation and Electrical Engineering，Jilin University，Changchun 130026，China

周峰，孟慶鑫，胡祥云等.2016.利用陣列感應(yīng)測(cè)井進(jìn)行儲(chǔ)層滲透率評(píng)價(jià).地球物理學(xué)報(bào)，59(11):4360-4371,

10.6038/cjg20161135.

Zhou F,Meng Q X,Hu X Y,et al.2016.Evaluation of reservoir permeability using array induction logging.Chinese J.Geophys.(in Chinese),59(11):4360-4371,doi:10.6038/cjg20161135.

國(guó)家自然科學(xué)基金(41674138，41304082，41304078)，中國(guó)石油科技創(chuàng)新基金(2015D-5006-0304)，油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室開(kāi)放基金(PRP/open-1302)資助.

周峰，男，1979年生，副教授，博士，研究方向?yàn)榫须姶欧?E-mail：zhoufeng617@163.com

*通信作者 胡祥云，男，1966年生，教授，博士生導(dǎo)師，研究方向?yàn)殡姶趴碧椒椒?E-mail：xyhu@163.com

10.6038/cjg20161135

P631 收稿日期2015-05-26，2016-09-20收修定稿

?? 中圖分類號(hào) P631

2015-09-02，2016-09-27收修定稿

A numerical study was conducted to investigate the feasibility of evaluating reservoir permeability with array induction logging.A mud invasion model was built up by coupling mud cake growth with multiple-phase fluid flow,and an array induction logging model was established based on the Born geometric factor theory.Joint forward simulations of mud invasion and array induction logging indicated that the responses of array induction logging can reflect the effect of mud invasion on the formation resistivity.Inversion based on the damped least square method revealed that the invasion depth can be acquired from array induction logging data.

We investigated the association between reservoir permeability and invasion depth,and found that in a reservoir with a permeability of 1 to 100 mD (1 mD= 0.987×10-3μm2),the reservoir permeability governs the invasion depth,and thus the permeability can be evaluated according to invasion depth.A two-dimensional numerical simulation showed that the inversed invasion depth curve had a similar fluctuation to the permeability variation.For a layered formation,a series of interpretation charts can be produced to evaluate the permeability of every layer with tolerable errors.The numerical investigation proves the feasibility of estimating reservoir permeability with array induction logging.