基于KL-E 的地質(zhì)力學(xué)模型參數(shù)反演及應(yīng)用

覃建華，楊琨，丁藝，張博寧，唐慧瑩

1.中國(guó)石油新疆油田分公司勘探開(kāi)發(fā)研究院，新疆 克拉瑪依 834000；2.中國(guó)石油新疆油田分公司開(kāi)發(fā)公司，新疆 克拉瑪依 834000；3.成都北方石油勘探開(kāi)發(fā)技術(shù)有限公司，四川 成都 610051；4.油氣藏地質(zhì)及開(kāi)發(fā)工程國(guó)家重點(diǎn)實(shí)驗(yàn)室·西南石油大學(xué)，四川 成都610500

引言

地應(yīng)力的分布特征對(duì)于了解油氣運(yùn)移、鉆井方案及水力壓裂設(shè)計(jì)、提出油氣藏生產(chǎn)的優(yōu)化設(shè)計(jì)方案具有重要意義[1-4]。要準(zhǔn)確預(yù)測(cè)地應(yīng)力值以及三維地應(yīng)力場(chǎng)的分布，就有必要明確儲(chǔ)層的地質(zhì)力學(xué)性質(zhì)、幾何特征及邊界條件[5]。

目前，已有多種方法可用于儲(chǔ)層地應(yīng)力的確定，例如，井壁崩塌法[6]、水力壓裂法[7]、聲發(fā)射法（AE）[8]以及震源分析法[9-10]等。此外，還可以根據(jù)經(jīng)驗(yàn)公式利用測(cè)井?dāng)?shù)據(jù)計(jì)算地應(yīng)力[11]。巖芯分析是目前確定應(yīng)力大小最直接、可靠的方法[12]，通常被作為硬數(shù)據(jù)對(duì)其他模型進(jìn)行標(biāo)定。但受費(fèi)用、測(cè)試樣品的干擾性及采樣點(diǎn)位置等影響，使得通過(guò)巖芯分析來(lái)描述整個(gè)儲(chǔ)層特征具有較大的挑戰(zhàn)性[13]。因此，有必要采用其他方法來(lái)獲得整個(gè)地應(yīng)力場(chǎng)的應(yīng)力分布。

側(cè)壓系數(shù)法假設(shè)地應(yīng)力水平分量與垂直分量之比是常數(shù)，同時(shí)考慮應(yīng)力隨深度的變化，但該方法得到的應(yīng)力分布評(píng)價(jià)較為粗略[14-15]。Gens 等[16]和Miranda 等[17]在隧道工程中引入了反分析方法。反分析方法是構(gòu)建關(guān)鍵參數(shù)（邊界參數(shù)、彈性模量）的一個(gè)函數(shù)（線性或非線性）來(lái)表示整個(gè)區(qū)域內(nèi)的位移或應(yīng)力分布，調(diào)整方程的系數(shù)直到與參考值相符[18-19]。該方法只適用于儲(chǔ)層性質(zhì)相對(duì)均勻的薄層，不適用于非均質(zhì)性較強(qiáng)，地質(zhì)條件復(fù)雜的地層。Ito 等[20]提出了一種基于有限元計(jì)算的地質(zhì)模型參數(shù)反演方法，該方法將應(yīng)力場(chǎng)視為上覆巖層應(yīng)力和構(gòu)造應(yīng)變誘導(dǎo)應(yīng)力的疊加。張社榮等[21]通過(guò)在線性回歸方程中考慮不同應(yīng)力分量之間的相互作用，改進(jìn)了該方法。近年來(lái)，一些更復(fù)雜的反演算法被應(yīng)用到了巖石力學(xué)問(wèn)題當(dāng)中。Nguyen 和Nestorovi[22]以隧道引起的沉降和水平位移為參考，采用非線性卡爾曼濾波算法對(duì)隧道地質(zhì)力學(xué)模型中的地質(zhì)材料參數(shù)進(jìn)行校準(zhǔn)。Li 等[23]提出了一種結(jié)合貝葉斯方法與多向量輸出算法（B–MSVM）的方法來(lái)反演地應(yīng)力模型的楊氏模量和側(cè)壓系數(shù)。Zhuang 等[24]采用多策略人工魚(yú)群優(yōu)化向量回歸算法，對(duì)隧道工程中圍巖力學(xué)參數(shù)進(jìn)行反算。然而，這種方法只能提供應(yīng)力分布，且地質(zhì)力學(xué)參數(shù)場(chǎng)沒(méi)有調(diào)整。上述方法在一定條件下對(duì)應(yīng)力場(chǎng)計(jì)算具有良好的適用性，但是并不對(duì)地質(zhì)力學(xué)參數(shù)場(chǎng)進(jìn)行更新。地震數(shù)據(jù)是解釋彈性參數(shù)（楊氏模量和泊松比）的可靠依據(jù)[25]。對(duì)于無(wú)法獲得可靠三維地震資料的情況，通常采用基于一維解釋資料的隨機(jī)算法（克里金法，高斯序貫?zāi)M）生成屬性場(chǎng)。在給定隨機(jī)參數(shù)和硬數(shù)據(jù)（測(cè)井或解釋數(shù)據(jù)）的情況下，可以生成任意數(shù)量的實(shí)現(xiàn)，但由于網(wǎng)格數(shù)量通常較大，人為調(diào)整每個(gè)網(wǎng)格的值非常困難[12]。Gao 等[26]采用逐次線性估計(jì)（SLE）算法對(duì)隧道開(kāi)挖問(wèn)題的單元地質(zhì)力學(xué)參數(shù)場(chǎng)進(jìn)行了統(tǒng)計(jì)校正，但這種方法的有效性受網(wǎng)格的質(zhì)量及數(shù)量等影響較大。為了加速對(duì)隨機(jī)場(chǎng)的同化，Karhunen–Loève 展開(kāi)（KL-E）已被多次用于滲流問(wèn)題的分析[27-29]。KL-E 利用特征函數(shù)和隨機(jī)變量的一系列線性組合產(chǎn)生不確定參數(shù)場(chǎng)，可以將屬性模型的輸入變量從數(shù)千個(gè)基于單元的屬性減少為一組數(shù)量較少的一維隨機(jī)變量，大幅度提高計(jì)算效率。本文將KL-E 算法用于地質(zhì)力學(xué)參數(shù)場(chǎng)的生成，以提高模型參數(shù)場(chǎng)的同化效率。

本文提出一個(gè)序貫地質(zhì)力學(xué)模型參數(shù)同化算法，同時(shí)對(duì)地質(zhì)力學(xué)參數(shù)場(chǎng)和模型邊界條件進(jìn)行調(diào)整。該方法利用Matlab 優(yōu)化工具箱中的全局優(yōu)化算法——遺傳算法（Genetic Algorithm,GA）實(shí)現(xiàn)隨機(jī)變量的同化，同時(shí)利用局部?jī)?yōu)化算法模式搜索（Pattern Search,PS）進(jìn)行邊界條件的反演。首先，采用KL E 生成地質(zhì)力學(xué)參數(shù)隨機(jī)場(chǎng)；然后，通過(guò)兩步序貫方法分別對(duì)邊界條件和參數(shù)場(chǎng)進(jìn)行同化。目前，計(jì)算中只考慮了楊氏模量屬性場(chǎng)和水平應(yīng)變邊界條件。文中給出的測(cè)試算例與四川盆地頁(yè)巖氣的實(shí)例均論證了該方法的有效性。計(jì)算結(jié)果表明，邊界條件的同化相對(duì)容易，而屬性場(chǎng)的收斂是相對(duì)緩慢和困難的。但通過(guò)對(duì)地質(zhì)力學(xué)參數(shù)場(chǎng)進(jìn)行同化計(jì)算，可進(jìn)一步提高模型計(jì)算應(yīng)力與一維應(yīng)力參考值的吻合度，從而提高整個(gè)地質(zhì)力學(xué)模型的可靠性。

1 地質(zhì)背景

本文分別用PS 和GA 對(duì)邊界條件和地質(zhì)力學(xué)參數(shù)場(chǎng)進(jìn)行同化。目前，只考慮楊氏模量的場(chǎng)分布，其分布由KL-E 生成?？紤]到一維解釋?xiě)?yīng)力或?qū)崪y(cè)應(yīng)力位置通常也有測(cè)井資料可用于解釋沿井眼的楊氏模量，本文在最后的實(shí)例計(jì)算中進(jìn)一步采用條件KL-E 算法，即在采用KL-E 算法的同時(shí)，滿足參考點(diǎn)參數(shù)值不變的前提條件，既能生成滿足統(tǒng)計(jì)參數(shù)的隨機(jī)場(chǎng)，又能保持有硬數(shù)據(jù)的位置參數(shù)值不變。目前，只進(jìn)行二維地質(zhì)力學(xué)參數(shù)場(chǎng)的同化。

1.1 地質(zhì)力學(xué)參數(shù)隨機(jī)場(chǎng)的生成

假設(shè)儲(chǔ)層地質(zhì)力學(xué)性質(zhì)的分布滿足高斯分布，就可利用KL-E 描述楊氏模量場(chǎng)。如果隨機(jī)場(chǎng)滿足高斯分布，則參數(shù)場(chǎng)分布可以分解為兩部分

式中：x=(x1,x2)—某點(diǎn)坐標(biāo)，m；

)—隨機(jī)場(chǎng)的平均值，單位與變量單位一致；

Y′(x)—隨機(jī)場(chǎng)的擾動(dòng)或空間波動(dòng)，單位與變量單位一致。

函數(shù)Y(x) 的協(xié)方差矩陣CY中關(guān)于位置x,y的協(xié)方值寫(xiě)為CY(x,y)，該函數(shù)有界、對(duì)稱且正定。CY(x,y)可以擴(kuò)展為

式中：

CY(x,y)—關(guān)于位置x,y 的協(xié)方值，單位為變量單位的平方；

y=(y1,y2)—另一點(diǎn)空間坐標(biāo)，m；

λn—協(xié)方差矩陣的特征值，為正數(shù)，無(wú)因次；

fn(x)，fn(y)—相互正交的協(xié)方差矩陣的特征函數(shù)，與待求變量單位一致。

特征值與特征函數(shù)之間滿足

式中：fi—特征向量。

本文采用的協(xié)方差計(jì)算公式為

式中：

σ2—待求變量的方差，單位為待求變量單位的平方；

η1—平行于x1方向的變量分布相關(guān)長(zhǎng)度，m；

η2—平行于x2方向的變量分布相關(guān)長(zhǎng)度，m。

通過(guò)將協(xié)方差矩陣分解為特征值和特征函數(shù)的組合，可以通過(guò)式（5）生成隨機(jī)場(chǎng)

式中：ξn—滿足方差為1、均值為0 正態(tài)分布的隨機(jī)數(shù)，無(wú)因次。

由于特征值的急劇衰減，只需使用式（5）中部分隨機(jī)變量逼近隨機(jī)場(chǎng)，這將大大節(jié)省計(jì)算時(shí)間。截?cái)囗?xiàng)為N的KL-E 可用于近似參數(shù)場(chǎng)分布

其中，特征值按最大值到最小值排列。截?cái)囗?xiàng)個(gè)數(shù)的確定采用保留特征值之和大于所有特征值之和80%的原則[29]。

為了考慮同化前獲得的硬數(shù)據(jù)，本文采用條件KL-E。條件KL-E 可以在給定位置的硬數(shù)據(jù)不變的情況下產(chǎn)生隨機(jī)場(chǎng)。本文采用了與Xue 等[30]相同的算法。

假設(shè)共有nY個(gè)測(cè)量點(diǎn)（給定硬數(shù)據(jù)所在位置），測(cè)量值分別為Y1(x1),Y2(x2),···,YnY(xnY)；坐標(biāo)為x1,x2,···,xnY。條件克里金平均值和協(xié)方差函數(shù)可以計(jì)算如下

式中：上標(biāo)c—條件；

μi(x)—每個(gè)測(cè)量點(diǎn)的權(quán)重函數(shù)，無(wú)因次，根據(jù)克里金公式計(jì)算為

確定后即可根據(jù)式（8）計(jì)算協(xié)方差矩陣并根據(jù)式（3）計(jì)算特征值與特征函數(shù)。最終條件隨機(jī)場(chǎng)表達(dá)式為

(x)—條件協(xié)方差矩陣特征函數(shù)，與變量單位一致。

為了減少計(jì)算時(shí)間，特征值和特征函數(shù)只在第一次實(shí)現(xiàn)計(jì)算一次，其他現(xiàn)實(shí)可直接采用。由于巖土力學(xué)性質(zhì)應(yīng)為正，因此，采用Y(x) >0 約束楊氏模量參數(shù)場(chǎng)取值。利用條件KL-E 進(jìn)行參數(shù)場(chǎng)同化，只需調(diào)整隨機(jī)參數(shù)ξn的值，而非每個(gè)網(wǎng)格的參數(shù)值，從而大大節(jié)省了同化參數(shù)的數(shù)量，提高了模型同化和收斂的效率。

1.2 地質(zhì)力學(xué)有限元模型

利用有限元模擬建立地質(zhì)力學(xué)模型并進(jìn)行計(jì)算。本文測(cè)試算例忽略了孔隙彈性效應(yīng)，假定孔隙壓力為常數(shù)且等于0，而在現(xiàn)場(chǎng)實(shí)例研究中將考慮孔隙彈性效應(yīng)和孔隙壓力的非均質(zhì)性。

為避免有限元模擬時(shí)可能出現(xiàn)的邊界效應(yīng)，先對(duì)目標(biāo)儲(chǔ)層網(wǎng)格進(jìn)行擴(kuò)展[31]。模型由上層（綠色）、側(cè)部（白色）、下層（紅色）和目標(biāo)層（深藍(lán)色）等4 個(gè)區(qū)域組成（圖1）。每個(gè)區(qū)域的網(wǎng)格如圖2 所示，油藏網(wǎng)格中的隨機(jī)顏色代表不同的物質(zhì)類型。由于ABAQUS 中使用的顏色是有限的，因此，某些具有不同屬性的網(wǎng)格將以相同的顏色顯示。本文只考慮目標(biāo)儲(chǔ)層的應(yīng)力和地質(zhì)力學(xué)性質(zhì)分布，其余區(qū)域用于消除模型邊界附近的應(yīng)力集中。

圖1 不同區(qū)域示意圖Fig.1 Schematic of different regions

圖2 不同區(qū)域網(wǎng)格Fig.2 The grid of different regions

有限元模型邊界條件施加：水平方向施加法向應(yīng)變邊界條件（εxx、εyy—x,y 方向的應(yīng)變量，即圖3中x,y 方向的邊界條件）。由于所建模型沒(méi)有擴(kuò)展到地表以減少網(wǎng)格數(shù)，因此，在模型頂部施加額外的面力（σzz，圖3 中的垂向邊界條件）以表征上層巖石的重力作用，同時(shí)固定另外3 個(gè)表面的法向位移（圖3）。

圖3 有限元模擬中使用的邊界條件Fig.3 The boundary conditions used in the FEM simulation

在同化過(guò)程中，上層、側(cè)部和下層網(wǎng)格的地質(zhì)力學(xué)屬性參數(shù)保持不變，僅目標(biāo)儲(chǔ)層中每個(gè)網(wǎng)格的性質(zhì)以及邊界條件隨著計(jì)算過(guò)程而改變，目標(biāo)層中網(wǎng)格屬性根據(jù)1.1 中介紹的條件KL–E 算法獲得。本文利用內(nèi)部編寫(xiě)的Python 腳本實(shí)現(xiàn)ABAQUS 有限元模型的網(wǎng)格屬性自動(dòng)賦值和后臺(tái)運(yùn)行。

1.3 序貫兩步同化算法

本文提出一個(gè)序貫兩步同化算法：首先，采用模式搜索（PS）算法對(duì)邊界條件（xx和yy）進(jìn)行同化；隨后，利用遺傳算法（GA）實(shí)現(xiàn)巖土力學(xué)參數(shù)（楊氏模量）的同化。若需要進(jìn)一步提高模型與參考值的匹配程度，可進(jìn)行另一輪反演，直到滿足收斂準(zhǔn)則為止。詳細(xì)工作流程如圖4 所示。

圖4 序貫兩步同化過(guò)程示意圖Fig.4 Schematic of the sequentially two-step assimilation procedure

優(yōu)化的目標(biāo)是找到邊界條件或地質(zhì)參數(shù)的值，以便將觀測(cè)位置處的模擬應(yīng)力和解釋?xiě)?yīng)力之間的差降到最小。目標(biāo)函數(shù)F(x)定義為

式中：

σmn,ir—通過(guò)一維地質(zhì)力學(xué)分析解釋得到的地應(yīng)力值，MPa；

σmn,i—三維地質(zhì)力學(xué)模型計(jì)算的應(yīng)力，MPa；

M—參考點(diǎn)數(shù)量，個(gè)。

水平應(yīng)變邊界條件范圍可通過(guò)式（12）估算[32]

式中：

σh—最小主應(yīng)力，MPa；

v—泊松比，無(wú)因次；

σv—垂向應(yīng)力，MPa；

α—比奧系數(shù)，無(wú)因次；

pp—孔隙壓力，MPa；

E—楊氏模量，GPa；

εh—最小水平應(yīng)變，無(wú)因次；

εH—最大水平應(yīng)變，無(wú)因次；

σH—最大主應(yīng)力，MPa。

隨機(jī)場(chǎng)生成需要的均值、方差與相關(guān)長(zhǎng)度等參數(shù)從參考點(diǎn)的解釋數(shù)據(jù)中獲得，KL-E 算法中隨機(jī)變量取值在[-2,2]。

2 結(jié)果與討論

2.1 測(cè)試算例

針對(duì)測(cè)試算例，應(yīng)用上述序貫兩步同化算法反演楊氏模量參數(shù)場(chǎng)和水平應(yīng)變邊界條件。測(cè)試算例中，目標(biāo)儲(chǔ)層網(wǎng)格數(shù)為20×20×1，網(wǎng)格尺寸為5 m×5 m×10 m。上、下層網(wǎng)格厚度均為100 m，網(wǎng)格厚度從目標(biāo)儲(chǔ)層厚度以1.5 倍的比例遞增。整個(gè)模型網(wǎng)格數(shù)為9 216（圖2）。固定4 個(gè)網(wǎng)格的楊氏模量參數(shù)：（5,50），（50,5），（50,100），（100,50），假設(shè)為參考測(cè)量點(diǎn)。采用條件KL-E 算法生成楊氏模量的隨機(jī)場(chǎng)。部分模型關(guān)鍵參數(shù)如表1 所示，邊界條件如表2 所示。這些參數(shù)從實(shí)際油藏模型中提取，在測(cè)試模型中，使用了更顯著的方差與更短的相關(guān)長(zhǎng)度以加強(qiáng)地質(zhì)力學(xué)屬性場(chǎng)的非均質(zhì)性。假設(shè)楊氏模量屬性場(chǎng)沿x與y 方向的相關(guān)長(zhǎng)度均為100 m。地質(zhì)力學(xué)參數(shù)隨機(jī)場(chǎng)的協(xié)方差矩陣中，前20 個(gè)特征值的總和超過(guò)了所有特征值和的85%，因此，在后續(xù)的隨機(jī)場(chǎng)生成過(guò)程中將只對(duì)前20 個(gè)特征向量對(duì)應(yīng)的隨機(jī)數(shù)進(jìn)行同化調(diào)整。設(shè)置遺傳算法每一代人口數(shù)為20，收斂條件為連續(xù)兩代的平均目標(biāo)函數(shù)值差異小于5%。

表1 測(cè)試算例屬性參數(shù)Tab.1 Property values for the synthetic case

表2 測(cè)試算例模型邊界條件Tab.2 Boundary conditions for the synthetic case

該測(cè)試算例將進(jìn)行兩次同化計(jì)算，即進(jìn)行兩輪邊界應(yīng)變和楊氏模量場(chǎng)的序貫同化。測(cè)試算例預(yù)設(shè)參數(shù)及各步驟的結(jié)果分別記錄在表3 與表4 中。楊氏模量參考場(chǎng)、初始楊氏模量和每一步的楊氏模量場(chǎng)如圖5 所示，對(duì)應(yīng)的應(yīng)力場(chǎng)如圖6～圖9 所示。本文定義張拉應(yīng)力為正。

圖5 楊氏模量場(chǎng)同化過(guò)程Fig.5 The assimilation process of Young′s modulus field

圖6 參考場(chǎng)下地應(yīng)力分布Fig.6 Distribution of in-situ stresses of reference field

表4 測(cè)試算例計(jì)算結(jié)果Tab.4 Results of synthetic case

根據(jù)數(shù)值計(jì)算結(jié)果，經(jīng)過(guò)25 代遺傳算法的同化，目標(biāo)值的變化基本停止。因此，在進(jìn)行場(chǎng)同化時(shí)設(shè)定遺傳算法在經(jīng)過(guò)25 代計(jì)算后停止。在進(jìn)行邊界條件的同化時(shí)，本文設(shè)定當(dāng)兩次迭代的相對(duì)差值小于1%時(shí)，PS 算法停止。

表3 表明，利用序貫兩步同化程序，可以有效地同化地質(zhì)力學(xué)模型的邊界條件和目標(biāo)值。此外，邊界條件的擬合速度快、精度高，表明邊界條件與地應(yīng)力具有很強(qiáng)的相關(guān)性。同時(shí)，首次擬合得到的目標(biāo)函數(shù)值同化程度最高，第二次擬合對(duì)目標(biāo)函數(shù)值的改變相對(duì)較小。

表3 測(cè)試算例預(yù)設(shè)參數(shù)Tab.3 Initial guess for the synthetic case

圖5 為楊氏模量場(chǎng)隨同化過(guò)程進(jìn)行的演化過(guò)程。由于已知4 參考點(diǎn)楊氏模量值，初始實(shí)現(xiàn)能夠部分獲取空間楊氏模量分布趨勢(shì)，但仍與參考場(chǎng)有較大差別，隨著同化過(guò)程的進(jìn)行，楊氏模量場(chǎng)分布與參考場(chǎng)不斷接近，即使在無(wú)參考點(diǎn)的區(qū)域仍然能夠捕捉到大部分場(chǎng)分布信息。圖6 為利用參考楊氏模量場(chǎng)與邊界條件計(jì)算得到的x方向與y 方向總應(yīng)力分布，從圖中可以看出，楊氏模量較大的區(qū)域?qū)?yīng)壓應(yīng)力較大，應(yīng)力分布趨勢(shì)受楊氏模量分布影響明顯。圖7 為初始楊氏模量場(chǎng)與邊界條件下的地應(yīng)力分布，可以看到，其分布與參考場(chǎng)有較大分別。圖8與圖9 分別為第一輪與第二輪同化后應(yīng)力分布，可以看到，當(dāng)完成第一輪同化后，應(yīng)力分布與參考場(chǎng)已較為接近。

圖7 初始場(chǎng)地應(yīng)力分布Fig.7 Distribution of in-situ stresses of the initial field

圖8 第一輪同化后地應(yīng)力分布Fig.8 Distribution of in-situ stresses after the first two-step assimilation

圖9 第二輪同化后地應(yīng)力分布Fig.9 Distribution of in-situ stresses after the second two-step assimilation

由圖5～圖9 可以看出，同化楊氏模量場(chǎng)和應(yīng)力場(chǎng)都能捕捉參考場(chǎng)的趨勢(shì)和主要特征。與隨機(jī)生成的初始場(chǎng)相比，同化后的模型具有與參考場(chǎng)更為相似的特征，在遺傳算法中，通過(guò)多輪迭代或增加每一代的種群數(shù)量，可以進(jìn)一步提高算法的精度，但同時(shí)也會(huì)增加計(jì)算時(shí)間。

2.2 實(shí)例研究

本文實(shí)例研究對(duì)象為四川盆地龍馬溪組某頁(yè)巖氣井區(qū)。首先，利用測(cè)井、室內(nèi)實(shí)驗(yàn)、水力壓裂以及鉆井等資料，對(duì)4 口水平壓裂井（圖10）進(jìn)行了一維地質(zhì)力學(xué)分析，獲得了沿井軌跡的地質(zhì)力學(xué)參數(shù)（楊氏模量、泊松比以及巖石硬度）和三向主應(yīng)力（Sxx,Syy,Szz）。

圖10 實(shí)際算例網(wǎng)格Fig.10 Grids of the practical case

三維地質(zhì)力學(xué)模型與測(cè)試算例相似，地質(zhì)力學(xué)模型包括目標(biāo)儲(chǔ)層、側(cè)部、上部和下部地層。模型尺度為2 800 m×2 800 m×10 m，單元尺度為dx=100 m，dy=100 m，dz=10 m。網(wǎng)格頂部位于地表以下1 000 m處，網(wǎng)格底部位于海拔-5 000 m處，由于本文僅考慮二維隨機(jī)場(chǎng)的同化，因此，選擇生產(chǎn)潛力較大的層作為目標(biāo)層，原始地質(zhì)模型中主力產(chǎn)層包含5 層細(xì)網(wǎng)格，這里只對(duì)其中物性最好（滲透率與孔隙度較高）的一層網(wǎng)格層進(jìn)行屬性場(chǎng)反演。目標(biāo)層共有8 個(gè)網(wǎng)格與井徑交叉，將沿井軌跡與交叉網(wǎng)格相交部分的一維解釋值進(jìn)行平均獲得網(wǎng)格平均值（圖11a），在這8 個(gè)網(wǎng)格處固定地質(zhì)力學(xué)參數(shù)，且認(rèn)定為應(yīng)力參考點(diǎn)，參考點(diǎn)網(wǎng)格編號(hào)如圖11b 所示。

圖11 測(cè)井?dāng)?shù)據(jù)在網(wǎng)格中的分布Fig.11 The distribution of logging data in grids

表5 列出了每個(gè)區(qū)域的材料特性，其中，楊氏模量與泊松比分布在x與y 方向的相關(guān)長(zhǎng)度均為1 000 m。需要通過(guò)匹配觀測(cè)網(wǎng)格中的垂直應(yīng)力來(lái)調(diào)整上層巖石的密度（圖12），從圖中可以看出，從一維分析得到的解釋垂直應(yīng)力和計(jì)算垂直應(yīng)力較為接近，各點(diǎn)編號(hào)在圖11 中標(biāo)出，擬合后上層巖石密度為2.55 kg/m3。由于地層實(shí)際情況與地質(zhì)力學(xué)模型存在差異，解釋和計(jì)算的垂直應(yīng)力之間仍存在不同，但相對(duì)差異較小。

圖12 垂向應(yīng)力解釋值與計(jì)算值對(duì)比Fig.12 Comparison of veritcal stress in interperted and calculated results

利用8 個(gè)參考點(diǎn)的一維解釋值可以計(jì)算楊氏模量和泊松比的均值、方差和相關(guān)長(zhǎng)度，實(shí)際算例不同區(qū)域參數(shù)如表5 所示。

表5 各區(qū)域材料特性Tab.5 Material properties in each region

研究區(qū)為高壓氣藏，壓力系數(shù)大于1.8。根據(jù)區(qū)塊多口直井的試壓資料，觀察到壓力系數(shù)與深度之間存在較好的相關(guān)性，因此，在本算例中直接根據(jù)實(shí)際的地層深度來(lái)計(jì)算孔隙壓力（圖13），計(jì)算公式如式（13）所示。

圖13 目的層孔隙壓力分布Fig.13 The distribution of pore pressure in target reservoir layer

式中：h—垂深，m；g—重力加速度，g=9.8 m/s2。

根據(jù)水平井的水力壓裂和成像測(cè)井以及微地震資料得出，最小水平應(yīng)力近似平行于y軸，因此，假定主應(yīng)力方向沿x，y 和z軸。目標(biāo)函數(shù)，即解釋?xiě)?yīng)力和計(jì)算應(yīng)力之間的差為

本算例中，參考點(diǎn)數(shù)量為8 個(gè)。利用式（12）估計(jì)邊界條件的初始猜測(cè)，并用條件KL-E 生成在固定參考點(diǎn)屬性值的楊氏模量場(chǎng)。表6 為實(shí)際算例的預(yù)設(shè)參數(shù)，表7 為同化過(guò)程中的計(jì)算結(jié)果及優(yōu)化算法選取參數(shù)。由表7 可以看出，第一輪兩步同化可以有效降低目標(biāo)值，而第二輪同化后目標(biāo)值變化不大，該結(jié)論與測(cè)試算例結(jié)論一致。與測(cè)試算例不同的是，實(shí)際算例優(yōu)化后的最終目標(biāo)值仍然較大。該現(xiàn)象可以從以下幾個(gè)方面解釋。首先，解釋的一維地質(zhì)力學(xué)性質(zhì)和應(yīng)力可能與實(shí)際值不同。此外，密度分布、泊松比和地質(zhì)模型的幾何形狀也可能對(duì)算法的準(zhǔn)確性產(chǎn)生影響。同時(shí)，從8 個(gè)點(diǎn)得到的隨機(jī)參數(shù)可能不足以反映真實(shí)場(chǎng)的特征。由于觀測(cè)到目標(biāo)值顯著降低，計(jì)算結(jié)果仍然證明了該方法在實(shí)際算例中的有效性。

表6 實(shí)際算例預(yù)設(shè)參數(shù)Tab.6 Initial guess of practical case

表7 實(shí)際算例計(jì)算結(jié)果Tab.7 Calculating results of practical case

初始場(chǎng)和同化場(chǎng)的應(yīng)力比較如圖14～圖16 所示。由于第二輪同化產(chǎn)生的目標(biāo)值變化較小，這里只展示第一輪同化后的結(jié)果。圖17 將數(shù)據(jù)同化過(guò)程中每個(gè)步驟后的計(jì)算應(yīng)力與解釋值進(jìn)行比較，圖中參考為一維地質(zhì)力學(xué)分析的解釋?xiě)?yīng)力，初始為通過(guò)邊界條件和楊氏模量分布初始猜測(cè)計(jì)算的應(yīng)力，同化BC 為第一輪邊界條件同化后獲得的結(jié)果，同化E 為第一輪楊氏模量同化后的結(jié)果。從圖17 可以看出，解釋?xiě)?yīng)力和初始猜測(cè)之間的差異相當(dāng)顯著。僅對(duì)邊界條件進(jìn)行同化不足以使模型結(jié)果與參考值匹配。楊氏模量場(chǎng)的同化有助于進(jìn)一步提高地質(zhì)力學(xué)模型與實(shí)際情況的一致性。

圖14 楊氏模量場(chǎng)同化結(jié)果Fig.14 Assimilation results of Young′s modulus field

圖15 初始地應(yīng)力場(chǎng)分布Fig.15 The initial distribution of in-situ stresses

圖16 一輪同化后地應(yīng)力方向Fig.16 The distribution of in-situ stresses after 1st round of assimilation

圖17 參考點(diǎn)沿x和y 方向的地應(yīng)力比較Fig.17 Comparison of in-situ stresses alongxand y directions at observation locations

3 結(jié)論

（1）利用KL-E 生成地質(zhì)力學(xué)參數(shù)場(chǎng)可以大大減少每次同化計(jì)算時(shí)需要更新的未知數(shù)。測(cè)試算例計(jì)算結(jié)果表明，采用本文提出的算法同化得到的地質(zhì)力學(xué)屬性參數(shù)場(chǎng)分布與實(shí)際情況較為接近。當(dāng)應(yīng)用于實(shí)際算例時(shí)，本文算法也可有效提高模型計(jì)算應(yīng)力與參考應(yīng)力的吻合度。

（2）相比屬性參數(shù)場(chǎng)，邊界條件的同化收斂速度更快。通常只需要兩次序貫同化邊界條件和地質(zhì)力學(xué)屬性參數(shù)場(chǎng)，目標(biāo)值將停止顯著變化。