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

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

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

吉林大學儀器科學與電氣工程學院，長春 130026

?

利用陣列感應(yīng)測井進行儲層滲透率評價

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

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

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

陣列感應(yīng)測井；滲透率評價；泥漿侵入

1 引言

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

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

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

2 泥漿侵入模型

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

2.1 多相流多組分模型

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

(1)

(2)

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

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

(3)

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

(4)

(5)

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

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

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

(6)

(7)

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

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

(8)

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

(9)

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

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

2.3 模型驗證

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

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

圖2 發(fā)表文獻中的侵入流速(a)和泥餅厚度(b)隨時間的變化曲線(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儲層厚度1m儲層外徑610m徑向網(wǎng)格數(shù)51儲層初始壓力25166kPa井底壓力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℃水相相對滲透率經(jīng)驗指數(shù)3油相相對滲透率經(jīng)驗指數(shù)3水相相對滲透率端點值30%油相相對滲透率端點值100%孔隙分布經(jīng)驗指數(shù)15毛細壓力系數(shù)2.74×10-2Pa·m徑向彌散度0m膠結(jié)指數(shù)2飽和度指數(shù)2巖性導電附加系數(shù)1

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

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

3.1 正演模型

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

(10)

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

(11)

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

(12)

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

3.2 反演方法

基于上述正演算法，應(yīng)用N個測井數(shù)據(jù)觀測值求取具有w個參數(shù)的正演模型來完成擬合反演，該過程的數(shù)學表達式寫為(李虎等，2013)：

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

(13)

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

(14)

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

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

(15)

其中J0為初始剛度陣，ε0為初始殘量.

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

(16)

(17)

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

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

3.3 方法驗證

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

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

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

圖4 侵入24 h間隔對應(yīng)的陣列感應(yīng)測井視電阻率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)測井響應(yīng)進行三參數(shù)反演，得到侵入帶電阻率、原狀地層電阻率和侵入深度的值如圖5所示，與泥漿侵入數(shù)值模擬進行對比可以看出反演結(jié)果較為準確地反映出受侵入影響儲層的電性結(jié)構(gòu)(圖3c)，對于本文所關(guān)心的侵入深度值而言，反演結(jié)果與數(shù)值模擬結(jié)果一致(圖3a)，證實了該反演方法的有效性.

4 滲透率評價應(yīng)用

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

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

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

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

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

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

4.2 半定量預(yù)測

以華北某油田測井和取心數(shù)據(jù)為例進行泥漿侵入的數(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，儲層和泥漿壓差為2 MPa，重力加速度為9.8 m·s-2，滲透率和孔隙度變化分別如圖9a、9b所示，這里的滲透率為水平滲透率，垂向滲透率為水平滲透率的1/10.

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

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

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

圖7 侵入24 h含水飽和度徑向分布(a)和儲層滲透率與侵入深度關(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時不同孔滲性儲層對應(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 反演的侵入深度(紅色)和儲層滲透率(藍色)對比Fig.11 Contrast between the curves of invasion depth (red) and permeability (blue)

圖12 二維軸對稱的分層地層模型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 圖版定量評估

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

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

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

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

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

5 結(jié)論

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

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

盡管目前還沒有一個理想的數(shù)學模型將泥漿侵入深度和儲層滲透率之間的相關(guān)性進行直接和準確的量化，但是本文的研究為利用陣列感應(yīng)測井進行滲透率定量評估提供可行性論證，也為將來借助陣列感應(yīng)測井數(shù)據(jù)對儲層滲透率進行直接準確的反演提供理論依據(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.

附中文參考文獻

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

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

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

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

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

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

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

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

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

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

(本文編輯 何燕)

高明亮，于生寶，鄭建波等.2016.基于IGA算法的電阻率神經(jīng)網(wǎng)絡(luò)反演成像研究.地球物理學報，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.

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

吉林大學儀器科學與電氣工程學院，長春 130026

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

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

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

作者簡介 高明亮，男，1987年生，博士，研究方向為直流電法非線性反演方法研究.E-mail：mlgao13@mails.jlu.edu.cn

*通訊作者 欒卉，女，1979年生，副教授，碩士生導師，研究方向為電磁勘探方法研究.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 引言

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

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

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

基于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ò)拓撲結(jié)構(gòu)如圖2a所示.神經(jīng)網(wǎng)絡(luò)包括輸入層、隱含層和輸出層.其中，輸入層神經(jīng)元的數(shù)量等于輸入變量的個數(shù)，輸出層神經(jīng)元的數(shù)量等于輸出變量的個數(shù)，隱含層神經(jīng)元的數(shù)量可以依據(jù)輸入信息量的大小和復雜度作相應(yīng)的改變.本文將X、Z位置和視電阻率作為神經(jīng)網(wǎng)絡(luò)的輸入，以真實電阻率作為神經(jīng)網(wǎng)絡(luò)的輸出，這樣的數(shù)據(jù)類型最適合神經(jīng)網(wǎng)絡(luò)的訓練和測試.具體討論參見第3節(jié).BP神經(jīng)網(wǎng)絡(luò)的主要思想就是通過樣本集訓練網(wǎng)絡(luò)，使得神經(jīng)網(wǎng)絡(luò)的輸出與期望輸出的誤差平方和達到最小.它是通過反復迭代使相對誤差函數(shù)在斜率下降的方向上計算網(wǎng)絡(luò)權(quán)閾值和偏差的變化而逐漸達到目標.每一次權(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ò)拓撲結(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ù)表達式如式(1)所示：

(1)

輸入層到隱含層單元J表達式如下：

(2)

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

隱含層單元J到輸出層aj表達式為

(3)

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

(4)

式中，ε為學習動量，η為學習率，t為訓練次數(shù).

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

(5)

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

2.2 免疫遺傳算法

免疫遺傳算法就是將免疫理論(Immune Algorithm，IA)和基本遺傳算法(Simple Genetic Algorithm，SGA)各自的優(yōu)點結(jié)合起來的一個多學科相互交叉、滲透的優(yōu)化算法.免疫遺傳算法可以看作一種新型融合算法，羅小平(2002)在遺傳算法中加入免疫算子，使遺傳算法變成具有免疫功能的新算法，是一種改進的遺傳算法，該算法既保留了遺傳算法的搜索特性，又利用了免疫算法的多機制求解多目標函數(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ù).對神經(jīng)網(wǎng)絡(luò)反演參數(shù)進行優(yōu)化設(shè)計時，應(yīng)在保證神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)穩(wěn)定的前提下，力求加快網(wǎng)絡(luò)收斂速度，以實現(xiàn)高精度反演.在優(yōu)化設(shè)計時，本文使用的目標函數(shù)與約束函數(shù)均是高度非線性的，屬于多目標優(yōu)化問題，適合于免疫遺傳算法的求解.

2.3.1 建立目標函數(shù)

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

(6)

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

(7)

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

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

(1) 編碼

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

(2) 初始抗體的生成

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

(3) 疫苗的提取

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

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

用設(shè)定的交叉概率Pc和變異概率Pm，按照標準遺傳算法進行交叉和變異操作.

(5) 接種疫苗

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

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

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

(8)

(7) 免疫機制選擇

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

(8) 終止條件的判斷

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

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

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

(c) 達到最大優(yōu)化時間.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

表3 三種反演算法性能對比

表4 三種反演算法試驗結(jié)果對比

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

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

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

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

表5 反演算法試驗結(jié)果對比

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

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

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

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

模型2為3個異常體的組合模型，主要用于檢驗IGA-BP和GABP、BP的反演性能的差異，模型2的基本參數(shù)設(shè)置如下：依然采用溫納-斯倫貝格裝置，接收電極37個，電極極距1 m，背景場的電阻率為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ò)反演測試中，神經(jīng)網(wǎng)絡(luò)輸入節(jié)點為 X、Z位置和視電阻率數(shù)據(jù)，網(wǎng)絡(luò)輸出為真實電阻率參數(shù)，網(wǎng)絡(luò)輸出的反演結(jié)果如圖7所示.

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

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

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

6 結(jié)論

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

致謝 感謝吉林大學地球信息探測儀器教育部重點實驗室提供了研究平臺，感謝吉林大學時間域電磁組的全體成員給予的幫助和指導，感謝審稿專家對該文的審定.

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．

附中文參考文獻

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

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

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

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

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

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

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

周明，孫樹棟.1999.遺傳算法原理與應(yīng)用.北京:國防工業(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)測井進行儲層滲透率評價.地球物理學報，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.

國家自然科學基金(41674138，41304082，41304078)，中國石油科技創(chuàng)新基金(2015D-5006-0304)，油氣資源與探測國家重點實驗室開放基金(PRP/open-1302)資助.

周峰，男，1979年生，副教授，博士，研究方向為井中電磁法.E-mail：zhoufeng617@163.com

*通信作者 胡祥云，男，1966年生，教授，博士生導師，研究方向為電磁勘探方法.E-mail：xyhu@163.com

10.6038/cjg20161135

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

?? 中圖分類號 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.