基于信賴域算法的試井曲線自動(dòng)擬合研究

王一 陳殿遠(yuǎn) 韓鑫 鄭志鋒 張騫

(中海石油(中國(guó))有限公司湛江分公司， 廣東 湛江 524057)

試井自動(dòng)擬合分析算法本質(zhì)上是一個(gè)解決地層或油藏未知參數(shù)的非線性函數(shù)問(wèn)題，即采用試井分析理論模型去匹配實(shí)際測(cè)試的壓力數(shù)據(jù)，從而求解最佳地層參數(shù)。國(guó)內(nèi)外學(xué)者一般是將通過(guò)試井解釋求解地層參數(shù)的問(wèn)題轉(zhuǎn)化為非線性最小二乘法問(wèn)題[1-2]。對(duì)于非線性最小二乘法問(wèn)題的求解，目前一般采用Newton算法(牛頓算法)和Gauss - Newton法(高斯-牛頓算法)[3-6]。同時(shí)，在此基礎(chǔ)上不斷優(yōu)化并發(fā)展了一些新算法，如Levenberg-Marquardt算法(萊文貝格-馬夸特算法，簡(jiǎn)稱LM算法)、Levenberg-Marquardt-Fletcher法(萊文貝格-馬夸特-弗萊徹算法，簡(jiǎn)稱LMF算法)[7]、QNew算法(擬牛頓算法)、FM算法(快速擬合算法)等，以及粒子群算法、遺傳算法[8-9]等智能算法。運(yùn)用這些算法時(shí)會(huì)遇到很多實(shí)際問(wèn)題，如目標(biāo)函數(shù)具發(fā)散性，或有些算法擬合過(guò)程中耗時(shí)過(guò)長(zhǎng)而無(wú)法滿足軟件編制的要求。為了避免此類問(wèn)題，研究人員進(jìn)行了多種嘗試。本次研究將嘗試在試井曲線擬合中引入一種自動(dòng)擬合算法 —— 信賴域算法。

1 信賴域算法理論

信賴域算法是一種求解非線性優(yōu)化問(wèn)題的數(shù)值方法，也是一種迭代算法，即從給定的初始解出發(fā)，通過(guò)逐步迭代，不斷改進(jìn)，直到獲得滿意的近似最優(yōu)解為止[10]。

首先，在每次迭代中，給出一個(gè)信賴域，這個(gè)信賴域一般是當(dāng)前迭代點(diǎn)xk的一個(gè)小鄰域；然后，在這個(gè)鄰域內(nèi)求解一個(gè)子問(wèn)題，得到試探步長(zhǎng)Sk；接著，通過(guò)某個(gè)評(píng)價(jià)函數(shù)來(lái)判斷是否接受該試探步，并確定下一次迭代的信賴域。如果試探步長(zhǎng)被接受，則有，xk+1=xk+Sk；否則，xk+1=xk。

新的信賴域大小取決于試探步長(zhǎng)的優(yōu)劣。如果試探步長(zhǎng)較優(yōu)，就在下一步使信賴域擴(kuò)大或保持不變；否則，應(yīng)縮小信賴域。定義當(dāng)前迭代點(diǎn)xk的鄰域定義為：

(1)

式中：Ωk為當(dāng)前迭代點(diǎn)xk的鄰域；Rn為對(duì)稱正定矩陣；Δk為第k步的信賴域半徑。

目標(biāo)函數(shù)在極值點(diǎn)附近，近似于1個(gè)二次函數(shù)。對(duì)于無(wú)約束優(yōu)化問(wèn)題，利用二次逼近，構(gòu)造式(2)所示信賴域子問(wèn)題：

(2)

式中：s=x-xk；gk是目標(biāo)函數(shù)f(x)在當(dāng)前迭代點(diǎn)xk處的梯度；Bk∈Rn×n對(duì)稱，是f(x)在xk處Hesse陣2f(xk)或者其近似值；q(k)(s)是目標(biāo)函數(shù)f(x)的二次逼近式，即二次模型函數(shù)。

假設(shè)Sk是式(2)的解，則目標(biāo)函數(shù)f(x)在第k步的實(shí)際下降量(真實(shí)下降量)Aredk為：

Aredk=f(xk)-f(xk+Sk)

(3)

二次模型函數(shù)q(k)(s)的下降量(預(yù)測(cè)下降量)Predk為：

Predk=q(k)(0)-q(k)(Sk)

(4)

定義比值rk：

(5)

rk衡量了二次模型與目標(biāo)函數(shù)的逼近程度，其值越趨近于1，表示接近程度越好。因此，我們也可以用rk來(lái)確定下次迭代的信賴域半徑。

2 信賴域算法實(shí)現(xiàn)步驟

信賴域算法的實(shí)現(xiàn)可分以下6個(gè)步驟。

Step2如果‖gk‖≤ε，則停止。

Step3近似求解子問(wèn)題，即式(2)，得到Sk。

Step4計(jì)算f(xk+Sk)和rk，令

Step5校正信賴域半徑，令

Δk+1∈(0,γ1Δk]，ifrk<η1

Δk+1∈[γ1Δk,Δk]，ifrk∈[η1，η2)

Step6產(chǎn)生Bk+1，校正q(k)，令k∶=k+1，轉(zhuǎn)Step 2。

3 算法代碼及效果分析

3.1 部分關(guān)鍵代碼

根據(jù)信賴域方法原理及實(shí)現(xiàn)步驟，編制相應(yīng)的試井曲線自動(dòng)擬合模塊，部分關(guān)鍵代碼如下：

public static List Auto_fitting(InArray data1, InArray data2, List> CSK1_parameter1_duration1_rate2)

{

∥從CSK1_parameter1_duration1_rate2中分別提取擬合參數(shù)CSK、常規(guī)參數(shù)、產(chǎn)量持續(xù)時(shí)間、關(guān)井前各個(gè)產(chǎn)量、擬壓力歷史數(shù)據(jù)恢復(fù)數(shù)據(jù)

Array CSKR = CSK1_parameter1_duration1_rate2[0]；

…∥此處符號(hào)代表省略部分代碼

…

∥從data1中提取時(shí)間數(shù)據(jù)和壓差數(shù)據(jù)

Arrayx1=data1[full, 0]；

Arrayy1=data1[full, 1]；

∥從data2中提取擬壓力變換表述

Array pressure1 =data2[full, 0]；

Array pseudopressure1 =data2[full, 1]；

…∥此處符號(hào)代表省略部分代碼

…

watch.Start()；∥開(kāi)始監(jiān)視代碼運(yùn)行時(shí)間

∥調(diào)用matlab部分

MLApp.MLApp matlab = ToolFunc.creatematlab()；

∥把對(duì)應(yīng)數(shù)據(jù)傳輸?shù)絤atlab運(yùn)算平臺(tái)

matlab.PutWorkspaceData(“x1”, “base”, x0)；

matlab.PutWorkspaceData(“y1”, “base”, y0)；

matlab.PutWorkspaceData(“pressure1”, “base”, pressure0)；

matlab.PutWorkspaceData(“pseudopressure1”, “base”, pseudopressure0)；

matlab.PutWorkspaceData(“para”, “base”, CSK1_parameter1_duration1_rate1_tmp)；

matlab.PutWorkspaceData(“CSKR”, “base”, CSKR1)；

matlab.PutWorkspaceData(“duration0”, “base”, duration0)；

matlab.PutWorkspaceData(“rate0”, “base”, rate0)；

matlab.PutWorkspaceData(“boud0”, “base”, boud1)；

matlab.PutWorkspaceData(“rate_times0”, “base”, rate_times0)；

watch.Stop()；∥停止監(jiān)視

TimeSpan timespan = watch.Elapsed；∥獲取當(dāng)前實(shí)例測(cè)量得出的總時(shí)間

System.Diagnostics.Debug.WriteLine(“執(zhí)行時(shí)間:{0}(s )”,timespan.TotalMilliseconds)；

∥執(zhí)行matlab的擬合函數(shù)

watch.Restart()；

∥針對(duì)不同試井模型，調(diào)用對(duì)應(yīng)的模型方法matlab.Execute(@“[bestCSK,resnorm]=Closed_Finite_Circle(x1,y1,pressure1,pseudopressure1,CSKR,para,duration0,rate0,rate_times0,boud0)”)；

watch.Stop()；∥停止監(jiān)視

timespan =watch.Elapsed；∥獲取當(dāng)前實(shí)例測(cè)量得出的總時(shí)間

System.Diagnostics.Debug.WriteLine(“函數(shù)擬合時(shí)間:{0}(毫秒)”, timespan.TotalMilliseconds)；

BestCSK = matlab.GetVariable(“bestCSK”, “base”)；

…∥此處符號(hào)代表省略部分代碼

…

}

3.2 效果分析

應(yīng)用均質(zhì)無(wú)限大、圓形封閉邊界理論模型對(duì)信賴域算法進(jìn)行效果分析，并與解決非線性最小二乘法問(wèn)題常用的LM算法進(jìn)行對(duì)比。

(1) 均質(zhì)無(wú)限大油藏模型分析。某均質(zhì)無(wú)限大油藏，原始地層壓力為30 MPa，其中某油井在無(wú)污染條件下以100 m3d的生產(chǎn)水平穩(wěn)定生產(chǎn)6 h。其他基本參數(shù)有：儲(chǔ)層有效厚度，10 m；體積系數(shù)，1.0；孔隙度，0.2；黏度，0.5 mPa·s；綜合壓縮系數(shù)，4×10-4MPa-1；井半徑，0.1 m。油井壓力與時(shí)間數(shù)據(jù)如表1所示。

表1 某油井壓力、時(shí)間數(shù)據(jù)

針對(duì)待擬合參數(shù)井筒儲(chǔ)集系數(shù)(C)、 表皮系數(shù)(S)、滲透率(K)，分別采用LM法和信賴域法，以壓力誤差平方和最小為目標(biāo)函數(shù)，對(duì)壓力降落雙對(duì)數(shù)曲線進(jìn)行自動(dòng)擬合。這2種方法獲得的最佳擬合參數(shù)差異非常小(見(jiàn)圖1、圖2)，其數(shù)值均在C=0.100 1 m3MPa，S=0.001，K=98.7×10-3μm2附近。

圖1 LM算法自動(dòng)擬合效果圖

(2) 均質(zhì)圓形封閉邊界模型分析。某封閉有限圓氣藏原始地層壓力為34.473 8 MPa，其中某氣井以200×104m3d的生產(chǎn)水平穩(wěn)定生產(chǎn)。其他基本參數(shù)有：儲(chǔ)層有效厚度，30 m；體積系數(shù)，0.003 817 4，孔隙度，0.2；黏度，0.026 528 mPa·s；綜合壓縮系數(shù)，0.019 192 MPa-1；井半徑，0.091 44 m；地層溫度，373.15 K。

針對(duì)需擬合參數(shù)井筒儲(chǔ)集系數(shù)(C)、表皮系數(shù)(S)、滲透率(K)、封閉有限圓半徑(R)，采用LM法和信賴域法進(jìn)行自動(dòng)擬合。這2種方法擬合所得的結(jié)果基本一致(見(jiàn)圖3、圖4)，其數(shù)值均在C=1.1 m3MPa，S=3，K=486.1×10-3μm2，R=1 512.7 m 附近。

圖2 信賴域算法自動(dòng)擬合效果圖

圖3 LM算法自動(dòng)擬合效果圖

圖4 信賴域算法自動(dòng)擬合效果圖

對(duì)比2個(gè)模型的計(jì)算結(jié)果(見(jiàn)表2)。這個(gè)模型的擬合誤差，即擬合計(jì)算壓力值與實(shí)測(cè)壓力值的相對(duì)誤差，為10-3。在此擬合精度下的試井解釋結(jié)果，可滿足生產(chǎn)解釋的要求，解釋結(jié)果與商業(yè)軟件解釋結(jié)果基本一致。但在運(yùn)行效率方面，LM算法耗時(shí)較長(zhǎng)，信賴域算法的效率明顯較高。

表2 兩種模型的算法結(jié)果對(duì)比

4 結(jié) 語(yǔ)

本次研究中，在試井曲線擬合時(shí)采用了信賴域算法。運(yùn)用信賴域算法，可對(duì)試井曲線進(jìn)行自動(dòng)擬合算法，實(shí)現(xiàn)對(duì)地層未知參數(shù)的最優(yōu)化求解。通過(guò)實(shí)例驗(yàn)證，認(rèn)為信賴域算法能夠在較快時(shí)間內(nèi)完成自動(dòng)擬合，且擬合誤差低于0.005%，能滿足試井解釋分析中參數(shù)求解的要求。