馮周，李寧，武宏亮，劉英明，夏守姬，原野

（1.中國石油勘探開發(fā)研究院；2.北京大學(xué)地球與空間科學(xué)學(xué)院）

0 引言

縫洞型碳酸鹽巖、火山巖油氣藏是目前中國油氣勘探的重點(diǎn)領(lǐng)域，但這類儲(chǔ)集層往往呈現(xiàn)出極強(qiáng)的非均質(zhì)性與各向異性，通過測(cè)井識(shí)別地層巖性和流體性質(zhì)存在較大困難。傳統(tǒng)的POR（單孔隙度程序）、CRA（復(fù)雜巖性程序）等處理方法一般采用固定的解釋模型分步計(jì)算地層參數(shù)，計(jì)算出的地層組分類型有限，難以滿足縫洞儲(chǔ)集層測(cè)井評(píng)價(jià)的需要。

20世紀(jì)80年代初，斯倫貝謝公司提出了測(cè)井資料最優(yōu)化處理方法[1]。與傳統(tǒng)處理方法不同，最優(yōu)化處理是將所有的測(cè)井信息、誤差及地區(qū)地質(zhì)經(jīng)驗(yàn)綜合成多維信息復(fù)合體，通過數(shù)學(xué)中的最優(yōu)化技術(shù)求取滿足地層條件的最優(yōu)解釋結(jié)果，為復(fù)雜儲(chǔ)集層測(cè)井評(píng)價(jià)提供了有效途徑。隨后，眾多國內(nèi)外學(xué)者對(duì)最優(yōu)化處理方法進(jìn)行了改進(jìn)和完善，并取得了較好的地質(zhì)應(yīng)用效果[2-3]。

測(cè)井最優(yōu)化處理的基礎(chǔ)是依據(jù)地層模型建立各類測(cè)井響應(yīng)方程，這些方程表征了測(cè)井響應(yīng)與地層特性參數(shù)之間的定量關(guān)系，其精度決定了最優(yōu)化處理結(jié)果與實(shí)際地層的符合程度。現(xiàn)有的最優(yōu)化處理方法中測(cè)井響應(yīng)方程都是基于均勻?qū)訝畹貙幽Ｐ徒⒌腫4]，在常規(guī)均質(zhì)儲(chǔ)集層中能取得較好的應(yīng)用效果，但在非均質(zhì)、各向異性縫洞儲(chǔ)集層中，受巖性變化及孔隙分布影響，地層測(cè)井響應(yīng)特別是電阻率響應(yīng)變化更為復(fù)雜，計(jì)算得到的地層組分含量特別是孔隙流體組分含量存在較大偏差。因此，本文提出針對(duì)縫洞儲(chǔ)集層的測(cè)井最優(yōu)化處理新方法，并通過模擬計(jì)算和實(shí)際測(cè)井資料處理對(duì)方法可靠性進(jìn)行驗(yàn)證。

1 測(cè)井響應(yīng)方程

設(shè)目標(biāo)縫洞地層由l種骨架礦物、n種黏土礦物、孔隙水及油氣組成。儀器探測(cè)范圍內(nèi)，各骨架礦物的體積分?jǐn)?shù)為Vmi，各黏土礦物的體積分?jǐn)?shù)為Vcj（包含一定量的束縛水），孔隙水與油氣的體積分?jǐn)?shù)分別為Vw和Vhyd。

對(duì)于常規(guī)測(cè)井，其響應(yīng)方程仍可根據(jù)體積模型統(tǒng)一表述為：

對(duì)于電阻率測(cè)井，其響應(yīng)方程通常根據(jù)由各類飽和度模型描述的地層電阻率與含油氣飽和度的關(guān)系來建立，傳統(tǒng)的最優(yōu)化處理方法中以Archie公式及其改進(jìn)模型為主，不能反映復(fù)雜縫洞孔隙結(jié)構(gòu)對(duì)地層電阻率的影響。李寧[5]于1989年以非均質(zhì)各向異性地層模型為基礎(chǔ)，通過完整的數(shù)學(xué)推導(dǎo)，給出了地層電阻增大率與含水飽和度之間的一般關(guān)系式：

在大量巖電實(shí)驗(yàn)測(cè)量和數(shù)據(jù)分析的基礎(chǔ)上，李寧[5]給出了（2）式在縫洞儲(chǔ)集層含水飽和度評(píng)價(jià)時(shí)的最佳截短形式：

（2）式表達(dá)了非均質(zhì)復(fù)雜儲(chǔ)集層電阻率與含水飽和度的精確定量關(guān)系，但在現(xiàn)有技術(shù)條件下，確定（2）式中的所有待定參數(shù)十分困難。（3）式用盡量少的待定參數(shù)和盡量簡單的表達(dá)形式獲得了盡可能高的計(jì)算精度，碳酸鹽巖、火山巖全直徑巖心實(shí)驗(yàn)及油田應(yīng)用證實(shí)了該式的正確性[6-8]。

根據(jù)地層電阻率與地層因素、電阻增大率的一般關(guān)系式，電阻率測(cè)井響應(yīng)方程可表示為：

將（3）式代入（4）式，則非均質(zhì)、各向異性縫洞儲(chǔ)集層原狀地層電阻率可表示為：

同理，沖洗帶地層電阻率可表示為：

（5）式中地層總孔隙度φ、含水飽和度Sw可用地層流體組分含量來表示：

2 測(cè)井最優(yōu)化問題求解

（5）式—（8）式為縫洞儲(chǔ)集層電阻率測(cè)井及其他常規(guī)測(cè)井響應(yīng)方程，其中各組分的響應(yīng)參數(shù)可通過理論計(jì)算或?qū)嶒?yàn)方法確定[9-11]。對(duì)給定的地層解釋模型，可由測(cè)井響應(yīng)方程列出各測(cè)井曲線的理論測(cè)井響應(yīng)值：

由最小二乘理論可建立測(cè)井最優(yōu)化目標(biāo)函數(shù)：

同時(shí)，地層組分還應(yīng)滿足一定的巖石物理及地質(zhì)約束條件，約束條件方程組一般形式可記為：

由（5）式、（6）式可見，理論電阻率響應(yīng)是關(guān)于地層組分含量的復(fù)雜非線性函數(shù)，則求解最優(yōu)化目標(biāo)函數(shù)取最小值時(shí)的v的問題是一個(gè)典型的帶約束條件的非線性最小二乘問題。以往的測(cè)井資料最優(yōu)化處理方法研究中，研究者對(duì)測(cè)井最優(yōu)化問題的求解算法進(jìn)行了一些探討與分析[12-15]，但主要針對(duì)線性最優(yōu)化問題，對(duì)非線性最優(yōu)化問題一般采用多目標(biāo)最優(yōu)化問題的通解算法，沒有針對(duì)測(cè)井最優(yōu)化目標(biāo)函數(shù)最小二乘形式的特點(diǎn)求解，因而求解效率和計(jì)算精度有限。因此，本文結(jié)合懲罰函數(shù)法和 Levenberg-Marquardt算法[16]對(duì)目標(biāo)函數(shù)進(jìn)行求解。

首先，利用懲罰函數(shù)法，將約束問題轉(zhuǎn)化為無約束問題，無約束問題目標(biāo)函數(shù)為：

當(dāng)v滿足約束條件（hr(v)小于等于0）時(shí)，懲罰項(xiàng)M P(v)為零；當(dāng)v不滿足約束條件時(shí)，M P(v)大于0，且隨M的增大而增大，當(dāng)M趨于無窮大時(shí)，無約束問題的最優(yōu)解v*將充分趨于約束區(qū)域的邊界，此時(shí)v*也是滿足約束條件的最優(yōu)解。

針對(duì)約束問題目標(biāo)函數(shù)的形式，懲罰函數(shù)可設(shè)為約束條件的平方和形式：

根據(jù)（10）式—（13）式，無約束問題目標(biāo)函數(shù)可表示為：

其中

對(duì)無約束非線性最小二乘問題，可利用Levenberg-Marquardt算法進(jìn)行求解，其每步的迭代增量可通過如下的正則方程組來計(jì)算：

（15）式中，阻尼因子μ反映當(dāng)前迭代點(diǎn)目標(biāo)函數(shù)實(shí)際下降量與 1階近似下降量的差異，在計(jì)算中自適應(yīng)調(diào)整。求解（15）式，即可得到每步的迭代增量，當(dāng)?shù)鉂M足精度要求且滿足（11）式時(shí)，即可認(rèn)為求得最優(yōu)解。

當(dāng)阻尼因子較小時(shí)，（15）式即為Gauss-Newton法最優(yōu)步長計(jì)算式；當(dāng)阻尼因子很大時(shí)，（15）式則為蛻化梯度下降法的最優(yōu)步長計(jì)算式。因此，本文提出的算法同時(shí)具有Gauss-Newton法的局部收斂性和梯度下降法的全局特性，具有較高的計(jì)算效率和精度。此外，該算法在計(jì)算過程中只需要求解目標(biāo)函數(shù)的1階導(dǎo)，運(yùn)算量較小，適用于測(cè)井響應(yīng)方程中復(fù)雜電阻率響應(yīng)方程形式。

3 模擬計(jì)算

為了進(jìn)一步驗(yàn)證本文提出的縫洞儲(chǔ)集層測(cè)井最優(yōu)化處理方法的穩(wěn)定性和計(jì)算精度，以縫洞型碳酸鹽巖儲(chǔ)集層中常見的灰?guī)r、云質(zhì)灰?guī)r地層模型（各地層組分含量見表1）為例進(jìn)行了模擬計(jì)算。根據(jù)各類測(cè)井曲線的響應(yīng)方程及響應(yīng)參數(shù)，可計(jì)算得到地層模型的理論測(cè)井響應(yīng)值，然后將理論測(cè)井響應(yīng)值及相應(yīng)的響應(yīng)參數(shù)值作為最優(yōu)化處理輸入信息，對(duì)每類地層模型分別進(jìn)行1 000組隨機(jī)初始值的計(jì)算。

圖1、圖2分別為灰?guī)r、云質(zhì)灰?guī)r地層模型地層組分含量的計(jì)算結(jié)果，結(jié)合表 1中的數(shù)據(jù)進(jìn)行分析可以發(fā)現(xiàn)：本文提出的最優(yōu)化處理方法受初始值影響較小，對(duì)于隨機(jī)選取的1 000組初始值，均能計(jì)算得到最優(yōu)解，計(jì)算結(jié)果與模型值基本一致，能夠滿足測(cè)井解釋的需要。

表1 灰?guī)r、云質(zhì)灰?guī)r地層模型各地層組分含量

圖1 灰?guī)r地層模型地層組分含量計(jì)算結(jié)果

圖2 云質(zhì)灰?guī)r地層模型地層組分含量計(jì)算結(jié)果

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

將本文提出的縫洞儲(chǔ)集層測(cè)井最優(yōu)化處理方法編寫成處理模塊，集成在中國石油新一代測(cè)井軟件CIFLog1.0中。利用該處理模塊對(duì)實(shí)際測(cè)井資料進(jìn)行了處理，并與傳統(tǒng)最優(yōu)化處理方法計(jì)算結(jié)果及巖心定量分析結(jié)果進(jìn)行了對(duì)比。以某碳酸鹽巖儲(chǔ)集層 L井測(cè)井資料處理為例進(jìn)行具體說明。

L井所在儲(chǔ)集層深度為6 052.6～6 077.4 m，巖性以灰?guī)r為主，取心巖心描述和物性分析表明該段儲(chǔ)集層溶蝕孔洞發(fā)育。為了進(jìn)一步確定儲(chǔ)集層的巖石物理特性，對(duì)該段儲(chǔ)集層全直徑巖心（見圖3a）進(jìn)行了驅(qū)替實(shí)驗(yàn)，并分別利用Archie公式與最優(yōu)飽和度方程對(duì)Sw-I實(shí)驗(yàn)測(cè)量結(jié)果進(jìn)行了擬合，擬合曲線見圖3b。由圖3b可知，在含水飽和度大于50%時(shí)，Archie公式擬合結(jié)果與最優(yōu)飽和度方程擬合結(jié)果比較接近，但在含水飽和度較低時(shí)，最優(yōu)飽和度方程擬合效果較好。

圖3 L井X號(hào)巖心照片及Sw-I擬合曲線

基于Sw-I方程擬合結(jié)果，分別使用傳統(tǒng)最優(yōu)化處理方法與本文方法對(duì) L井測(cè)井資料進(jìn)行處理，其中傳統(tǒng)最優(yōu)化處理方法采用Archie公式擬合結(jié)果，本文方法采用最優(yōu)飽和度方程擬合結(jié)果，其他地層響應(yīng)參數(shù)均一致。

圖 4為利用兩種最優(yōu)化處理方法計(jì)算得到的地層方解石、白云石含量與巖心分析結(jié)果，圖 5為采用兩種最優(yōu)化處理方法計(jì)算得到的白云石含量與巖心分析結(jié)果間誤差分布情況。由圖4、圖5可知，本文提出的縫洞儲(chǔ)集層最優(yōu)化處理方法的計(jì)算結(jié)果與巖心分析結(jié)果更為接近，計(jì)算精度更高。

圖4 兩種最優(yōu)化處理方法得到的礦物含量與巖心分析結(jié)果對(duì)比

圖5 兩種最優(yōu)化處理方法得到的白云石含量與巖心分析結(jié)果間誤差分布

將利用兩種最優(yōu)化處理方法計(jì)算得到的地層流體組分含量與巖心分析結(jié)果（見圖6）進(jìn)行對(duì)比后發(fā)現(xiàn)：在儲(chǔ)集層段，采用本文方法處理得到的含油氣飽和度與密閉取心結(jié)果吻合度較好，采用傳統(tǒng)方法處理得到的含油氣飽和度與密閉取心結(jié)果相比偏高，兩種方法的含油氣飽和度相差約10%。以儲(chǔ)集層有效厚度20 m、孔隙度6%計(jì)算，假設(shè)油氣藏含油氣面積為10 km2，含油氣飽和度相差10%相當(dāng)于油氣儲(chǔ)量相差1.2×106m3，會(huì)對(duì)整個(gè)油氣藏資源評(píng)價(jià)及開發(fā)方案設(shè)計(jì)產(chǎn)生重要影響。圖 7為采用兩種最優(yōu)化處理方法計(jì)算得到的含油氣飽和度與巖心分析結(jié)果間誤差分布情況，可以看出，本文方法的計(jì)算結(jié)果與巖心分析結(jié)果更接近。

圖6 兩種最優(yōu)化處理方法得到的流體組分含量與巖心分析結(jié)果對(duì)比

圖7 兩種最優(yōu)化處理方法得到的含油氣飽和度與密閉取心結(jié)果間誤差分布

5 結(jié)論

本文提出了針對(duì)非均質(zhì)、各向異性縫洞儲(chǔ)集層的測(cè)井最優(yōu)化處理新方法。以非均質(zhì)、各向異性縫洞儲(chǔ)集層模型為基礎(chǔ)，明確了各測(cè)井曲線響應(yīng)方程形式，最大程度上表達(dá)了縫洞儲(chǔ)集層地層響應(yīng)特征。建立了縫洞儲(chǔ)集層測(cè)井最優(yōu)化目標(biāo)函數(shù)，并結(jié)合懲罰函數(shù)法和Levenberg-Marquardt算法求解目標(biāo)函數(shù)，算法具有計(jì)算量小、求解效率高的特點(diǎn)。

以灰?guī)r、云質(zhì)灰?guī)r地層模型為例進(jìn)行了模擬計(jì)算。結(jié)果表明：本文提出的最優(yōu)化處理方法受初始值影響較小，對(duì)于隨機(jī)選取的1 000組初始值，均能計(jì)算得到最優(yōu)解且計(jì)算精度高，適用于測(cè)井資料評(píng)價(jià)與解釋。

采用本文提出的縫洞儲(chǔ)集層測(cè)井最優(yōu)化處理方法與傳統(tǒng)最優(yōu)化處理方法分別進(jìn)行了實(shí)際測(cè)井資料處理，并與巖心分析結(jié)果進(jìn)行了對(duì)比。結(jié)果表明：本文方法計(jì)算精度更高，計(jì)算結(jié)果與巖心分析結(jié)果較接近。

符號(hào)注釋：

Vmi——各骨架礦物的體積分?jǐn)?shù)，%；Vcj——各黏土礦物的體積分?jǐn)?shù)，%；Vw，Vhyd——孔隙水與油氣的體積分?jǐn)?shù)，%；Llog——理論測(cè)井響應(yīng)；Pmi——各骨架礦物的測(cè)井響應(yīng)參數(shù)；Pcj——各黏土礦物的測(cè)井響應(yīng)參數(shù)；Pw，Phyd——孔隙水與油氣的測(cè)井響應(yīng)參數(shù)；l——骨架礦物種類；n——黏土礦物種類；I——地層電阻增大率；Sw——含水飽和度，%；nb——組合薄片數(shù)量；lq——第q個(gè)薄片上導(dǎo)電礦物的單元總數(shù)；pq，hqk，θqk，b1，b2，b3——待定系數(shù)；F——地層因素；a——巖性系數(shù)；m——膠結(jié)指數(shù)；φ——地層總孔隙度，%；φcj——各黏土礦物的孔隙度，%；Rt——地層電阻率，?·m；Rxo——沖洗帶地層電阻率，?·m；Rw——地層水電阻率，?·m；Rmf——鉆井液濾液電阻率，?·m；Sxo——沖洗帶含水飽和度，%；tcs——各測(cè)井曲線理論測(cè)井響應(yīng)值；fs(v) ——根據(jù)地層解釋模型建立的各測(cè)井曲線測(cè)井響應(yīng)方程；v——各地層組分含量組成的向量；F(v) ——測(cè)井最優(yōu)化目標(biāo)函數(shù)；tms——各測(cè)井曲線實(shí)際測(cè)井響應(yīng)值；ws——不同測(cè)井方法標(biāo)準(zhǔn)化處理因子；nl——測(cè)井曲線數(shù)目；hr(v) ——地層組分約束條件方程；nc——約束條件數(shù)目；cr——約束條件系數(shù)矩陣C的第r行；br——約束條件邊界；P(v) ——懲罰函數(shù)；M——懲罰因子；Φ(v，M) ——無約束問題目標(biāo)函數(shù)；J——矩陣R的Jacobi矩陣；E——單位矩陣；μ——阻尼因子；h——每步的迭代增量；Pe——光電吸收截面指數(shù)，b/e；GR——自然伽馬，API；Rlls——淺側(cè)向電阻率，?·m；Rlld——深側(cè)向電阻率，?·m；ρcore——巖心分析密度，g/cm3；Δt——聲波時(shí)差，μs/m；ρ——巖性密度，g/cm3；φN——中子孔隙度，%；Wc——計(jì)算方解石含量（干重），%；Wcc——巖心分析方解石含量（干重），%；Wd——計(jì)算白云石含量（干重），%；Wcd——巖心分析白云石含量（干重），%；φc——巖心分析孔隙度，%；Sc——密閉取心飽和度，%；So——計(jì)算飽和度，%。

[1]Mayer C.GLOBAL, a new approach to computer-processed log interpretation[R].SPE 9341-MS, 1980.

[2]Quirein J, Kimminau S, La Vigne J, et al.A coherent framework for developing and applying multiple formation evaluation models[C]//Transactions of the SPWLA 27th Annual Logging Symposium.Houston: SPWLA, 1986.

[3]馮國慶, 陳軍, 張烈輝, 等.最優(yōu)化測(cè)井解釋的遺傳算法實(shí)現(xiàn)[J].天然氣工業(yè), 2002, 22(6): 48-51.Feng Guoqing, Chen Jun, Zhang Liehui, et al.Optimum logging interpretation realizing based on genetic algorithm[J].Natural Gas Industry, 2002, 22(6): 48-51.

[4]鄧少貴, 莫宣學(xué), 盧春利, 等.縫-洞型地層縫洞的雙側(cè)向測(cè)井響應(yīng)數(shù)值模擬[J].石油勘探與開發(fā), 2012, 39(6): 706-712.Deng Shaogui, Mo Xuanxue, Lu Chunli, et al.Numerical simulation of the dual laterolog response to fractures and caves in fracturedcavernous formation[J].Petroleum Exploration and Development,2012, 39(6): 706-712.

[5]李寧.電阻率-孔隙度、電阻率-含油(氣)飽和度關(guān)系的一般形式及其最佳逼近函數(shù)類型的確定(I)[J].地球物理學(xué)報(bào), 1989, 32(5): 580-591.Li Ning.General forms of the resistivity-porosity and resistivityoil/gas saturation relations, as well as the determination of their optimum approximating function types (I)[J].Chinese Journal of Geophysics, 1989, 32(5): 580-591.

[6]李寧, 陶宏根, 劉傳平.酸性火山巖測(cè)井解釋理論方法與應(yīng)用[M].北京: 石油工業(yè)出版社, 2009.Li Ning, Tao Honggen, Liu Chuanping.Theoretical method of acidic volcanic logging interpretation and its application[M].Beijing:Petroleum Industry Press, 2009.

[7]李寧.中國海相碳酸鹽巖測(cè)井解釋概論[M].北京: 科學(xué)出版社,2013.Li Ning.Introduction to logging interpretation of marine carbonate in China[M].Beijing: Science Press, 2013.

[8]邵先杰, 孫玉波, 孫景民, 等.煤巖參數(shù)測(cè)井解釋方法: 以韓城礦區(qū)為例[J].石油勘探與開發(fā), 2013, 40(5): 559-565.Shao Xianjie, Sun Yubo, Sun Jingmin, et al.Logging interpretation of coal petrologic parameters: A case study of Hancheng mining area[J].Petroleum Exploration and Development, 2013, 40(5): 559-565.

[9]高楚橋, 譚廷棟.常見測(cè)井響應(yīng)參數(shù)的理論計(jì)算[J].石油地球物理勘探, 1997, 32(6): 818-825.Gao Chuqiao, Tan Tingdong.Theoretical computation of ordinary parameters of logging responses[J].Oil Geophysical Prospecting,1997, 32(6): 818-825.

[10]斯倫貝謝公司.測(cè)井解釋常用巖石礦物手冊(cè)[M].吳慶巖, 張愛軍,譯.北京: 石油工業(yè)出版社, 1998.Schlumberger.Handbook of common rock minerals in logging interpretation[M].Wu Qingyan, Zhang Aijun, Trans.Beijing: Petroleum Industry Press, 1998.

[11]寧伏龍, 吳能友, 李實(shí), 等.基于常規(guī)測(cè)井方法估算原位水合物儲(chǔ)集層力學(xué)參數(shù)[J].石油勘探與開發(fā), 2013, 40(4): 507-512.Ning Fulong, Wu Nengyou, Li Shi, et al.Estimation of in-situ mechanical properties of gas hydrate-bearing sediments by well logging[J].Petroleum Exploration and Development, 2013, 40(4): 507-512.

[12]肖立志, 鐘興水.GLOBAL測(cè)井解釋方法Incoherence函數(shù)性質(zhì)研究[J].石油學(xué)報(bào), 1990, 11(2): 49-57.Xiao Lizhi, Zhong Xingshui.Study on the incoherence function in GLOBAL program[J].Acta Petrolei Sinica, 1990, 11(2): 49-57.

[13]高楚橋, 鐘興水.計(jì)算巖石礦物成分初始值的極值函數(shù)法[J].江漢石油學(xué)院學(xué)報(bào), 1992, 14(4): 25-30.Gao Chuqiao, Zhong Xingshui.A minimization method to estimate the initial value of matrix content[J].Journal of Jianghan Petroleum Institute, 1992, 14(4): 25-30.

[14]雍世和, 孫建孟.測(cè)井?dāng)?shù)字處理中最優(yōu)化方法的選擇[J].石油大學(xué)學(xué)報(bào): 自然科學(xué)版, 1988, 12(4/5): 11-26.Yong Shihe, Sun Jianmeng.Selection of the optimization method for digital well logging data processing[J].Journal of the University of Petroleum, China: Edition of Natural Science, 1988, 12(4/5): 11-26.

[15]謝剛, 胡振平, 羅利, 等.基于約束最小二乘理論的復(fù)雜巖性測(cè)井識(shí)別方法[J].測(cè)井技術(shù), 2007, 31(4): 354-356.Xie Gang, Hu Zhenping, Luo Li, et al.Logging identification of complicated lithology with constrained least square algorithm[J].Well Logging Technology, 2007, 31(4): 354-356.

[16]Madsen K, Nielsen H B, Tingleff O.Methods for non-linear least squares problems[M].Kgs.Lyngby, Copenhagen: Technical University of Denmark, 2004.