• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

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

    2016-11-24 02:05:21周峰孟慶鑫胡祥云EvertSlob潘和平馬火林
    地球物理學(xué)報(bào) 2016年11期
    關(guān)鍵詞:BP神經(jīng)網(wǎng)絡(luò)

    周峰,孟慶鑫,胡祥云,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.

    猜你喜歡
    BP神經(jīng)網(wǎng)絡(luò)
    基于神經(jīng)網(wǎng)絡(luò)的北京市房?jī)r(jià)預(yù)測(cè)研究
    商情(2016年43期)2016-12-23 14:23:13
    一種基于OpenCV的車牌識(shí)別方法
    基于遺傳算法—BP神經(jīng)網(wǎng)絡(luò)的乳腺腫瘤輔助診斷模型
    一種基于改進(jìn)BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)T/R組件溫度的方法
    基于BP神經(jīng)網(wǎng)絡(luò)的光通信系統(tǒng)故障診斷
    科技視界(2016年26期)2016-12-17 17:57:49
    提高BP神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)速率的算法研究
    考試周刊(2016年21期)2016-12-16 11:02:03
    就bp神經(jīng)網(wǎng)絡(luò)銀行選址模型的相關(guān)研究
    基于DEA—GA—BP的建設(shè)工程評(píng)標(biāo)方法研究
    基于BP神經(jīng)網(wǎng)絡(luò)的旅行社發(fā)展方向研究
    商情(2016年39期)2016-11-21 09:30:36
    復(fù)雜背景下的手勢(shì)識(shí)別方法
    国产精品人妻久久久久久| 亚洲久久久久久中文字幕| 国产精品精品国产色婷婷| 我的老师免费观看完整版| 亚洲va日本ⅴa欧美va伊人久久| 欧美国产日韩亚洲一区| 性色avwww在线观看| 五月玫瑰六月丁香| 精品人妻1区二区| 国产免费男女视频| 成人av在线播放网站| 99精品久久久久人妻精品| 岛国在线免费视频观看| 亚洲av免费在线观看| 天堂√8在线中文| 欧美黄色片欧美黄色片| 国产黄a三级三级三级人| 男女床上黄色一级片免费看| 91av网一区二区| 97人妻精品一区二区三区麻豆| 大型黄色视频在线免费观看| 成年版毛片免费区| 亚洲午夜理论影院| 午夜视频国产福利| 亚洲中文日韩欧美视频| 亚洲欧美精品综合久久99| 成人鲁丝片一二三区免费| 高清毛片免费观看视频网站| a在线观看视频网站| 麻豆国产av国片精品| 色播亚洲综合网| 毛片一级片免费看久久久久 | 精品久久久久久久末码| 69人妻影院| 色播亚洲综合网| 精品欧美国产一区二区三| 美女 人体艺术 gogo| 国产成人aa在线观看| 别揉我奶头~嗯~啊~动态视频| 午夜精品久久久久久毛片777| 美女cb高潮喷水在线观看| 禁无遮挡网站| 国产淫片久久久久久久久 | 日韩欧美国产在线观看| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 好看av亚洲va欧美ⅴa在| 免费av观看视频| 免费无遮挡裸体视频| 午夜福利欧美成人| 日韩欧美三级三区| 国产亚洲av嫩草精品影院| 一本精品99久久精品77| 悠悠久久av| 久久精品影院6| 精品一区二区三区视频在线| 少妇的逼水好多| 国产三级中文精品| 搡老妇女老女人老熟妇| 精品久久久久久久久久久久久| 亚洲av五月六月丁香网| 国产欧美日韩一区二区三| 老鸭窝网址在线观看| 亚洲男人的天堂狠狠| 欧美不卡视频在线免费观看| avwww免费| 女人十人毛片免费观看3o分钟| 色综合欧美亚洲国产小说| 免费看光身美女| 婷婷色综合大香蕉| 嫩草影院入口| 精品久久久久久久久久免费视频| 一区二区三区高清视频在线| 男女床上黄色一级片免费看| 亚洲自偷自拍三级| 国产精品人妻久久久久久| 可以在线观看毛片的网站| 亚洲国产日韩欧美精品在线观看| 亚洲成av人片在线播放无| 日韩国内少妇激情av| 日本成人三级电影网站| 伊人久久精品亚洲午夜| 国产视频一区二区在线看| 亚洲成人久久性| 好男人电影高清在线观看| 深夜a级毛片| 99国产精品一区二区蜜桃av| 国产中年淑女户外野战色| 看免费av毛片| 国产白丝娇喘喷水9色精品| 成人国产一区最新在线观看| 久久国产乱子伦精品免费另类| 91九色精品人成在线观看| 免费av不卡在线播放| 亚洲真实伦在线观看| 久久99热这里只有精品18| 欧美性感艳星| 别揉我奶头~嗯~啊~动态视频| 久久伊人香网站| 亚洲国产日韩欧美精品在线观看| 欧美不卡视频在线免费观看| 九九在线视频观看精品| 三级男女做爰猛烈吃奶摸视频| 老熟妇仑乱视频hdxx| 五月伊人婷婷丁香| 99在线人妻在线中文字幕| 亚洲av一区综合| 亚洲欧美清纯卡通| 女人被狂操c到高潮| 中文字幕人妻熟人妻熟丝袜美| 日韩精品中文字幕看吧| 最新中文字幕久久久久| 久久午夜亚洲精品久久| netflix在线观看网站| 亚洲成人精品中文字幕电影| 午夜老司机福利剧场| 99热这里只有是精品50| 一区二区三区高清视频在线| 色综合欧美亚洲国产小说| 国产免费男女视频| 美女黄网站色视频| 成年免费大片在线观看| 欧美最新免费一区二区三区 | 午夜福利视频1000在线观看| 波多野结衣巨乳人妻| 午夜福利18| 欧美成人性av电影在线观看| 中文字幕人妻熟人妻熟丝袜美| 国产精品一区二区性色av| 九九热线精品视视频播放| 女人被狂操c到高潮| 白带黄色成豆腐渣| 国产精品一区二区三区四区久久| 亚洲欧美激情综合另类| 亚洲色图av天堂| 国产精品,欧美在线| 1024手机看黄色片| 久久精品国产亚洲av涩爱 | 丁香六月欧美| www日本黄色视频网| 国产精品久久视频播放| 色综合婷婷激情| 久久国产乱子伦精品免费另类| 毛片一级片免费看久久久久 | 成人毛片a级毛片在线播放| 成人永久免费在线观看视频| 精品99又大又爽又粗少妇毛片 | 精品久久久久久久人妻蜜臀av| 欧美不卡视频在线免费观看| 亚洲中文字幕一区二区三区有码在线看| 人人妻人人澡欧美一区二区| 精品乱码久久久久久99久播| 一区二区三区四区激情视频 | 亚洲av中文字字幕乱码综合| 欧美日韩黄片免| 97超视频在线观看视频| 欧美黄色淫秽网站| 人妻夜夜爽99麻豆av| 国产精品女同一区二区软件 | 在线天堂最新版资源| 99在线人妻在线中文字幕| 欧美丝袜亚洲另类 | 久久精品夜夜夜夜夜久久蜜豆| 久久久久九九精品影院| 中文资源天堂在线| 丰满人妻一区二区三区视频av| 三级毛片av免费| 欧美日韩瑟瑟在线播放| 国产av不卡久久| 欧美3d第一页| 我的老师免费观看完整版| 免费大片18禁| 淫妇啪啪啪对白视频| 午夜福利在线在线| 久久九九热精品免费| 欧美高清成人免费视频www| 久久精品夜夜夜夜夜久久蜜豆| 午夜精品在线福利| 亚洲av五月六月丁香网| 国产黄a三级三级三级人| 亚洲精品456在线播放app | 亚洲熟妇中文字幕五十中出| 婷婷丁香在线五月| 国产乱人视频| 18美女黄网站色大片免费观看| 99久久精品一区二区三区| 国产男靠女视频免费网站| 国产伦精品一区二区三区视频9| 欧美乱色亚洲激情| 色5月婷婷丁香| 一区二区三区高清视频在线| 老熟妇仑乱视频hdxx| 97超视频在线观看视频| 久久6这里有精品| 国产精品亚洲av一区麻豆| 免费在线观看日本一区| 久久精品91蜜桃| 麻豆国产97在线/欧美| 亚洲精品久久国产高清桃花| 亚洲男人的天堂狠狠| 黄色配什么色好看| ponron亚洲| 久久久久久久久久黄片| 狠狠狠狠99中文字幕| 精品无人区乱码1区二区| 在线国产一区二区在线| av天堂在线播放| 超碰av人人做人人爽久久| 国产成人欧美在线观看| 人妻制服诱惑在线中文字幕| 午夜福利欧美成人| 性插视频无遮挡在线免费观看| 国产精品久久久久久久久免 | 国产高清有码在线观看视频| 午夜精品在线福利| 欧美3d第一页| 老司机深夜福利视频在线观看| 精品久久久久久,| 一级毛片久久久久久久久女| 欧美中文日本在线观看视频| 国产探花极品一区二区| 婷婷色综合大香蕉| 丝袜美腿在线中文| 亚洲国产精品久久男人天堂| 精品福利观看| 国产免费av片在线观看野外av| 欧美成人a在线观看| 黄色丝袜av网址大全| 亚洲成人免费电影在线观看| 国产私拍福利视频在线观看| 国产精华一区二区三区| 成人永久免费在线观看视频| 国产av不卡久久| 内地一区二区视频在线| 美女大奶头视频| av女优亚洲男人天堂| 两性午夜刺激爽爽歪歪视频在线观看| www.999成人在线观看| 亚洲欧美精品综合久久99| 一进一出抽搐gif免费好疼| 一区二区三区四区激情视频 | 69人妻影院| 美女被艹到高潮喷水动态| 亚洲美女视频黄频| 特大巨黑吊av在线直播| 婷婷精品国产亚洲av在线| 俺也久久电影网| 99国产精品一区二区三区| 长腿黑丝高跟| 国产又黄又爽又无遮挡在线| 午夜福利欧美成人| 成人鲁丝片一二三区免费| 麻豆久久精品国产亚洲av| 欧美在线一区亚洲| 久久精品国产亚洲av涩爱 | 白带黄色成豆腐渣| 日本免费a在线| 精品国内亚洲2022精品成人| 国产91精品成人一区二区三区| 国产亚洲欧美98| 1000部很黄的大片| 青草久久国产| 欧美最新免费一区二区三区 | 中文字幕av在线有码专区| 久久精品久久久久久噜噜老黄 | 国产伦一二天堂av在线观看| 亚洲国产欧美人成| 夜夜躁狠狠躁天天躁| 男女下面进入的视频免费午夜| 最近最新中文字幕大全电影3| 好男人在线观看高清免费视频| 桃色一区二区三区在线观看| 中文字幕av在线有码专区| av在线蜜桃| 亚洲国产欧美人成| 欧美色欧美亚洲另类二区| 日本熟妇午夜| 国产欧美日韩一区二区三| 久久久色成人| 怎么达到女性高潮| 国产三级中文精品| 亚洲av一区综合| 久久精品国产亚洲av天美| 日韩中文字幕欧美一区二区| 变态另类丝袜制服| 色视频www国产| 波野结衣二区三区在线| 国产成人福利小说| 色综合欧美亚洲国产小说| 波野结衣二区三区在线| 美女高潮喷水抽搐中文字幕| 亚洲狠狠婷婷综合久久图片| 天堂动漫精品| 非洲黑人性xxxx精品又粗又长| 99热这里只有是精品50| 亚洲专区中文字幕在线| 亚洲av一区综合| 日本在线视频免费播放| 永久网站在线| 听说在线观看完整版免费高清| 成人亚洲精品av一区二区| 永久网站在线| 在线观看美女被高潮喷水网站 | 国产在线精品亚洲第一网站| 大型黄色视频在线免费观看| 观看美女的网站| 日韩欧美在线二视频| 亚洲第一区二区三区不卡| 18禁在线播放成人免费| 国产精品永久免费网站| 久久精品人妻少妇| 日日干狠狠操夜夜爽| 男女那种视频在线观看| 久久久久亚洲av毛片大全| 欧美成人a在线观看| 亚洲午夜理论影院| 色综合亚洲欧美另类图片| 亚洲人成伊人成综合网2020| 国产白丝娇喘喷水9色精品| 九色成人免费人妻av| 精品久久久久久成人av| 91字幕亚洲| 欧美zozozo另类| 中文字幕久久专区| 18禁黄网站禁片免费观看直播| 久久性视频一级片| 欧美高清成人免费视频www| 亚洲成人中文字幕在线播放| 在线观看午夜福利视频| 99精品在免费线老司机午夜| 久久热精品热| 日本a在线网址| 真实男女啪啪啪动态图| 久久久久亚洲av毛片大全| 亚洲精品影视一区二区三区av| 人妻丰满熟妇av一区二区三区| 欧美日韩综合久久久久久 | 国内少妇人妻偷人精品xxx网站| 免费观看人在逋| 国产 一区 欧美 日韩| 午夜福利在线在线| 色精品久久人妻99蜜桃| 美女被艹到高潮喷水动态| 丁香欧美五月| 美女cb高潮喷水在线观看| 欧美日本视频| 国产私拍福利视频在线观看| 亚洲欧美激情综合另类| 欧美最新免费一区二区三区 | 亚洲一区高清亚洲精品| 亚洲人成网站在线播放欧美日韩| 久久久久久久午夜电影| 久久草成人影院| 99久久久亚洲精品蜜臀av| 91在线精品国自产拍蜜月| 99久久成人亚洲精品观看| 日韩人妻高清精品专区| 午夜激情福利司机影院| 国产乱人视频| 99热这里只有是精品在线观看 | 亚洲精品乱码久久久v下载方式| 成人三级黄色视频| 精品午夜福利视频在线观看一区| 少妇被粗大猛烈的视频| 99热这里只有精品一区| 国产高清激情床上av| 中文字幕免费在线视频6| 少妇丰满av| 成人亚洲精品av一区二区| 国产一区二区在线观看日韩| 日本三级黄在线观看| 国产高清激情床上av| 国产久久久一区二区三区| 婷婷六月久久综合丁香| 观看免费一级毛片| 天堂√8在线中文| 99精品在免费线老司机午夜| 亚洲无线在线观看| 久久久久免费精品人妻一区二区| 亚洲av第一区精品v没综合| 亚洲va日本ⅴa欧美va伊人久久| 天美传媒精品一区二区| 久久久久久久亚洲中文字幕 | 精品午夜福利视频在线观看一区| 婷婷色综合大香蕉| 人妻久久中文字幕网| 五月伊人婷婷丁香| 国内久久婷婷六月综合欲色啪| 国产蜜桃级精品一区二区三区| 高清在线国产一区| 精品乱码久久久久久99久播| 国产精品电影一区二区三区| 噜噜噜噜噜久久久久久91| 日韩欧美在线乱码| 亚洲av第一区精品v没综合| 久久久久国产精品人妻aⅴ院| 91久久精品国产一区二区成人| 国产乱人伦免费视频| 亚洲精品影视一区二区三区av| 欧美日韩综合久久久久久 | 在线看三级毛片| 午夜影院日韩av| 亚洲精品456在线播放app | 天天躁日日操中文字幕| 99视频精品全部免费 在线| 毛片女人毛片| 久久久久久九九精品二区国产| 少妇人妻精品综合一区二区 | 18禁黄网站禁片午夜丰满| 神马国产精品三级电影在线观看| 日本黄色片子视频| 午夜a级毛片| 精品久久久久久久久久久久久| 丰满人妻一区二区三区视频av| 天堂网av新在线| 国产一区二区亚洲精品在线观看| 好看av亚洲va欧美ⅴa在| 黄色视频,在线免费观看| 国产色爽女视频免费观看| 久久亚洲精品不卡| 夜夜夜夜夜久久久久| 波多野结衣高清无吗| 精品欧美国产一区二区三| 变态另类成人亚洲欧美熟女| 又黄又爽又刺激的免费视频.| 亚洲美女黄片视频| 在线观看午夜福利视频| 欧美不卡视频在线免费观看| 性色av乱码一区二区三区2| 中文字幕av在线有码专区| 国产一区二区激情短视频| 国产aⅴ精品一区二区三区波| 欧美在线黄色| 99久久精品国产亚洲精品| 身体一侧抽搐| 欧美成人a在线观看| 99热这里只有是精品50| 久久亚洲精品不卡| 成年免费大片在线观看| 欧美一区二区精品小视频在线| 中文字幕熟女人妻在线| 啦啦啦观看免费观看视频高清| 色综合站精品国产| 欧美激情国产日韩精品一区| 老鸭窝网址在线观看| 日韩中文字幕欧美一区二区| 51国产日韩欧美| 色在线成人网| 99视频精品全部免费 在线| 亚洲欧美日韩高清专用| 免费看a级黄色片| 精品日产1卡2卡| 亚洲激情在线av| 亚洲欧美日韩高清在线视频| 一个人免费在线观看的高清视频| 国产毛片a区久久久久| 99riav亚洲国产免费| 一进一出抽搐gif免费好疼| 免费在线观看成人毛片| 日本 av在线| 精品久久久久久久久久免费视频| 久久久精品欧美日韩精品| 欧美激情国产日韩精品一区| 日韩欧美精品免费久久 | 国产精品嫩草影院av在线观看 | 1000部很黄的大片| 午夜免费激情av| 亚洲美女搞黄在线观看 | 国产精品自产拍在线观看55亚洲| 国产精品一区二区免费欧美| 男女做爰动态图高潮gif福利片| 国产精品久久久久久久久免 | 日韩欧美免费精品| 窝窝影院91人妻| 亚洲自偷自拍三级| 中文字幕免费在线视频6| 亚洲成人精品中文字幕电影| 麻豆国产av国片精品| 成人高潮视频无遮挡免费网站| 国产精品日韩av在线免费观看| 国产探花在线观看一区二区| 日韩欧美三级三区| 别揉我奶头 嗯啊视频| 一a级毛片在线观看| 又爽又黄无遮挡网站| 日本与韩国留学比较| 亚洲精品456在线播放app | 亚洲熟妇熟女久久| 真实男女啪啪啪动态图| 91麻豆精品激情在线观看国产| 蜜桃久久精品国产亚洲av| 午夜a级毛片| 日韩欧美在线二视频| 男女下面进入的视频免费午夜| 老司机福利观看| av欧美777| 日韩欧美在线二视频| 一区福利在线观看| 欧美高清性xxxxhd video| 久久亚洲精品不卡| 亚洲不卡免费看| 国产精品av视频在线免费观看| 精品午夜福利视频在线观看一区| av国产免费在线观看| 99久国产av精品| 亚洲成人中文字幕在线播放| 人妻久久中文字幕网| 国产不卡一卡二| 精品99又大又爽又粗少妇毛片 | 国产av不卡久久| 午夜精品一区二区三区免费看| 国产欧美日韩精品亚洲av| 中文字幕久久专区| 九色成人免费人妻av| 99国产精品一区二区三区| 日本一二三区视频观看| 久久人人精品亚洲av| 99久久成人亚洲精品观看| 欧美性猛交黑人性爽| 99久久九九国产精品国产免费| 久9热在线精品视频| 亚洲成人中文字幕在线播放| 人妻丰满熟妇av一区二区三区| 性色av乱码一区二区三区2| 久久精品91蜜桃| 国产成人av教育| 久久亚洲精品不卡| 少妇的逼水好多| 亚洲精华国产精华精| 久久精品人妻少妇| 国内久久婷婷六月综合欲色啪| 琪琪午夜伦伦电影理论片6080| 亚洲av五月六月丁香网| 一进一出抽搐动态| 黄色女人牲交| 国产精品不卡视频一区二区 | 97碰自拍视频| 精品久久久久久久久久久久久| 一个人免费在线观看的高清视频| 小说图片视频综合网站| 国产亚洲欧美98| 欧美黄色片欧美黄色片| 欧美日本视频| 国产亚洲精品久久久com| 日韩成人在线观看一区二区三区| 真实男女啪啪啪动态图| 一夜夜www| 亚洲第一区二区三区不卡| 色视频www国产| 亚洲av免费在线观看| 久久精品国产亚洲av涩爱 | 97人妻精品一区二区三区麻豆| 久久人妻av系列| av在线天堂中文字幕| 亚洲成人免费电影在线观看| 国产单亲对白刺激| 精品乱码久久久久久99久播| 午夜日韩欧美国产| 一本久久中文字幕| 91av网一区二区| 夜夜躁狠狠躁天天躁| 日本免费a在线| 悠悠久久av| 内地一区二区视频在线| 亚洲av电影在线进入| 国产精品亚洲一级av第二区| 夜夜夜夜夜久久久久| 亚洲18禁久久av| 黄色一级大片看看| 少妇人妻精品综合一区二区 | 亚洲男人的天堂狠狠| 欧美性猛交黑人性爽| 国产亚洲精品久久久久久毛片| 亚洲成a人片在线一区二区| 免费电影在线观看免费观看| 国产精品亚洲美女久久久| 好男人在线观看高清免费视频| 久久精品久久久久久噜噜老黄 | 午夜影院日韩av| 欧美一级a爱片免费观看看| 两个人视频免费观看高清| 91久久精品电影网| www.熟女人妻精品国产| 俄罗斯特黄特色一大片| 国产精品一区二区性色av| 国产精品久久电影中文字幕| 精品久久久久久久久av| 久久国产精品影院| 女人被狂操c到高潮| 99在线人妻在线中文字幕| 黄色日韩在线| 首页视频小说图片口味搜索| 免费人成视频x8x8入口观看| 成人国产一区最新在线观看| 国产aⅴ精品一区二区三区波| 日本a在线网址| 又粗又爽又猛毛片免费看| 最近最新中文字幕大全电影3| 日本a在线网址| 成人国产一区最新在线观看| 国产 一区 欧美 日韩| 两个人视频免费观看高清| 嫩草影视91久久| 国产高清三级在线| 两个人视频免费观看高清| 日韩欧美国产一区二区入口| 亚洲在线自拍视频| 99久久无色码亚洲精品果冻| 2021天堂中文幕一二区在线观| 免费观看精品视频网站| 一区二区三区免费毛片| 久久国产乱子免费精品| 赤兔流量卡办理| 欧美乱妇无乱码| 久久精品国产亚洲av香蕉五月| 国产成+人综合+亚洲专区|