彭 杰，倪小威，徐思慧，湯 鵬

(1.長(zhǎng)江大學(xué) 油氣資源與勘探技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430100； 2.長(zhǎng)江大學(xué) 地球物理與石油資源學(xué)院，湖北 武漢 430100； 3.塔里木油田分公司油氣田產(chǎn)能建設(shè)事業(yè)部，新疆 庫(kù)爾勒 841000； 4.塔里木油田分公司勘探事業(yè)部，新疆 庫(kù)爾勒 841000； 5.中國(guó)石油塔里木油田分公司安全環(huán)保與工程監(jiān)督中心，新疆 庫(kù)爾勒 841000)

地球物理反演是根據(jù)觀(guān)測(cè)數(shù)據(jù)、利用數(shù)學(xué)方法來(lái)估算地球物理參數(shù)的分布狀態(tài)的一種有效手段。在電法測(cè)井反演中，通過(guò)對(duì)測(cè)井資料進(jìn)行快速有效的一維反演，可以初步反映地層原始特征，同時(shí)也可以為測(cè)井質(zhì)量好壞提供重要的參照[1]。但是這種實(shí)時(shí)反演技術(shù)對(duì)算法要求較高，要有較高的收斂速度。由于在反演過(guò)程中每一次迭代都會(huì)增大數(shù)值的計(jì)算量，這樣也會(huì)造成反演效率較低[2]。對(duì)一維反演算法的優(yōu)化處理可以從2個(gè)方向考慮：①重新尋找新的算法；②對(duì)傳統(tǒng)的一維算法進(jìn)行改進(jìn)，使其具有符合實(shí)際生產(chǎn)需求的精確度和速度。對(duì)測(cè)井?dāng)?shù)據(jù)反演地層電阻率的非線(xiàn)性問(wèn)題，國(guó)內(nèi)外學(xué)者進(jìn)行了許多探索。

本文以陣列側(cè)向測(cè)井?dāng)?shù)據(jù)為研究依據(jù)，一般的陣列側(cè)向反演算法有牛頓法、最小二乘法等[3-4]，這類(lèi)算法具有較快的收斂速度，但是反演結(jié)果受初始值影響大，如何進(jìn)行合理初值設(shè)定又是一個(gè)問(wèn)題。隨著反演算法的不斷發(fā)展，部分學(xué)者將模擬退火、遺傳算法及差分進(jìn)化算法引入了電測(cè)井資料的反演處理中[5-7]。這類(lèi)算法具備較強(qiáng)的全局搜索能力，同時(shí)反演結(jié)果幾乎不受初值的影響，但此類(lèi)算法往往收斂速度較慢，且編程復(fù)雜，不利用現(xiàn)場(chǎng)的推廣使用。隨著仿生學(xué)科的發(fā)展，1987年Kennedy和Eberhart基于對(duì)鳥(niǎo)群覓食規(guī)律的研究提出了粒子群優(yōu)化算法，粒子群算法由于結(jié)構(gòu)簡(jiǎn)單、易于編程實(shí)現(xiàn)、收斂速度快、搜索范圍大等優(yōu)點(diǎn)，已在不同領(lǐng)域得到了廣泛應(yīng)用[8]，同時(shí)也得到了油田現(xiàn)場(chǎng)工作者的青睞。并且在滿(mǎn)足反演精度的要求下，能夠替代反演過(guò)程中某些正演計(jì)算方法，目前主要以圖版數(shù)字化方法來(lái)代替嚴(yán)格正演計(jì)算[9]。

為了解決基本粒子群算法無(wú)法平衡局部搜索(決定收斂速度)與全局搜索能力(決定尋優(yōu)成功率及精度)的問(wèn)題，提出了自適應(yīng)非線(xiàn)性慣性權(quán)重，使得算法在迭代前期具有較強(qiáng)的全局搜索能力，利于尋找全局最優(yōu)解，而在后期慣性權(quán)重較大，加強(qiáng)算法的收斂速度，同時(shí)非線(xiàn)性變化的權(quán)重也加強(qiáng)了種群粒子間的信息交流，更利于全局最優(yōu)解的尋找。同時(shí)為進(jìn)一步提高反演效率，利用預(yù)先計(jì)算的陣列側(cè)向測(cè)井偽幾何因子數(shù)據(jù)表格代替嚴(yán)格正演計(jì)算，在反演過(guò)程中調(diào)用線(xiàn)性插值算法數(shù)據(jù)，也加快了反演的計(jì)算速度。采用此反演方法進(jìn)行三參數(shù)反演，反演結(jié)果與嚴(yán)格動(dòng)態(tài)泥漿侵入正演結(jié)果符合度較好。本文方法可為現(xiàn)場(chǎng)陣列側(cè)向測(cè)井?dāng)?shù)據(jù)的實(shí)時(shí)反演處理提供一定技術(shù)支持。

1 陣列側(cè)向電極系結(jié)構(gòu)

儀器的測(cè)井電極系結(jié)構(gòu)如圖1所示：主電極A0、A1與A1′分別為對(duì)應(yīng)的屏蔽電極，A2與A2′、A3與A3′、A4與A4′、A5與A5′也都為相對(duì)應(yīng)的屏幕電極，共5對(duì)；M1與M1′、M2與M2′、M3與M3′、M4與M4′、M5與M5′、M6與M6′、M7與M7′都是相對(duì)應(yīng)的監(jiān)督電極，共有7對(duì)，它們都呈主電極中心對(duì)稱(chēng)分布。儀器測(cè)量了4種探測(cè)深度的電阻率曲線(xiàn)，它們由淺到深的依次為MLR1、MLR2、MLR3、MLR4。

圖1 陣列側(cè)向測(cè)井電極系結(jié)構(gòu)Fig.1 Structure of array laterolog electrode system

2 有限元正演算法理論及偽幾何因子

求取側(cè)向測(cè)井正演響應(yīng)的實(shí)質(zhì)在于求取一光滑的電位函數(shù)μ對(duì)側(cè)向儀器形成的穩(wěn)流場(chǎng)進(jìn)行描述[10]，電位函數(shù)μ應(yīng)滿(mǎn)足式(1)：

(1)

式中，R為模型中不同區(qū)域的電阻率。在井眼中，R為鉆井液電阻率；在地層中，R為原狀地層電阻率。

給式(1)添加相對(duì)應(yīng)的邊界條件，規(guī)定在儀器的絕緣區(qū)域內(nèi)電位函數(shù)μ的法向?qū)?shù)為0；距儀器源無(wú)限遠(yuǎn)的區(qū)域，電位μ也為0。這樣就變成了求正演響應(yīng)函數(shù)的定解問(wèn)題。

對(duì)定解問(wèn)題的求解往往采取等效的思想將定解問(wèn)題轉(zhuǎn)化為泛函φ的極值問(wèn)題[11]：

(2)

式中，IE為電極系的電流；μE為電極系的電位；E為電極個(gè)數(shù)；積分區(qū)間包含儀器表面和無(wú)窮遠(yuǎn)邊界包圍的空間，對(duì)所有電極進(jìn)行求和。

利用有限元素法求取陣列側(cè)向測(cè)井正演響應(yīng)的詳細(xì)過(guò)程見(jiàn)文獻(xiàn)[12]。

通常引入偽幾何因子理論來(lái)描述[13]側(cè)向類(lèi)測(cè)井儀器的徑向探測(cè)特性。即視電阻率可以看作是沖洗帶、井眼及原狀地層電阻率的加權(quán)和，而不考慮地層縱向上的電阻率變化，即一維模型。沖洗帶電阻率的權(quán)值就是偽幾何因子λ，λ表征了陣列側(cè)向測(cè)井正演響應(yīng)受泥漿侵入的程度大小。經(jīng)過(guò)井眼校正后，排除該因素影響，那么視電阻率就等于沖洗帶及原始地層的綜合貢獻(xiàn)值，偽幾何因子λ可用式(3)來(lái)描述。

(3)

式中，RMLRi是隨著i變化的視電阻值；Rt為地層真實(shí)電阻率；Rxo為沖洗帶電阻率。

研究λ在不同Rt/Rxo的條件下隨泥漿侵入深度變化的特征，并制作MLR1—MLR4對(duì)應(yīng)的偽幾何因子數(shù)據(jù)表格。表1、表2分別為MLR1、MLR4的偽幾何因子數(shù)據(jù)，在實(shí)際反演過(guò)程中，利用線(xiàn)性表格插值代替嚴(yán)格正演計(jì)算。表中Di為泥漿侵入深度。

表1 MLR1偽幾何因子數(shù)據(jù)Tab.1 Pseudo-geometric factor data of MLR1

表2 MLR4偽幾何因子數(shù)據(jù)Tab.1 Pseudo-geometric factor data of MLR4

3 反演原理及算法

3.1 非線(xiàn)性反演原理

電阻率反演實(shí)質(zhì)就是通過(guò)實(shí)際測(cè)井資料來(lái)推斷模型電阻率大小。通過(guò)比較M個(gè)未知參數(shù)與N個(gè)實(shí)際數(shù)據(jù)，當(dāng)它們誤差在所設(shè)定的閾值以?xún)?nèi)，則可以把所求的M個(gè)參數(shù)看作實(shí)際地層參數(shù)[14]。研究基于陣列側(cè)向測(cè)井儀器的三參數(shù)反演，地層是水平均勻?qū)訝畹?，模型表示如下?/p>

(4)

式中，f(x)為所求目標(biāo)函數(shù)；xi為所求模型參數(shù)；yj為對(duì)應(yīng)j種方法的實(shí)際測(cè)井?dāng)?shù)據(jù)；Φj為對(duì)應(yīng)第j種方法的正演響應(yīng)算子，本文中，Φj滿(mǎn)足：

Φj=λj·Rxo+(1-λj)·Rt

(5)

式中，λj為正演數(shù)據(jù)線(xiàn)性插值結(jié)果。其中，j=1，2，3，4時(shí)分別對(duì)應(yīng)曲線(xiàn)為MLR1，MLR2，MLR3，MLR4。

3.2 基本粒子群算法

(6)

(7)

式中，c1、c2為加速因子，通常在[0 2]之間取值；r1、r2為[0 1]之間的隨機(jī)數(shù)。

將式(7)代入式(4)中，當(dāng)目標(biāo)函數(shù)小于預(yù)先設(shè)定的閾值或者迭代達(dá)到最大次數(shù)，輸出目標(biāo)函數(shù)最小的情況下對(duì)應(yīng)的位置屬性作為反演結(jié)果。

3.3 改進(jìn)粒子群算法

(1)采用混沌序列產(chǎn)生初始種群。混沌信號(hào)在全空間內(nèi)具有隨機(jī)性和遍歷性，這一特性能在求解最優(yōu)化問(wèn)題時(shí)，有效地避免算法陷入局部極值[16]。常用的產(chǎn)生混沌序列的方法是一維logistic映射，即：

rk+1=μrk(1-rk)

(8)

式中，rk為隨機(jī)數(shù)；k為迭代次數(shù)；μ為控制參數(shù)，當(dāng)μ=4時(shí)，能夠?qū)崿F(xiàn)[0 1]的完全映射，從而達(dá)到完全的混沌狀態(tài)，提高算法的全局尋優(yōu)能力。

將式(8)引入到初始種群產(chǎn)生中，即：

(9)

(2)自適應(yīng)非線(xiàn)性慣性權(quán)重?；玖Ｗ尤核惴o(wú)法平衡算法的前期尋優(yōu)及后期收斂能力，前人[17]為解決這個(gè)問(wèn)題，在速度位置的更新公式中引入了慣性權(quán)重ω，分別提出了固值慣性權(quán)重、線(xiàn)性慣性權(quán)重、非線(xiàn)性慣性權(quán)重等方法[18-19]。通常對(duì)于全局搜索比較好的優(yōu)化方法是在算法搜索前期就增強(qiáng)它的廣度探索能力以便得到優(yōu)秀的種子解，后期搜索算法應(yīng)具有較好的深度捜索能力以便加快算法的收斂速度[20]，即要求慣性權(quán)重前期大、后期小。本文基于一種自適應(yīng)非線(xiàn)性慣性權(quán)重，在滿(mǎn)足慣性權(quán)重前期大、后期小的特點(diǎn)的同時(shí)，強(qiáng)調(diào)慣性權(quán)重的非線(xiàn)性變化，保證粒子間信息的充分交流。慣性權(quán)重更新公式如下[21]：

(10)

(11)

(12)

式中，ωmax=0.9；ωmin=0.4。

(3)自適應(yīng)加速因子。加速因子c1決定粒子的自我學(xué)習(xí)能力，加速因子c2決定粒子的社會(huì)學(xué)習(xí)能力[22]。在反演初期，要求算法具備較強(qiáng)的社會(huì)學(xué)習(xí)能力，以擴(kuò)大種群多樣性，增加尋找全局最優(yōu)解的可能；在反演后期，要求算法具備較強(qiáng)的自我學(xué)習(xí)能力，增強(qiáng)局部搜索能力，加快算法收斂速度。為滿(mǎn)足該特點(diǎn)，將加速因子c1、c2作如下定義[23]：

(13)

(14)

式中，c1w=2.5；c1e=0.5；c2w=0.5；c2e=2.5。

本文改進(jìn)粒子群算法流程如圖2所示。

圖2 改進(jìn)算法流程Fig.2 Improved algorithm flow chart

4 算法評(píng)價(jià)

算法的評(píng)價(jià)指標(biāo)包括時(shí)間復(fù)雜度和空間復(fù)雜度，本文對(duì)遺傳類(lèi)算法性能的評(píng)價(jià)主要從尋優(yōu)成功率、平均收斂代數(shù)、平均最優(yōu)適應(yīng)度、最優(yōu)個(gè)體進(jìn)化曲線(xiàn)及抗噪性能等這些指標(biāo)進(jìn)行分析[24]。

表3記錄了改進(jìn)粒子群算法、基本粒子群算法、馬奎特算法在30次隨機(jī)計(jì)算中尋優(yōu)成功率、平均收斂代數(shù)、平均最優(yōu)適應(yīng)度等指標(biāo)的對(duì)比結(jié)果。從表3可知，改進(jìn)粒子群算法比后2種算法尋優(yōu)成功率更高。平均收斂代數(shù)越少代表收斂速度越快，改進(jìn)粒子群算法是基本粒子群算法收斂速度的2倍。本文是基于預(yù)測(cè)值與實(shí)際值的差值構(gòu)建的適應(yīng)度函數(shù)，最優(yōu)適應(yīng)度越小，表示預(yù)測(cè)值越接近真實(shí)值，算法性能越優(yōu)。改進(jìn)粒子群算法比基本粒子群算法性能更優(yōu)，稍遜于馬奎特算法。

表3 算法指標(biāo)對(duì)比Tab.3 Comparison of algorithm indexes

改進(jìn)粒子群算法及另外2種算法最優(yōu)個(gè)體進(jìn)化曲線(xiàn)如圖3所示。由圖3可知，改進(jìn)粒子群算法比基本粒子群算法迭代次數(shù)少很多，收斂速度顯然更快；馬奎特算法比粒子群類(lèi)算法有更大的初值適應(yīng)度，說(shuō)明初始值對(duì)馬奎特算法影響嚴(yán)重，對(duì)粒子群類(lèi)算法影響很小。

圖3 三種算法最優(yōu)個(gè)體進(jìn)化曲線(xiàn)對(duì)比Fig.3 Comparison of optimal individual evolution curves of three algorithms

對(duì)抗噪性能的評(píng)價(jià)，一般在測(cè)井響應(yīng)中加入高斯白噪聲，信噪比用式(15)表示[25]。

(15)

構(gòu)建徑向階躍三參數(shù)正演模型，模型參數(shù)設(shè)置見(jiàn)表4。

表4 地層模型Tab.4 Strata model

先對(duì)6個(gè)模型分別加入5%、10%、20%的高斯白噪聲，然后再進(jìn)行反演處理，所得結(jié)果見(jiàn)表5。

表5 不同信噪比下地層模型反演結(jié)果對(duì)比Tab.5 Comparison of formation model inversion results under different signal-to-noise ratios

加入0%、5%、10%、20%的高斯白噪聲進(jìn)行反演處理后，平均相對(duì)誤差分別為2.63%、3.10%、4.53%、8.71%。可以看出，改進(jìn)粒子群算法抗噪性能較好，當(dāng)加入20%的高斯白噪聲時(shí)，依然可以較好地反演出原狀地層信息，且反演后平均相對(duì)誤差小于10%。

5 反演效果驗(yàn)證

實(shí)際泥漿侵入地層，造成地層電阻率分布并不是階躍突變的，而是漸變的[26]?；跐B流—對(duì)流原理，可求得在不同時(shí)間內(nèi)泥漿侵入地層后地層的電阻率徑向分布特征，具體求解步驟可見(jiàn)文獻(xiàn)[27]。但在實(shí)際反演處理中，若考慮泥漿的動(dòng)態(tài)侵入，往往造成正演影響因素繁多，正演計(jì)算復(fù)雜且耗時(shí)長(zhǎng)，反演結(jié)果甚至往往不收斂，故提出了簡(jiǎn)易的泥漿侵入徑向階躍介質(zhì)模型，以此來(lái)簡(jiǎn)化反演過(guò)程中的正演計(jì)算。但這種反演模型的適用性及精度需與利用嚴(yán)格滲流—對(duì)流理論正演計(jì)算的結(jié)果進(jìn)行對(duì)比，分析其誤差，當(dāng)誤差在可控范圍內(nèi)時(shí)，此反演方法方可投入實(shí)際運(yùn)用。不同泥漿侵入條件下，陣列側(cè)向測(cè)井三參數(shù)反演結(jié)果與基于滲流—對(duì)流理論的正演結(jié)果對(duì)比如圖4、圖5所示。

圖4 泥漿高侵下反演效果對(duì)比Fig.4 Comparison of inversion results under high mud penetration

圖5 泥漿低侵下反演效果對(duì)比Fig.5 Comparison of inversion results under low mud penetration

實(shí)現(xiàn)基于滲流—對(duì)流理論的正演計(jì)算的基本參數(shù)設(shè)置如下：地層孔隙度10%，地層滲透率為5×10-3μm2，井眼與地層壓力差為0.1 MPa。在泥漿高侵條件下，地層水礦化度為18 g/L，泥漿礦化度為5 g/L；在泥漿低侵條件下，地層水礦化度為18 g/L，泥漿礦化度為36 g/L。

從圖4、圖5可知，陣列側(cè)向測(cè)井三參數(shù)反演的侵入帶電阻率基本代表了漸變侵入帶電阻率的加權(quán)平均和，可以在一定程度上反映侵入帶電阻率特征。最為重要的是，改進(jìn)粒子群算法可以準(zhǔn)確地反演出地層真實(shí)電阻率，且比基本粒子群算法評(píng)價(jià)效果更好；同時(shí)，改進(jìn)粒子群算法、基本粒子群算法都可以比較準(zhǔn)確地反演出侵入帶半徑的大小，說(shuō)明此方法可以在一定程度上反演出泥漿侵入后的地層徑向電阻率分布剖面，具備現(xiàn)場(chǎng)應(yīng)用價(jià)值。但此方法也有一定缺陷，在泥漿高侵過(guò)程中，實(shí)際地層電阻率分布會(huì)出現(xiàn)低阻環(huán)帶特征，但三參數(shù)反演模型并不能反映出此特征。針對(duì)這種情況，后期精細(xì)解釋時(shí)，可考慮進(jìn)行陣列側(cè)向測(cè)井四參數(shù)、五參數(shù)反演，提高反演模型適用性。

6 結(jié)論

(1)預(yù)先計(jì)算陣列側(cè)向測(cè)井偽幾何因子正演數(shù)據(jù)，在反演過(guò)程中直接進(jìn)行線(xiàn)性插值調(diào)用，相較于嚴(yán)格正演計(jì)算，在保證反演精度的前提下有效提高了反演效率。

(2)提出了一種改進(jìn)的粒子群算法，經(jīng)過(guò)與基本粒子群算法、馬奎特算法的綜合對(duì)比分析，發(fā)現(xiàn)該改進(jìn)算法在尋優(yōu)成功率、收斂速度、抗噪性等方面具備較明顯的優(yōu)勢(shì)。

(3)將反演結(jié)果與基于滲流理論的正演計(jì)算結(jié)果進(jìn)行對(duì)比驗(yàn)證，證明了反演模型的適用性及反演算法的準(zhǔn)確性，同時(shí)根據(jù)反演模型缺陷給出了合理建議。