王志強(qiáng)，李航宇，劉樹陽，徐建春，范晨

(中國(guó)石油大學(xué)(華東) 石油工程學(xué)院，山東 青島，266580)

在全球氣候變化日益嚴(yán)峻的背景下，溫室氣體的排放與利用正成為各國(guó)家與機(jī)構(gòu)重點(diǎn)關(guān)注的話題[1]。2020年，全球與能源相關(guān)的二氧化碳排放量為315億t[2]。為了應(yīng)對(duì)碳排放問題，我國(guó)明確提出CO2排放力爭(zhēng)于2030 年前達(dá)到峰值，努力爭(zhēng)取2060年前實(shí)現(xiàn)碳中和的“雙碳目標(biāo)”[3]。

在石油化工領(lǐng)域，二氧化碳捕集、利用與封存(CCUS)被認(rèn)為是安全、有效地實(shí)現(xiàn)碳“負(fù)排放”的路徑之一[4]。CCUS 是指將CO2從工業(yè)排放源中分離后或直接加以利用或地質(zhì)封存，以實(shí)現(xiàn)CO2減排的工業(yè)過程[5]。政府間氣候變化專門委員會(huì)表明：在2 ℃時(shí)，不采用CCUS，溫室氣體減排成本將增加138%。國(guó)際能源署表明CCUS 在各類減排技術(shù)中的貢獻(xiàn)度將達(dá)到15%[6]。

2021 年，全球的商業(yè)CCUS 設(shè)施大量增加，達(dá)到135 個(gè)，年捕集能力達(dá)到1.493 億t[7]。在全球商業(yè)化CCUS 項(xiàng)目中，將CO2直接注入咸水層中，實(shí)現(xiàn)永久性封存被認(rèn)為是極具潛力的方式之一[8]。中國(guó)深部咸水層的CO2封存能力約為31 350 億t，其中陸地咸水層和海域咸水層埋存占比分別為74%和26%。油藏和氣藏的封存能力分別為48.24億t和51.94億t[9]。

目前工業(yè)上采取的CO2注入方式是依靠廢棄油氣井或新建注入井，將液態(tài)或超臨界態(tài)CO2通過井筒注入目標(biāo)層位[10]。CO2注入能力是實(shí)施咸水層碳封存項(xiàng)目的重要參考依據(jù)，也是影響注入方案設(shè)計(jì)、注入設(shè)備配套等問題的關(guān)鍵參數(shù)。從經(jīng)濟(jì)性出發(fā)，CO2地質(zhì)封存項(xiàng)目往往需要以最大的注入速率進(jìn)行，因此，CO2注入井井底可能達(dá)到所允許的最大注入壓力。此外，CO2的注入能力還會(huì)受到地質(zhì)參數(shù)[11-13]和操作參數(shù)[14-15]的綜合影響。然而，目前針對(duì)礦場(chǎng)尺度CO2注入能力的評(píng)估大多采用傳統(tǒng)的數(shù)值模擬方法，計(jì)算時(shí)間較長(zhǎng)，模擬時(shí)考慮的影響因素也比較單一，沒有建立復(fù)雜多變量條件下的CO2注入能力的預(yù)測(cè)體系。

近年來，隨著人工智能技術(shù)的應(yīng)用與推廣，油田生產(chǎn)也進(jìn)入信息化、智能化時(shí)代。BP(back propagation)神經(jīng)網(wǎng)絡(luò)[16]是一種按誤差反向傳播訓(xùn)練的多層前饋網(wǎng)絡(luò)，具備很強(qiáng)的非線性映射能力，網(wǎng)絡(luò)的中間層數(shù)和各層的神經(jīng)元個(gè)數(shù)可以根據(jù)情況任意調(diào)節(jié)，具備結(jié)構(gòu)簡(jiǎn)單、操作性強(qiáng)、收斂速度快等優(yōu)點(diǎn)，被廣泛應(yīng)用于非線性建模、函數(shù)逼近、模式識(shí)別等方面。在油氣領(lǐng)域，BP 神經(jīng)網(wǎng)絡(luò)也大量應(yīng)用于產(chǎn)能預(yù)測(cè)、錄井解釋和注采方案優(yōu)化等方面[17-19]。張延旭等[20]采用BP 神經(jīng)網(wǎng)絡(luò)法預(yù)測(cè)了油藏封存CO2的效果，埋存系數(shù)的預(yù)測(cè)結(jié)果與實(shí)際數(shù)值模擬結(jié)果相對(duì)誤差在8%以內(nèi)，但是沒有對(duì)激活函數(shù)、隱藏層神經(jīng)元數(shù)量及學(xué)習(xí)率進(jìn)行優(yōu)化，所得模型的泛化能力較差；VO-THANH等[21]使用廣義回歸神經(jīng)網(wǎng)絡(luò)模型(GRNN)預(yù)測(cè)咸水層CO2封存過程中的殘余捕集與溶解捕集指數(shù)，取得了較好的效果，決定系數(shù)分別達(dá)到0.999 5 和0.999 8，但是該模型的輸入?yún)?shù)忽略了模型尺寸對(duì)注入能力的影響。

綜上所述，目前人工智能算法在咸水層CO2封存中注入能力評(píng)估的應(yīng)用較少，尚未建立考慮多因素、復(fù)雜條件的注入能力預(yù)測(cè)模型。為了更加準(zhǔn)確、快速地預(yù)測(cè)咸水層CO2封存的注入能力，本文首先考慮不同地質(zhì)參數(shù)及操作參數(shù)的復(fù)雜條件，并進(jìn)行數(shù)值模擬；其次，基于BP基神經(jīng)網(wǎng)絡(luò)的方法形成模擬-預(yù)測(cè)耦合技術(shù)，對(duì)激活函數(shù)、隱藏層神經(jīng)元數(shù)量及學(xué)習(xí)率進(jìn)行優(yōu)化，構(gòu)建了可以替代傳統(tǒng)數(shù)值模擬的非線性回歸模型；最后，根據(jù)目標(biāo)封存區(qū)塊的地質(zhì)屬性與操作參數(shù)，實(shí)現(xiàn)CO2注入能力的快速精準(zhǔn)預(yù)測(cè)，以便為咸水層CO2封存方案的制定提供技術(shù)支撐。

1 基于BP神經(jīng)網(wǎng)絡(luò)的模擬-預(yù)測(cè)耦合原理

1.1 BP神經(jīng)網(wǎng)絡(luò)基本原理

圖1所示為BP神經(jīng)網(wǎng)絡(luò)示意圖。由圖1所見：BP 神經(jīng)網(wǎng)絡(luò)模型拓?fù)浣Y(jié)構(gòu)包括輸入層、隱藏層和輸出層。算法流程包括信號(hào)的正向傳播和誤差的反向傳播2個(gè)過程。正向傳播時(shí)，輸入層的數(shù)值通過激活函數(shù)加權(quán)求和計(jì)算，傳播到隱藏層，再以相同的方式傳播到輸出層。若在輸出層得到的預(yù)測(cè)值不滿足精度要求，則將誤差信號(hào)沿原來的路徑返回，通過梯度下降法實(shí)現(xiàn)權(quán)重更新，然后繼續(xù)下一步的正向傳播，通過該方式不斷降低誤差，最終達(dá)到所需精度。

圖1 BP神經(jīng)網(wǎng)絡(luò)示意圖Fig.1 Schematic diagram of BP neural network

對(duì)于具有單一輸出變量的訓(xùn)練集D={(xd，y)}m，訓(xùn)練可以轉(zhuǎn)化為下面的優(yōu)化問題：

式中：d為輸入數(shù)據(jù)維度；m為訓(xùn)練樣本容量；ν為輸入層到隱藏層的連接權(quán)；ω為隱藏層到輸出層的連接權(quán)；γ為隱藏層閾值；θ為輸出層閾值；yk為第k個(gè)數(shù)據(jù)的真實(shí)輸出值；y^k為第k個(gè)數(shù)據(jù)的預(yù)測(cè)輸出值。

1.2 模擬-預(yù)測(cè)耦合技術(shù)流程

圖2 所示為模擬-預(yù)測(cè)耦合技術(shù)的基本流程。模擬-預(yù)測(cè)耦合技術(shù)的基本流程如下：

圖2 模擬-預(yù)測(cè)耦合技術(shù)的基本流程Fig.2 Basic flow of simulation-prediction coupling technology

1) 使用CMG-GEM 商業(yè)數(shù)值模擬軟件，建立考慮咸水層CO2封存機(jī)理的數(shù)值模型，獲取特定參數(shù)下的CO2注入能力。

2) 改變地質(zhì)模型參數(shù)及操作參數(shù)，得到不同特征下的CO2注入能力，構(gòu)建用于機(jī)器學(xué)習(xí)的樣本空間。

3) 建立BP神經(jīng)網(wǎng)絡(luò)模型，通過參數(shù)優(yōu)化生成CO2注入能力的代理模型，進(jìn)而可以預(yù)測(cè)任意地質(zhì)模型參數(shù)及操作參數(shù)下的注入能力。

2 學(xué)習(xí)樣本庫的生成

2.1 數(shù)值模擬模型

2.1.1 數(shù)值模擬軟件選擇及驗(yàn)證

CO2在咸水層中封存的主要機(jī)理[22]為：構(gòu)造封存，束縛空間封存(或殘余封存)，溶解封存及礦化封存。目前，可有效模擬咸水層CO2封存過程中的各種機(jī)理的軟件包括CMG-GEM、TOUGH2/ECO2N和ECLIPSE等[23]。何斌等[12]使用TOUGH2/ECO2N 軟件研究鄂爾多斯盆地深部咸水層的CO2注入能力。本文首先依據(jù)何斌等[12]提供的參數(shù)建立相同的地質(zhì)模型，然后對(duì)比相同條件下的CO2注入速率，結(jié)果如圖3 所示。由圖3 可見：CMGGEM 可較好地模擬咸水層CO2的注入情況，與TOUGH2/ECO2N 模擬結(jié)果基本一致。另外，CMG-GEM 軟件可以更加方便地實(shí)現(xiàn)敏感性參數(shù)組合與方案設(shè)計(jì)。因此，本文采用CMG-GEM 軟件完成模擬-預(yù)測(cè)耦合技術(shù)中的模擬部分。

圖3 CMG-GEM數(shù)值模擬軟件驗(yàn)證Fig.3 Verification of numerical simulation software CMG-GEM

2.1.2 地質(zhì)模型建立及網(wǎng)格劃分

CO2的目標(biāo)儲(chǔ)層為咸水層，一般由高孔隙度、高滲透率的砂巖層組成。在咸水層上方有一定厚度的非滲透性蓋層，一般為泥頁巖層，依靠其極低的滲透率和較高的毛管力，對(duì)CO2形成較強(qiáng)的阻力。咸水層CO2地質(zhì)封存模型示意圖如圖4 所示。由圖4 可見：在注氣井的對(duì)角線處設(shè)置1 口生產(chǎn)井，將井底壓力設(shè)置為原始地層壓力。當(dāng)生產(chǎn)井附近壓力上升時(shí)，即可通過生產(chǎn)井排出一部分咸水，為CO2的順利注入提供空間。

圖4 咸水層CO2地質(zhì)封存模型示意圖Fig.4 Schematic diagram of CO2 storage in saline aquifer

基礎(chǔ)方案中模型長(zhǎng)×寬×高為1 000 m×1 000 m×400 m，其中蓋層厚度為100 m，儲(chǔ)層厚度為 300 m?；A(chǔ)方案模型的網(wǎng)格長(zhǎng)×寬×高為50 m× 50 m×10 m。因此，x，y和z方向上的網(wǎng)格數(shù)量分別為20，20和40個(gè)，總網(wǎng)格數(shù)為16 000個(gè)?；A(chǔ)方案模型網(wǎng)格劃分如圖5所示?；A(chǔ)方案中頂界面深度為1 200 m，儲(chǔ)層面積為1 km2。

圖5 基礎(chǔ)方案模型網(wǎng)格劃分Fig.5 Meshing of base model

2.1.3 初始條件與邊界條件

基礎(chǔ)模型初始條件如表1所示，模型外邊界條件為封閉邊界。

表1 基礎(chǔ)方案模型初始條件Table 1 Initial conditions of base model

2.2 代理模型樣本庫的建立

沈平平等[24]建立了咸水層CO2注入能力的經(jīng)驗(yàn)公式：

式中：q為井筒中的流速，m3/d；K為儲(chǔ)層的滲透率，10-3μm2；h為儲(chǔ)層厚度，m；Δp為儲(chǔ)層與井筒間的壓差，MPa；μ為注入相的黏度，mPa·s；re為等效泄流半徑，m；rw為井筒半徑，m。

式(2)表明，注入能力受到儲(chǔ)層滲透率、厚度、儲(chǔ)層及井筒間的壓差、注入相的黏度、等效泄流半徑及井筒半徑的影響。何斌等[12]采用井筒-儲(chǔ)層耦合的方法，發(fā)現(xiàn)增加孔隙度可以提高CO2的注入能力，但效果遠(yuǎn)沒有增加滲透率顯著。

因此，為確定不同地質(zhì)參數(shù)和操作參數(shù)對(duì)注入能力的影響，以儲(chǔ)層孔隙度、儲(chǔ)層滲透率、儲(chǔ)層面積、儲(chǔ)層厚度、垂向與水平滲透率比kv/kh、注入壓力與儲(chǔ)層壓力的壓差為自變量，根據(jù)實(shí)際咸水層CO2封存項(xiàng)目及數(shù)值模擬經(jīng)驗(yàn)設(shè)置取值范圍，如表2所示。

表2 代理模型輸入?yún)?shù)取值Table 2 Value of input parameters of proxy model

在設(shè)計(jì)模擬方案時(shí)，各參數(shù)每個(gè)水平的全組合數(shù)為4 320，數(shù)量較大，因此，使用拉丁超立方方法[25]進(jìn)行抽樣。通過該方法得到220組數(shù)值模擬方案，如表3 所示，作為下一步BP 神經(jīng)網(wǎng)絡(luò)訓(xùn)練和預(yù)測(cè)的樣本空間。

通過現(xiàn)有咸水層注入方案和數(shù)值模擬的分析，BP 神經(jīng)網(wǎng)絡(luò)的輸出參數(shù)選擇為保持恒定注入壓力下5 a的累計(jì)注氣量。各模擬方案的運(yùn)行結(jié)果如表3所示。

3 代理模型的求解、驗(yàn)證與應(yīng)用

3.1 數(shù)據(jù)預(yù)處理

由表3 可知，各輸入?yún)?shù)之間存在不同的量綱，參數(shù)的跨度范圍過大，嚴(yán)重影響神經(jīng)網(wǎng)絡(luò)模型的學(xué)習(xí)效果，因此，采用線性歸一化方法對(duì)輸入?yún)?shù)進(jìn)行處理，將所有結(jié)果映射到[0，1]之間。轉(zhuǎn)換函數(shù)如下：

表3 代理模型樣本空間Table 3 Database of proxy model

式中：x′為歸一化參數(shù)；x為實(shí)際值；max(x)和min(x)分別為變量x的最大值和最小值。

3.2 參數(shù)優(yōu)化

為得到適合該樣本空間的最優(yōu)BP 神經(jīng)網(wǎng)絡(luò)，需要對(duì)激活函數(shù)、隱藏層神經(jīng)元數(shù)量及學(xué)習(xí)率進(jìn)行優(yōu)化。神經(jīng)網(wǎng)絡(luò)模型參數(shù)優(yōu)化取值如表4所示。

表4 神經(jīng)網(wǎng)絡(luò)模型參數(shù)優(yōu)化取值Table 4 Value for parameter optimization of neural network model

為避免測(cè)試集與訓(xùn)練集隨機(jī)選取導(dǎo)致數(shù)據(jù)劃分的敏感性，通常采用“十折交叉驗(yàn)證”的方法對(duì)數(shù)據(jù)進(jìn)行分組。“十折交叉驗(yàn)證”方法原理如圖6所示。首先將樣本空間分為10組，輪流取1組作為驗(yàn)證集，其余作為訓(xùn)練集，得到相應(yīng)的損失函數(shù)；然后，將10 輪損失函數(shù)的平均值作為最終的目標(biāo)函數(shù)E。損失函數(shù)分布越集中，目標(biāo)函數(shù)值越小，說明該預(yù)測(cè)效果越好。

圖6 十折交叉驗(yàn)證原理示意圖Fig.6 Schematic diagram of tenfold cross-validation principle

誤差函數(shù)選取平均相對(duì)誤差R1為

式中：n為樣本數(shù)量；yi為實(shí)際樣本標(biāo)簽值；為代理模型反歸一化后的輸出值。

1) 激活函數(shù)。激活函數(shù)首先對(duì)輸入信息進(jìn)行非線性變換，然后將變換后的輸出信息作為輸入信息傳給下一層神經(jīng)元。不同的激活函數(shù)有不同的優(yōu)缺點(diǎn)及適用情形。對(duì)于表4中給定的每一種激活函數(shù)，通過“十折交叉驗(yàn)證”方法均可得到10個(gè)損失函數(shù)的結(jié)果，將其繪制在箱線圖上，如圖7所示。每種激活函數(shù)最終的目標(biāo)函數(shù)值即為箱線圖中的“均值”。由圖7可見：satlin函數(shù)對(duì)應(yīng)的損失函數(shù)相對(duì)誤差分布最集中，而且最終的目標(biāo)函數(shù)值E也比較低，沒有異常值點(diǎn)。

圖7 BP神經(jīng)網(wǎng)絡(luò)激活函數(shù)優(yōu)選Fig.7 Optimization of activation function of BP neural network

2) 隱藏層神經(jīng)元數(shù)量。合理地選擇隱藏層神經(jīng)元數(shù)量，對(duì)模型的性能有很大的影響。隱藏層神經(jīng)元數(shù)量較少會(huì)導(dǎo)致欠擬合，數(shù)量過多又會(huì)導(dǎo)致過擬合。設(shè)置隱藏層神經(jīng)元數(shù)量范圍為1～20，對(duì)于每1個(gè)隱藏層神經(jīng)元數(shù)量，繪制相應(yīng)的損失函數(shù)箱線圖，具體做法與優(yōu)化激活函數(shù)的方法相同。最終繪制的損失函數(shù)箱線圖如圖8所示。由圖8可見：當(dāng)隱藏層神經(jīng)元數(shù)量為10 時(shí)，損失函數(shù)相對(duì)誤差分布比較集中，整體水平較低，沒有異常值點(diǎn)。

圖8 BP神經(jīng)網(wǎng)絡(luò)隱藏層神經(jīng)元數(shù)量?jī)?yōu)選Fig.8 Optimization of number of neurons in hidden layer of BP neural network

3) 學(xué)習(xí)率。合適的學(xué)習(xí)率能夠使目標(biāo)函數(shù)在合適的時(shí)間內(nèi)收斂到局部最小值。如果學(xué)習(xí)率較小，那么消耗大量時(shí)間或因梯度消失而無法收斂；如果學(xué)習(xí)率較高，權(quán)重更新幅度太大，導(dǎo)致參數(shù)在極值點(diǎn)兩端發(fā)散。學(xué)習(xí)率的典型范圍是1×10-5～1。不同學(xué)習(xí)率對(duì)應(yīng)的最終目標(biāo)函數(shù)E的分布如圖9所示。從圖9 可見：當(dāng)學(xué)習(xí)率為0.4 時(shí)，最終目標(biāo)函數(shù)值E最低。

圖9 BP神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)率優(yōu)化Fig.9 Optimization of learning rate of BP neural network

通過上述優(yōu)化，最終選取BP 神經(jīng)網(wǎng)絡(luò)的最優(yōu)激活函數(shù)為satlin，最優(yōu)的隱藏層神經(jīng)元數(shù)量為10個(gè)，最優(yōu)的學(xué)習(xí)率為0.4。

3.3 代理模型的驗(yàn)證

使用上述最優(yōu)的參數(shù)構(gòu)建BP 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)，對(duì)樣本空間重新進(jìn)行訓(xùn)練與測(cè)試。隨機(jī)選取10%的數(shù)據(jù)作為測(cè)試集，剩余90%的數(shù)據(jù)作為訓(xùn)練集。使用Pearson(皮爾遜)相關(guān)系數(shù)R評(píng)價(jià)測(cè)試結(jié)果的精度。R越接近于1，說明二者相關(guān)性越高，代理模型的預(yù)測(cè)精度越高。

Pearson(皮爾遜)相關(guān)系數(shù)R可表示為

式中：N為樣本數(shù)量；X為實(shí)際樣本標(biāo)簽值；Y為BP神經(jīng)網(wǎng)絡(luò)代理模型輸出值。

訓(xùn)練集、測(cè)試集與全部數(shù)據(jù)的真實(shí)值與代理模型輸出結(jié)果對(duì)比如圖10 所示。由圖10 可以看出：訓(xùn)練集、測(cè)試集與全部數(shù)據(jù)的皮爾遜相關(guān)系數(shù)分別為0.974 17，0.953 79和0.968 91。這表明該代理模型對(duì)CO2封存能力的預(yù)測(cè)具有較大的可靠性與準(zhǔn)確性。對(duì)于任意給定的CO2封存目標(biāo)區(qū)，即可通過該代理模型進(jìn)行CO2注入能力的快速精準(zhǔn)預(yù)測(cè)，為施工方案的制定提供技術(shù)支撐。

圖10 真實(shí)值與代理模型輸出結(jié)果對(duì)比Fig.10 Comparison between true value and prediction value

3.4 代理模型的應(yīng)用

選取挪威Sleipner、加拿大Quest 及美國(guó)Illinois咸水層CO2封存項(xiàng)目作為測(cè)試，項(xiàng)目的基本參數(shù)如表5所示。

使用上述優(yōu)化后的BP 神經(jīng)網(wǎng)絡(luò)，將表5 中的各輸入?yún)?shù)進(jìn)行歸一化處理后，代入神經(jīng)網(wǎng)絡(luò)模型進(jìn)行預(yù)測(cè)，然后再將得到的輸出值反歸一化，得到代理模型的CO2注入能力預(yù)測(cè)值。從表5可以看出：對(duì)于Sleipner、Quest 和Illinois 項(xiàng)目，該代理模型預(yù)測(cè)的相對(duì)誤差分別為6.85%，5.57%和3.26%，充分說明了該模型的可靠性與準(zhǔn)確性。另外，采用CMG-GEM 軟件模擬基礎(chǔ)模型共耗時(shí)568.7 s，而使用構(gòu)建好的代理模型進(jìn)行訓(xùn)練、測(cè)試及預(yù)測(cè)共耗時(shí)9.5 s，計(jì)算速率得到大幅度提升。

表5 實(shí)際碳封存項(xiàng)目測(cè)試Table 5 Test of actual CO2 storage projects

4 結(jié)論

1) 用于CO2注入能力預(yù)測(cè)的BP 神經(jīng)網(wǎng)絡(luò)的最優(yōu)激活函數(shù)為satlin，最優(yōu)的隱藏層神經(jīng)元數(shù)量為10，最優(yōu)的學(xué)習(xí)率為0.4。

2) 使用最優(yōu)的BP 神經(jīng)網(wǎng)絡(luò)模型參數(shù)，訓(xùn)練集、測(cè)試集與全部數(shù)據(jù)的皮爾遜相關(guān)系數(shù)均超過0.95，驗(yàn)證了模型的準(zhǔn)確性。

3) 使用訓(xùn)練好的模型對(duì)Sleipner，Quest 和Illinois 碳封存項(xiàng)目進(jìn)行測(cè)試，結(jié)果顯示預(yù)測(cè)年注氣量的相對(duì)誤差分別為6.85%，5.57%和3.26%，充分說明了該模型的可靠性。相比較傳統(tǒng)的數(shù)值模擬方法，該代理模型可大幅度提升計(jì)算速率。