致密砂巖儲層自適應巖石物理建模方法

張金強 劉振峰 劉喜武 韓 磊

(①中國石化石油勘探開發研究院，北京 102206；②頁巖油氣富集機理與有效開發國家重點實驗室 北京 102206；③國家能源頁巖油研發中心 北京 102206)

0 引言

地震儲層預測技術應用于致密砂巖油氣藏時，由于儲層本身的物理性質、特征(如儲層與圍巖的縱波阻抗接近或重疊區域大)，常用的疊后縱波阻抗反演技術難以取得好的預測效果[1-5]。在此情況下需要進行疊前反演，聯合應用儲層的縱、橫波阻抗、泊松比等彈性參數對儲層及流體進行預測[6-10]。在疊前反演中，最為重要和基礎的工作是通過對已鉆井的測井資料進行巖石物理建模，獲得已鉆井的橫波信息。巖石物理建模是致密砂巖儲層地震預測中的基礎，不可回避，因此得到越來越多的重視。

對砂泥巖巖石物理研究的方法可以分為兩類。一類以巖石物理測試數據為基礎，建立砂泥巖彈性參數與巖石物性、礦物組分間的關系，并最終以經驗公式的形式顯式地表達這種關系，其中以Nur等[11-12]、Krief等[13]、Han等[14]的方法為代表；另一類則是以等效介質理論為基礎，從測井及解釋數據出發，建立地層彈性參數與地層物性參數間的聯系，其中以Xu等[15-16]提出的巖石物理模型(Xu-White模型)最具代表性，至今仍被廣泛應用[17-18]。

在第一類方法中，Nur等[11-12]提出了臨界孔隙度模型，建立了巖石骨架速度與巖石基質速度和孔隙度間的關系。Krief等[13]以Raymer等[19]的數據為基礎提出干巖石的體積模量等于基質體積模量乘以一個孔隙度的冪函數，在孔隙度遠小于臨界孔隙度的條件下，該模型與臨界孔隙度模型結果相近。Han等[14]在研究中建立了縱波速度與孔隙度以及泥質含量間的關系。第一類方法由于是建立在某個地區有限的巖心測試數據基礎上，使用時必須利用本區的測試數據或測井數據對相應的經驗公式進行校驗，并調整相應的經驗系數。

第二類方法中，Xu等[15-16]使用了Voigt-Reuss-Hill(V-R-H)模型、Kuster-Toks?z(K-T)模型和微分等效介質(DEM)模型，通過測井解釋的礦物組分和孔隙度等參數重構地層的縱、橫波速度；王震宇等[20]在致密砂巖的巖石物理建模中使用自洽模型計算巖石的體積模量；劉振峰[21]則考慮致密砂巖儲層可能發育裂縫的特點，將各向異性理論引入到致密砂巖巖石物理建模中。

第二類方法由于不受測試數據的限制，巖石物理建模需要的輸入數據易于獲得，因而得到廣泛應用，成為疊前反演必需的技術環節和手段。但是這類常規巖石物理建模方法應用于致密砂巖仍面臨如何準確獲取巖石基質彈性參數的挑戰。對于眾多的巖石物理建模方法，無論使用何種模型，流程如何變化，都需要利用巖石的礦物組分估算巖石基質的彈性模量。現有的巖石物理建模方法通常是將組成巖石的礦物體積百分含量和對應的各礦物彈性模量看作已知直接代入，或是對每種礦物的彈性模量利用經驗進行選擇。這些方法在實際應用中存在許多困難。對于致密砂巖，應用這種方法更困難，其中最主要的限制是對于礦物成分的解釋。

通常測井資料解釋中將致密砂巖的組成礦物解釋為砂質成分和泥質成分兩種，并通過放射性測井等資料確定砂質成分和泥質成分的體積百分比。但是這種劃分極其粗略，砂質成分既可以由純的石英礦物組成，也可以由石英與長石等礦物組成，各種礦物的體積百分比也因不同的成因條件而千差萬別。在此條件下，試圖以單純石英礦物的體積模量替代砂質成分的體積模量必然導致結果的偏差。對于泥質成分，由于黏土礦物的復雜性，問題更為突出。因此，現有的巖石物理建模方法不適用致密砂巖。

張金強等[22]針對碳酸鹽巖儲層的某些特定情形提出了利用統計方法獲取巖石基質模量的方法，但該方法同樣不適用于致密砂巖的復雜礦物條件。為此，本文通過聯合應用臨界孔隙度模型和Gassmann方程，提出了一種自適應估算泥質與砂質成分彈性模量的方法。

1 自適應巖石物理建模方法

1.1 常規巖石物理建模方法回顧

通過對常規巖石物理建模方法流程的回顧，可以更為清楚地認識在致密砂巖儲層領域所面臨的問題與挑戰。以應用Xu-White模型進行巖石物理建模為例，其基本流程包括以下關鍵環節。

(1)巖石基質彈性模量計算。假設巖石基質由N種礦物組成，每種礦物的體積含量為Vi，則巖石基質的彈性模量可以由下式計算

(1)

式中Xm、Xi分別為巖石基質和第i種礦物的彈性模量，既可以是體積模量，也可以是剪切模量，或者由二者線性組合形成其他彈性模量。式(1)是經典的V-R-H模型，利用Voigt界限和Reuss界限的平均計算混合物的體積模量。

(2)巖石骨架彈性模量計算。在得到巖石基質的彈性模量后，利用K-T模型和微分等效介質理論，將孔隙加入巖石基質。根據K-T模型，有

(2)

(3)

式中：Kd、Km、K′分別為巖石骨架、基質和孔隙充填物的體積模量；μd、μm、μ′ 則為相應的剪切模量；φl為巖石第l種孔隙充填物的體積百分含量；αl為對應的橫縱比；F(αl)=T2(αl)-T1(αl)/3；T1、T2是αl的函數，具體形式見附錄A。在這一環節中，將孔隙充填物的體積模量與剪切模量都設為零。K-T模型的應用條件為φ/α?1，即孔隙度與孔隙橫縱比之比遠小于1。實際應用中，這一條件較難滿足，為此利用微分等效介質理論進行處理，即將每種基質的孔隙分成若干份，依次將孔隙加入基質中，每次加入一份孔隙后基質的彈性模量進行更新。

(3)利用Gassmann方程計算孔隙充填流體后飽和流體巖石的彈性模量及速度。根據Gassmann方程，流體飽和巖石的彈性模量是巖石骨架、巖石基質以及孔隙流體體積模量的函數

(4)

μs=μd

(5)

式中：Ks、μs為流體飽和巖石的體積模量和剪切模量；φ為總孔隙度；Kf為流體的體積模量。故流體飽和巖石的縱波速度和橫波速度分別為

(6)

在常規巖石物理建模的三個關鍵環節中，人們更多地關注了第二個環節。如針對微分等效介質理論對于不同孔隙添加順序結果會出現不同的缺陷，有學者提出了應用自洽理論進行建模的思路，也獲得了一定的應用。但是，在實踐中第一個環節面臨更大的挑戰。在第一個環節，首先假設組成巖石基質的礦物成分及其體積百分比是已知的，這一條件較為苛刻。實際測井數據解釋僅可以獲得泥質和砂質的百分含量，不能獲得每種礦物的百分含量，難以滿足高精度巖石物理建模要求。因為砂質成分可以由石英、長石或其他碎屑顆粒組成，黏土可以由高嶺石、云母、綠泥石等復雜礦物組成。為應對這一困難，常規巖石物理建模通常采用兩種方法。第一種方法認為砂質成分僅由石英礦物顆粒組成，并通過交會圖的方法獲得干黏土的模量參數；第二種方法，根據地質分析資料(如X射線衍射分析)確定砂質和泥質的平均礦物組成，并在整個工區內應用。第一種方法所引入的誤差是顯而易見的，第二種方法從某種程度上減少了第一種方法造成的誤差，但是對于陸相致密砂巖儲層，由于具有沉積相帶橫向變化快的特點，難以利用有限的地質分析資料確定砂質和泥質成分的礦物組成，因此仍會在建模過程中引入較大誤差。另外，實際應用中，地質分析資料也存在缺失或獲取困難等問題。

1.2 致密砂巖自適應巖石物理建模方法

假設地層在一定深度范圍內(可以是一個沉積段、小層，也可以人工劃分)穩定沉積，測井解釋的砂質成分、泥質成分分別由特定不變的礦物組分組成，這樣可以將二者等效為兩種介質。若能通過已有數據反推得到兩種等效介質的彈性模量，就解決了上文所述巖石物理建模中面臨的實際問題。具體步驟如下。

(1)巖石基質彈性模量的估計。根據臨界孔隙度模型，有

(7)

(8)

式中φc為臨界孔隙度，對于致密砂巖一般約為0.4。

利用式(6)可根據縱、橫波速度和密度數據求取流體飽和巖石的彈性模量Ks、μs，再根據式(5)可獲得巖石骨架模量μd，然后利用式(8)可得到μm。

將式(7)代入式(4)，經簡單運算，可得

Kfx2]Km-KfKs(x-φ)=0

(9)

式中x=φ/φc。式(9)是關于Km的二次方程。解式(9)可得到方程的兩個根，取其中正值為Km的解。分析式(9)中Km各次項的系數，可以看出，當孔隙度φ等于零時，各次項的系數也等于零，此時方程無解。為防止這種情況發生，對于解釋孔隙度極小的曲線段，可根據物理上的含義，令Km=Ks。

(2)砂質及泥質成分彈性模量上限、下限的求取。在得到巖石基質的彈性模量Km、μm后，結合測井解釋的泥質含量曲線，可以通過最小二乘法求取砂質及泥質成分的等效彈性模量。第一步中應用臨界孔隙度模型，而臨界孔隙度模型是以成分較純的砂巖數據為基礎得出的經驗公式，未考慮復雜組分和孔隙結構對巖石骨架模量的影響，因此會給計算引入誤差。鑒于此，分別利用Voigt上限公式和Reuss下限公式進行反演。應用Voigt上限公式求取的砂質、泥質礦物的彈性模量作為下限，應用Reuss下限公式求取的彈性模量作為上限。以應用Voigt上限公式求取砂、泥質體積模量的過程為例說明具體求解方法。

設在一定深度范圍內，有L個樣點，其基質的體積模量可表示為Km∈R1×L。相應的測井解釋泥質含量為Vsh∈R1×L，設砂質和泥質成分的體積模量分別為Ksand、Ksh，則根據Voigt上限公式，有

Km=Ksand(I-Vsh)+KshVsh

(10)

式中I為元素全部為1的向量。則求解Ksand、Ksh就轉化為標準的最小二乘問題

(11)

(3)利用最優化方法確定砂質、泥質成分的彈性模量。由于泥質、砂質成分彈性模量與巖石速度間是復雜的非線性關系，需要用全局尋優的算法進行求解。綜合考慮計算的全局收斂性和計算成本，本文選用模擬退火算法。正演采用式(1)～式(6)所述流程，同時定義目標函數為

(12)

考慮到步驟(1)、步驟(2)中應用經驗公式可能引入的計算誤差，可以對求取的砂質和泥質成分的彈性模量范圍進行適當放大，保證方法的容錯性。

(4)砂質、泥質成分彈性模量的微調。利用以上三步可以大致得到全局最優點附近的參數值，可以采用一定的策略在結果附近對反演的參數進行微調。一個較為簡便的策略是：分別保持砂質、泥質成分的縱波模量不變，小幅度調整砂質、泥質成分的體積模量和剪切模量。其中，縱波模量的定義為

(13)

按照這個策略，如果小幅增大泥質成分剪切模量(如+0.01GPa)，就要相應地減小泥質成分體積模量(-0.013GPa)，反之亦然。這個過程可以對砂質成分和泥質成分分別進行。

通過以上四步，解決了砂質、泥質成分彈性模量選取的難題。然后，可以利用常規方法完成相應的建模過程。

2 自適應巖石物理建模實際算例

利用所提方法，對鄂爾多斯盆地南緣M靶區的X1井進行巖石物理建模。主要目的層為延長組長8油層組。該段地層主要由河流—三角洲相中砂巖、細砂巖構成，交錯層理發育，夾灰黑色、灰色和黑色泥巖。自然電位曲線呈平緩波狀夾齒狀或漏斗狀低值異常，自然伽馬值整體偏高，砂巖段一般在70～100API之間；電阻率曲線呈齒狀中高阻。在砂巖段密度有降低的趨勢。圖1為X1井長8油層組的測井曲線，可以看出，長8油層組主要由砂巖構成，泥質含量不高。該井有偶極橫波測井資料，橫波時差曲線與縱波時差曲線形態相關度高。繪制該井段縱橫波速度交會圖(圖2a)，可以看出泥巖(高GR值)具有較低的縱、橫波速度。同時，在縱波速度與縱橫波速度比交會圖(圖2b)上，泥巖與砂巖的表現出一定的分異，泥巖整體表現出較高的縱橫波速度比，但存在一定程度的混疊，且規律性不強，難以進行合理的、有明確地質意義的解釋。

圖1 X1井測井曲線

圖2 X1井實測數據交會分析

應用本文提出的方法對X1井進行巖石物理建模。首先，根據式(6)，計算巖石的剪切模量和體積模量，然后利用式(5)和式(8)得到巖石基質的剪切模量，進而利用式(9)求取巖石基質的體積模量。圖3為原位巖石彈性模量(紅色曲線)與計算的基質彈性模量(藍色曲線)間的對比，可以看出，砂巖段基質彈性模量明顯大于巖石的彈性模量，而對于泥巖段，由于本區泥巖壓實強烈，孔隙度低，故基質彈性模量與巖石彈性模量非常接近；泥巖段(深度約為2033～2047m)基質的體積模量和剪切模量明顯低于相鄰砂巖段。

圖3 原位巖石與計算的基質的彈性模量對比

第二，以巖石基質的彈性模量和解釋的泥質含量曲線作為輸入，分別應用Voigt公式和Reuss公式，并利用最小二乘法求取泥質和砂質成分的彈性模量的下限和上限。在該算例中，求得的泥質成分與砂質成分彈性模量的上、下限值如表1所示。

第三，利用模擬退火算法求取砂質與泥質成分的彈性模量。考慮到步驟二中存在的誤差會導致求取的上、下限值存在誤差，在模擬退火算法的求解中，將上限上浮20%，下限下調20%，這樣放大了求解區間，增加了方法的容錯性。利用模擬退火算法得到泥質和砂質成分的體積模量分別為19.34、46.55GPa，剪切模量為9.42、25.16GPa。圖4為實測縱橫波速度曲線與本文模擬退火算法求得的彈性模量進行巖石物理建模及應用表1中的上限和下限進行巖石物理建模的結果對比。可以看出，利用模擬退火算法求得的砂、泥質彈性模量進行巖石物理建模，其精度明顯高于直接利用上、下限進行巖石物理建模的結果。

圖4 實測縱橫波速度與巖石物理建模結果對比

表1 砂質、泥質成分彈性模量上、下限

最后，對上一步求得的結果進行微調，即

(14)

式中：K0、μ0為模擬退火求得的砂、泥質體積模量和剪切模量組成的向量；n為自然數；向量b可取四種形式：(0.1，0.1)T、(-0.1，-0.1)T、(0.1，-0.1)T、(-0.1，-0.1)T，本文取第三種。在此條件下，衡量建模結果精度的目標函數會出現一定程度的降低，當n=7時，目標函數降到最低。據此得到泥質和砂質成分的體積模量分別為20.27、45.61GPa，剪切模量為8.72、25.86GPa。最終巖石物理建模成果及誤差曲線如圖5所示。由圖可見，縱波速度、橫波速度的建模結果與實測結果吻合度很高，從誤差曲線看，縱波建模精度高于橫波建模，但整體都較高，誤差基本在5%以下。盡管建模結果與實測結果整體上非常接近，但由于利用了測井曲線中誤差較小(受井眼環境影響小)的放射性測井曲線解釋的泥質含量曲線進行建模，建模結果修正了實測結果中可能存在的受環境影響的誤差，因此在建模后的交會圖上不同巖性的區分更為明顯(圖6)。

圖5 巖石物理建模最終結果與實測速度曲線的對比

在以泥質含量為色標的建模后的縱波速度與縱橫波速度比交會圖(圖6a)上，泥巖段(高泥質含量)表現為明顯的低速、高縱橫波速度比特點，且隨著泥質含量的增高，有縱橫波速度比升高、縱波速度降低的趨勢。砂巖段與泥巖段在交會圖上分異明顯，表現為中—高速度、低縱橫波速度比的特點，并且隨著孔隙度增加縱橫波速度比有降低的趨勢(圖6b)。

通過對X1井的巖石物理建模確定了泥質與砂質成分的體積模量和剪切模量等建模參數，應用這些建模參數對工區其他無橫波測井資料的井進行巖石物理建模以估計缺失的橫波測井數據。圖7為部分井巖石物理建模后的縱波速度與實測速度的對比，二者吻合很好，展示了方法的穩定性和可靠性。圖8為工區內五口井巖石物理建模后的縱波速度與縱橫波速度比的交會圖。對比圖8與圖2b、圖6可以看出，無橫波資料井建模后得到交會圖形態與實測井的形態相似，表現出明顯的規律性。泥巖(高伽馬值)表現出較高的縱橫波速度比和較低的縱波速度，并隨砂質含量增加(伽馬值變小)速度增大，縱橫波速度比降低。而砂巖段整體表現為低縱橫波速度比的特點，隨孔隙度減小速度增高，縱橫波速度比亦增高。

圖6 巖石物理建模后縱波速度與縱橫波速度比交會圖

圖7 工區內另外三口井的巖石物理建模后的縱波速度曲線(藍色)與實測曲線(紅色)的對比

圖8 工區內五口井巖石物理建模后縱波速度與縱橫波速度比交會圖

以上規律對利用地震資料進行優質儲層預測至關重要。

3 結論

在現實實踐中利用巖石物理建模方法估算缺失的橫波測井資料仍面臨諸多挑戰。對于礦物組成復雜的致密砂巖儲層，精確地估計砂、泥成分的彈性模量是成功應用巖石物理建模方法的基礎，也是實踐中不可回避的難點。本文給出了自適應的巖石物理建模方法以克服這一困難，即：通過聯立Gassman方程與臨界孔隙度模型進行巖石基質模量的估算，再通過最小二乘方法估計砂、泥質成分彈性模量的上、下限，然后通過全局尋優算法最終確定砂、泥質成分彈性模量。應用該方法對鄂爾多斯南部M工區進行巖石物理建模，結果可更清晰地展示儲層物性與儲層彈性參數間的相互關系，并方便、可靠地估算了缺失的橫波資料。

本文方法具有簡便、實用的特點，除要求一口井有橫波測井資料外，其他井只要有常規測井資料和常規解釋資料即可。

附錄A 函數T1和T2的具體表達式

在K-T模型中出現的T1和T2，是巖石基質、巖石橢球狀包含物彈性模量及橢球體橫縱比(橢球體短軸與長軸之比)的函數。具體形式為

(A-1)

(A-2)

式中

(A-3)

[f+θ-R(f-θ+2θ2)]

(A-4)

(A-5)

(A-6)

(A-7)

H6=1+A[1+f-R(f+θ)]+B(1-θ)(3-4R)

(A-8)

(A-9)

+B(1-θ)(3-4R)

(A-10)

H9=A[(R-1)f-Rθ]+Bθ(3-4R)

(A-11)

(A-12)

(A-13)

(A-14)

(A-15)

(A-16)

在應用K-T模型將孔隙充填物加入巖石基質過程中，令孔隙充填物的模量為零，這樣有A=-1、B=0。