基于混合模擬退火算法的陣列側(cè)向測(cè)井實(shí)時(shí)反演研究

馮 進(jìn)， 倪小威， 楊 清， 管 耀， 劉迪仁

（1. 中海石油（中國(guó)）有限公司深圳分公司研究院，廣東深圳 518054；2. 中國(guó)石油塔里木油田分公司油氣田產(chǎn)能建設(shè)事業(yè)部，新疆庫(kù)爾勒 841000；3. 油氣資源與勘探技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室（長(zhǎng)江大學(xué)），湖北武漢 430100；4. 長(zhǎng)江大學(xué)地球物理與石油資源學(xué)院，湖北武漢 430100）

電阻率測(cè)井資料往往受鉆井液侵入和井眼、圍巖等環(huán)境因素的影響，不能準(zhǔn)確反映鉆井液侵入剖面信息，故需要對(duì)其進(jìn)行反演或校正處理[1-6]。早期為解決電阻率失真問(wèn)題，油田現(xiàn)場(chǎng)多采用圖版校正法[7]。這種方法操作比較簡(jiǎn)便，但需要多次進(jìn)行插值處理，且無(wú)法反映全部地層情況，準(zhǔn)確性和普適性都較低[8]。

隨著計(jì)算機(jī)技術(shù)的發(fā)展，線性反演技術(shù)被引入電測(cè)井資料處理[9-10]。其中，馬奎特算法是一種線性反演算法，其收斂速度極快，有利于測(cè)井資料的實(shí)時(shí)反演，但在實(shí)際應(yīng)用中其反演結(jié)果受初始值影響極大，若初始值設(shè)置不合理，很難準(zhǔn)確反演出地層的真實(shí)電阻率值[11]。為了提高反演精度，針對(duì)同時(shí)存在多種電測(cè)井資料的油井，提出了聯(lián)合反演的技術(shù)思路，通過(guò)增加反演過(guò)程中的地層信息來(lái)提高反演精度，此類(lèi)方法具備一定的理論可行性[12]，但實(shí)際上極少有油井進(jìn)行2 種以上不同原理的電測(cè)井作業(yè)，因此其普適性并不高。陣列型電測(cè)井儀器（如陣列側(cè)向測(cè)井儀器、陣列感應(yīng)測(cè)井儀器）的探測(cè)深度深、分辨率高[13-14]，能夠提供多條測(cè)井信息，可以很好地解決聯(lián)合反演時(shí)存在的問(wèn)題。隨著反演技術(shù)的進(jìn)一步發(fā)展，非線性反演算法被引入測(cè)井反演領(lǐng)域[15-20]。模擬退火算法是一種典型的非線性反演算法，具備較強(qiáng)的全局搜索能力，在計(jì)算時(shí)間充足的情況下基本上能找到全局最優(yōu)解；但模擬退火算法在處理復(fù)雜的非線性?xún)?yōu)化問(wèn)題時(shí)收斂速度過(guò)慢，難以滿(mǎn)足測(cè)井資料的實(shí)時(shí)反演需求。

為解決傳統(tǒng)模擬退火算法與馬奎特算法存在的問(wèn)題，筆者將2 種反演方法結(jié)合，提出了混合模擬退火算法：利用幅度差法優(yōu)化了模擬退火算法的初始值選擇，然后利用退火策略進(jìn)行了迭代，最后將迭代過(guò)一定次數(shù)的結(jié)果作為馬奎特算法的初始值，利用馬奎特算法進(jìn)行迭代尋優(yōu)?；旌夏M退火算法既保留了模擬退火算法的全局尋優(yōu)能力，也很好地利用了馬奎特算法的后期收斂能力，不僅能滿(mǎn)足實(shí)時(shí)反演的要求，而且還可以進(jìn)一步提高反演精度。

1 陣列側(cè)向測(cè)井反演模型

實(shí)際測(cè)井過(guò)程中，陣列側(cè)向測(cè)井儀器可獲得由淺至深的4 條視電阻率曲線MLR1、MLR2、MLR3和MLR4。視電阻率響應(yīng)可以認(rèn)為是沖洗帶半徑rxo、沖洗帶電阻率 Rxo、地層電阻率 Rt、井眼和圍巖等地層參數(shù)的非線性函數(shù)[21-23]。實(shí)時(shí)反演過(guò)程中，井眼校正時(shí)忽略圍巖的影響，故視電阻率可以被認(rèn)為是三參數(shù)（沖洗帶半徑 rxo、沖洗帶電阻率 Rxo和地層電阻率 Rt）的非線性函數(shù)。陣列側(cè)向測(cè)井儀器在水平層狀介質(zhì)中的電阻率反演模型可表示為：

式中： f(X) 為待求目標(biāo)函數(shù)； xi表示待求模型參數(shù)；i 為反演種群規(guī)模； yj為第j 種測(cè)井方法的實(shí)際測(cè)井?dāng)?shù)據(jù)； Φj為第j 種測(cè)井方法的正演響應(yīng)算子，可通過(guò)模式匹配法或者有限元法得到。其中， j=1,2,3,4，分別對(duì)應(yīng)MLR1，MLR2，MLR3 和MLR4（下同）。

2 混合模擬退火算法

2.1 模擬退火算法原理及步驟

模擬退火算法是在金屬退火機(jī)制上演化而成的一種非線性反演算法，從概率意義上來(lái)說(shuō)，模擬退火算法總能找到全局最小點(diǎn)[24]。其控制因素主要分為搜索空間 Ω、能量函數(shù) f(x)、 狀態(tài)轉(zhuǎn)移規(guī)則 p和冷卻進(jìn)度表 T(x)等4 類(lèi)。

模擬退火算法的一般步驟是：首先隨機(jī)給定初始解，然后借助控制參數(shù)產(chǎn)生的一系列Mapkob 鏈、新解產(chǎn)生裝置和接受準(zhǔn)則，重復(fù)“產(chǎn)生新解—計(jì)算適應(yīng)度函數(shù)—判斷是否接受新解—接受/拋棄新解”的步驟，不斷進(jìn)行迭代，直到目標(biāo)函數(shù)適應(yīng)度達(dá)到最小。

2.1.1 初始值的選取

選取的初始值是否合適會(huì)直接影響模擬退火算法的收斂速度[25]。初始值若偏離真實(shí)值較小，模擬退火時(shí)只需進(jìn)行少數(shù)迭代即可尋找到最優(yōu)解；若初始值偏離真實(shí)值過(guò)大，模擬退火算法迭代次數(shù)會(huì)明顯增加。為此，引入基于電阻率幅度差信息的沖洗帶半徑初始值選取策略，定義幅度差系數(shù)s1～s6：

式中： Mj為 視電阻率； s1為MLR4 與MLR1 的幅度差系 數(shù)； s2為MLR4 與MLR2 的 幅 度 差 系 數(shù)； s3為MLR4 與MLR3 的幅度差系數(shù)； s4為MLR3 與MLR1的幅度差系數(shù)； s5為MLR3 與MLR2 的幅度差系數(shù)；s6為MLR2 與MLR1 的幅度差系數(shù)。

利用有限元法計(jì)算得到的s1～s6隨沖洗帶半徑的變化關(guān)系如圖1所示。

圖 1 幅度差系數(shù)隨沖洗帶半徑的變化關(guān)系Fig. 1 The relationship between the coefficient of amplitude difference and with the radius of the flushing zone

基于圖1 分析沖洗帶半徑 rxo與 s1， s2， s3， s4， s5和s6的 多元回歸關(guān)系，得到 rxo的回歸關(guān)系式：

將沖洗帶半徑多元回歸結(jié)果與理論計(jì)算值進(jìn)行相關(guān)性分析，可以看出R2高達(dá)0.984 4，說(shuō)明了多元回歸的準(zhǔn)確性（見(jiàn)圖2）。

圖 2 多元回歸結(jié)果與理論值相關(guān)性分析Fig. 2 Correlation analysis between the multiple regression results of flushing zone radius and the theoretical values

對(duì)于沖洗帶電阻率及地層電阻率，分別將微球型聚焦測(cè)井結(jié)果RMSF 及MLR4 作為其初始值。

2.1.2 新解的求取

在初始值的基礎(chǔ)上進(jìn)行新解的求取， rxo， Rxo和Rt分別按照以下公式進(jìn)行迭代更新：

式中： Tk為退火溫度； T0為初始溫度，一般數(shù)值較大； k 為迭代次數(shù)； Rtnew為地層電阻率新解， Ω·m；Rxonew為 沖洗帶電阻率新解， Ω·m ； rxonew為沖洗帶半徑新解，mm； Rtmax， Rtmin為地層電阻率的最大、最小值， Ω·m ； Rxomax， Rxomin為沖洗帶電阻率的最大、最小值， Ω·m ； rxomax， rxomin為沖洗帶半徑的最大、最小值，mm； rand 為0 到1 之間的隨機(jī)數(shù)； sign(x)為符號(hào)函數(shù)，當(dāng)x 大于等于0 時(shí)為+1，當(dāng)x 小于0 時(shí)為-1。

2.1.3 接受準(zhǔn)則

采用Metropolis 準(zhǔn)則，判斷是否接受模擬退火算法產(chǎn)生的新解[26]：

式中：p（Tk）為接受新解的概率；f（A），f（B）分別為原始解、新解對(duì)應(yīng)的系統(tǒng)能量（即適應(yīng)度）。

若新解對(duì)應(yīng)的系統(tǒng)能量 f(B)小于前一狀態(tài)系統(tǒng)能量 f(A)，則直接接受新解；反之，若新解系統(tǒng)能量f(B)( 大于或者等)于前一狀態(tài)系統(tǒng)能量 f(A)，則以的概率接受新解；接受的新解作為下一次迭代的初始值。

2.1.4 終止條件

采用溫度終止準(zhǔn)則[27]，若系統(tǒng)溫度大于預(yù)先設(shè)定的系統(tǒng)最低溫度，則繼續(xù)進(jìn)行迭代；若系統(tǒng)溫度小于或等于系統(tǒng)最低溫度，則算法終止。

2.2 混合模擬退火算法

采用模擬退火算法對(duì)式（1）進(jìn)行求解，迭代若干次（一般15 次左右），將迭代后的新解作為初始值，利用馬奎特算法進(jìn)行反演處理，即可得到三參數(shù)反演結(jié)果。

將利用模擬退火算法得到的初始值帶入正演模型，將陣列側(cè)向測(cè)井響應(yīng)在該初始值附近線性化：

式中： Mja為視電阻率響應(yīng)， Ω·m。

視電阻率擬值矩陣和測(cè)井?dāng)?shù)據(jù)矩陣分別表示為：

式中： ΔP 為未知參量矩陣； [J]為Jacobi 矩陣。

式（18）采用阻尼最小二乘法進(jìn)行求解，即可實(shí)現(xiàn)三參數(shù)反演?；旌夏M退火算法的基本步驟為：

1）設(shè)定初始溫度、最低溫度和退火因子等參數(shù)數(shù)值；

2）基于電阻率幅度差信息策略生成初始值，并計(jì)算初始值對(duì)應(yīng)的適應(yīng)度；

3）將初始值代入式（9）—式（12），進(jìn)行模擬退火操作，利用Metropolis 準(zhǔn)則判斷是否進(jìn)行個(gè)體更新；

4）利用模擬退火算法迭代一定的代數(shù)（具體根據(jù)實(shí)際問(wèn)題而定），將新解作為初始值代入馬奎特算法；

5）當(dāng)?shù)螖?shù)達(dá)到最大值或者適應(yīng)度值小于預(yù)先設(shè)定的閾值時(shí)，將此時(shí)的新解作為反演結(jié)果輸出；若不滿(mǎn)足收斂條件，則跳轉(zhuǎn)至步驟2）重新進(jìn)行初始計(jì)算。

3 混合模擬退火算法準(zhǔn)確性評(píng)價(jià)

算法的性能評(píng)價(jià)主要包括尋優(yōu)成功率、收斂平均代數(shù)、平均最優(yōu)適應(yīng)度和最優(yōu)個(gè)體進(jìn)化曲線等方面[28-29]。馬奎特算法、模擬退火算法和混合模擬退火算法的尋優(yōu)成功率、收斂平均代數(shù)及平均最優(yōu)適應(yīng)度的對(duì)比分析如表1 所示，混合模擬退火算法反演時(shí)長(zhǎng)與儀器測(cè)量時(shí)長(zhǎng)的對(duì)比結(jié)果如表2 所示。

表 1 3 種算法的性能對(duì)比Table 1 Comparison of the performances with three algorithms

表 2 混合模擬退火算法反演時(shí)長(zhǎng)與儀器測(cè)量時(shí)長(zhǎng)對(duì)比Table 2 Comparison of the inversion time of hybrid simulated annealing algorithm with the instrument measurement time

從表1 可知，混合模擬退火算法尋優(yōu)成功率遠(yuǎn)高于馬奎特算法，且收斂速度遠(yuǎn)高于模擬退火算法，不僅很好地兼顧了馬奎特算法、模擬退火算法的優(yōu)點(diǎn)，而且避免了二者的缺點(diǎn)。從表2 可以看出，混合模擬退火算法的反演速度快于儀器測(cè)量速度，滿(mǎn)足實(shí)時(shí)反演的要求。

3 種算法的最優(yōu)個(gè)體進(jìn)化曲線如圖3 所示。從圖3 可以看出，混合模擬退火算法較模擬退火算法及馬奎特算法的初始適應(yīng)度更小，說(shuō)明基于電阻率幅度差信息的初始值優(yōu)化策略的有效性。混合模擬退火算法適應(yīng)度小于0.005，反演精度高。

圖 3 算法最優(yōu)個(gè)體進(jìn)化曲線Fig. 3 Optimal individual evolution curve with three algorithms

4 應(yīng)用實(shí)例

選取W1P-7 井1 880.00～2 020.00 m 層段進(jìn)行反演處理，處理結(jié)果如圖4 所示。W1P-7 井巖電參數(shù)a=1，b=1，m=1.72，n=1.87。圖4 中：第一道為巖性曲線道；第二道為為電阻率資料道，其中RMSF 為微球型聚焦測(cè)井視電阻率，MLR1C—MLR4C 為陣列側(cè)向測(cè)井經(jīng)過(guò)井眼校正之后的視電阻率曲線。從巖性曲線道可以看出，1 982.00～1 992.00 m 層段的伽馬曲線表現(xiàn)出明顯的砂巖特征；同時(shí)，該層段的電阻率曲線也表現(xiàn)出明顯的幅度差特征，可以確定1 982.00～1 992.00 m 層段為滲透層。第三道為反演電阻率道，第四道為根據(jù)電阻率反演結(jié)果計(jì)算的含水飽和度。根據(jù)計(jì)算的含水飽和度，可以非常直觀地顯示出在處理層段的頂部、底部分別存在一油層和一水層，分別命名為W1 層和W2 層。根據(jù)試油資料，井深1 990.00 m 處累計(jì)泵抽0.83 h，泵抽地層流體達(dá)24.49 L，流體性質(zhì)為油，證明W1 層為油層；井深2 002.01 m 處累計(jì)泵抽2.70 h，泵抽地層流體達(dá)83.86 L，流體性質(zhì)為地層水，證明W2 層為水層。試油結(jié)果與反演處理結(jié)果一致，表明了混合模擬退火算法的正確性與實(shí)用性。

圖 4 W1P-7 井反演處理結(jié)果Fig. 4 Inversion processing results of Well W1P-7

5 結(jié)論與建議

1）提出基于電阻率幅度差信息初始值選取策略的混合模擬退火反演算法，理論模型與實(shí)際資料驗(yàn)證表明，該算法不僅保存了傳統(tǒng)實(shí)時(shí)反演方法的速度優(yōu)勢(shì)（單點(diǎn)反演耗時(shí)僅需0.2 s 左右），還進(jìn)一步提高了反演精度，同時(shí)三參數(shù)反演結(jié)果與試油結(jié)果相匹配。

2）以往的電阻率測(cè)井反演算法往往以單一算法為主，算法收斂性和尋優(yōu)能力不能同時(shí)得到保證，混合模擬退火算法兼顧了收斂速度和尋優(yōu)能力，反演精度更高。

3）本文的研究有一定局限性，沖洗帶半徑初值生成公式只適用于單一陣列側(cè)向測(cè)井儀器，實(shí)際應(yīng)用中應(yīng)根據(jù)陣列側(cè)向測(cè)井儀器種類(lèi)結(jié)合本文方法重新生成初始值公式。此外，支持向量機(jī)、差分進(jìn)化算法等非線性反演算法比模擬退火算法的尋優(yōu)能力更強(qiáng)，建議進(jìn)一步研究此類(lèi)算法與線性反演算法的混合反演應(yīng)用，以得到更好的反演效果。