王依誠，姜漢橋，于馥瑋，成寶洋，徐飛，李俊鍵

中國石油大學(xué)(北京)油氣探測(cè)與工程國家重點(diǎn)實(shí)驗(yàn)室，北京 102249

0 引言

巖心數(shù)字化是智能油田發(fā)展的重要一環(huán)，國內(nèi)各大油田都在進(jìn)行建設(shè)自動(dòng)化的立體巖心庫，巖心數(shù)字化后，研究人員不受時(shí)間、地域限制，隨時(shí)對(duì)巖心進(jìn)行分析，提高了生產(chǎn)決策效率[1-2]?？紫抖燃皾B透率是研究流體巖石滲流規(guī)律的主要表征參數(shù)，也是儲(chǔ)集層產(chǎn)能挖潛和提高采收率的關(guān)鍵[3-6]，在巖心數(shù)字化建設(shè)過程中具有重要意義。

巖心孔滲參數(shù)獲取主要有兩個(gè)手段：(1)通過巖心驅(qū)替實(shí)驗(yàn)測(cè)定，此類方法需要特定的實(shí)驗(yàn)條件和手段，且對(duì)巖心有一定的損傷； (2)通過數(shù)字巖心建模并進(jìn)行數(shù)值模擬，目前主要以X射線CT掃描為結(jié)果，建立真實(shí)三維多孔介質(zhì)模型，并通過格子玻爾茲曼方法或者有限元模擬方法進(jìn)行模擬求解。2008年朱益華等[7]從現(xiàn)場(chǎng)采集的彩色鑄體剖面圖中提取孔隙信息重構(gòu)了一種接近真實(shí)巖石的3D孔隙介質(zhì)數(shù)字巖心并用格子玻爾茲曼方法研究孔隙巖石的滲流特性。2015年宋睿等[8]基于CT圖像通過有限元方法建立了原始巖樣孔隙形態(tài)的結(jié)構(gòu)化網(wǎng)絡(luò)模型，并求解了巖樣滲透率數(shù)據(jù)。數(shù)值模擬方法受到計(jì)算資源的限制，其模擬尺寸一般小于500×500×500網(wǎng)格，且計(jì)算非常耗時(shí)。2013年Amabeoku 等[9]采用并行運(yùn)算機(jī)制，利用有限元方法將碳酸鹽巖數(shù)字巖心的模擬尺寸做到2000×2000×2000網(wǎng)格，將模擬運(yùn)算時(shí)間從兩周縮減為一天。

目前，缺少一種對(duì)于巖心孔隙度及滲透率的快速預(yù)測(cè)方法，這嚴(yán)重阻礙了數(shù)字化巖心庫的建設(shè)過程，機(jī)器學(xué)習(xí)算法是一類從數(shù)據(jù)中自動(dòng)分析獲得規(guī)律，并利用規(guī)律對(duì)未知數(shù)據(jù)進(jìn)行預(yù)測(cè)的算法，已廣泛應(yīng)用于數(shù)據(jù)挖掘、計(jì)算機(jī)視覺、自然語言處理、生物特征識(shí)別等領(lǐng)域。而巖心孔隙度及滲透率參數(shù)本質(zhì)上是表征巖心多孔介質(zhì)的流動(dòng)規(guī)律的特征，因此在問題類型上與機(jī)器學(xué)習(xí)算法的應(yīng)用領(lǐng)域一致。近年已有學(xué)者將機(jī)器學(xué)習(xí)應(yīng)用到流體力學(xué)領(lǐng)域，2002年Milano等[10]通過神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)湍流過程中近壁速度場(chǎng)；2019年Rahman等[11]基于LSTM(長(zhǎng)短期記憶模型)提高湍流流動(dòng)模擬中的數(shù)值模擬質(zhì)量。

本文主要通過建立實(shí)際巖心CT掃描結(jié)果建立數(shù)字巖心孔滲特征對(duì)應(yīng)數(shù)據(jù)庫，通過機(jī)器學(xué)習(xí)算法對(duì)數(shù)據(jù)庫進(jìn)行訓(xùn)練，實(shí)現(xiàn)數(shù)字巖心孔滲參數(shù)的快速預(yù)測(cè)，為油田數(shù)字化建設(shè)提供有效的指導(dǎo)作用。

1 基于機(jī)器學(xué)習(xí)的數(shù)字巖心孔滲預(yù)測(cè)方法

基于機(jī)器學(xué)習(xí)的數(shù)字巖心的孔滲預(yù)測(cè)方法以歷史數(shù)據(jù)為基礎(chǔ)，通過實(shí)驗(yàn)測(cè)定方法、測(cè)井結(jié)果、三維孔隙網(wǎng)絡(luò)數(shù)值模擬結(jié)果等為數(shù)字巖心添加上數(shù)據(jù)標(biāo)簽，形成孔隙預(yù)測(cè)學(xué)習(xí)數(shù)據(jù)庫，通過模型對(duì)數(shù)據(jù)庫的學(xué)習(xí)，當(dāng)再有新的數(shù)據(jù)時(shí)，只需要通過CT掃描，生成原始CT掃描數(shù)據(jù)，即可快速完成巖心的孔滲參數(shù)預(yù)測(cè)，無需再進(jìn)行實(shí)驗(yàn)測(cè)定或者數(shù)值模擬，能夠極大提高參數(shù)獲取效率，具體工作流程如圖1所示。

1.1 卷積神經(jīng)網(wǎng)絡(luò)

巖心通過X射線CT掃描獲取的數(shù)據(jù)為巖心二維灰度圖像序列，通過選取合適的尺寸將二維切片圖像 序列通過重建算法重構(gòu)得到最終的三維數(shù)字巖心體，即為機(jī)器學(xué)習(xí)的輸入數(shù)據(jù)，而輸出數(shù)據(jù)則為巖心的孔隙度及滲透率參數(shù)，為一組連續(xù)數(shù)值結(jié)果。根據(jù)輸入輸出數(shù)據(jù)特性，本文主要通過三維卷積神經(jīng)網(wǎng)絡(luò)構(gòu)建數(shù)字巖心孔滲預(yù)測(cè)模型。卷積神經(jīng)網(wǎng)絡(luò)(Convolutional Neural Network, CNN)是一種前饋神經(jīng)網(wǎng)絡(luò)，它的人工神經(jīng)元可以響應(yīng)一部分覆蓋范圍內(nèi)的周圍單元，對(duì)于大型圖像處理有出色表現(xiàn)。1998年LeCun等[12-14]人確立了CNN的現(xiàn)代結(jié)構(gòu)，CNN能夠得出原始圖像的有效表征，這使得CNN能夠直接從原始像素中，經(jīng)過極少的預(yù)處理，識(shí)別視覺上的規(guī)律。典型的卷積神經(jīng)網(wǎng)絡(luò)包括輸入層、卷積層、池化層、全連接層和輸出層(如圖 2)。

1.1.1 卷積層

卷積層主要作用是識(shí)別輸入數(shù)據(jù)的特征，一個(gè)卷積核可以提取一個(gè)特征，不同卷積核可以提取圖像的不同特征。CNN的最大特點(diǎn)是局部感知和權(quán)值共享[15-18]。一般認(rèn)為人對(duì)外界的認(rèn)知是從局部到全局的，而圖像的空間聯(lián)系也是局部的像素聯(lián)系較為緊密，而距離較遠(yuǎn)的像素相關(guān)性則較弱。因而，每個(gè)神經(jīng)元其實(shí)沒有必要對(duì)全局圖像進(jìn)行感知，只需要對(duì)局部進(jìn)行感知，然后在更高層將局部的信息綜合起來就得到了全局的信息。為了檢測(cè)圖像不同位置的同一特征，某個(gè)局部區(qū)域的權(quán)值可以共享用于另一個(gè)區(qū)域(即權(quán)值共享性)。局部感知和權(quán)值共享保證模型準(zhǔn)確提取特征的同時(shí)，減少權(quán)值參數(shù)的個(gè)數(shù)，從而大幅度減少計(jì)算復(fù)雜度。卷積層的計(jì)算公式為：

其中，σ表示激活函數(shù)，I表示輸入圖像灰度矩陣，ω表示卷積核，°表示卷積操作，b表示偏置項(xiàng)。

卷積層提取特征的過程如圖3所示，原始輸入像素為10×10,選擇5×5卷積核，進(jìn)行卷積運(yùn)算，最終得到(10-5+1) ×(10-5+1)=6×6特征圖。

圖1 基于機(jī)器學(xué)習(xí)的數(shù)字巖心孔滲預(yù)測(cè)流程Fig. 1 The work flow of digital core pore permeability prediction based on machine learning

圖2 卷積神經(jīng)網(wǎng)絡(luò)模型的基本結(jié)構(gòu)Fig. 2 Architecture for convolutional neural network

1.1.2 池化層

池化層主要是在卷積后來簡(jiǎn)化卷積層輸出維度來降低計(jì)算復(fù)雜度，是一種統(tǒng)計(jì)聚合操作，常見的池化方法有最大值池化和平均池化。以最大池化為例，圖4左圖所展示的池化前特征圖，選擇池化過濾器為2×2，步長(zhǎng)為2，其左邊部分中左上角2×2的矩陣中6最大，右上角2×2的矩陣中8最大，左下角2×2的矩陣中3最大，右下角2×2的矩陣中4最大，所以其池化結(jié)果如圖4右所示。

圖3 二維卷積操作過程Fig. 3 Architecture for convolutional neural network

1.1.3 全連接層

區(qū)別于上述局部連接的概念, 全連接表示每個(gè)神經(jīng)元都與相鄰層的所有神經(jīng)元互相連接。全連接層實(shí)質(zhì)上就是一個(gè)多層感知器，根據(jù)研究目的可選擇不同的感知器。如需要進(jìn)行分類，則輸出層使用Softmax激勵(lì)函數(shù)；如果任務(wù)是回歸，則輸出層不需要使用Softmax函數(shù)，直接對(duì)前一層進(jìn)行累加即可。對(duì)于整個(gè)CNN模型來說, 卷積層和池化層稱為特征提取器, 全連接層就是目標(biāo)轉(zhuǎn)換函數(shù)。

1.2 三維卷積神經(jīng)網(wǎng)絡(luò)

二維卷積在圖像處理有著不凡的效果，然而在處理視頻、CT數(shù)據(jù)上無法反應(yīng)時(shí)間、空間上的特征信息[19]，需要通過三維卷積神經(jīng)網(wǎng)絡(luò)(3D CNN)來完成此類任務(wù)。三維卷積與二維卷積主要在卷積、池化時(shí)多了一個(gè)維度，具體區(qū)別見圖5。

圖4 池化過程示意圖Fig. 4 Architecture for Convolutional Neural Network

圖5 三維卷積示意圖Fig. 5 Architecture for Convolutional Neural Network

2 實(shí)例應(yīng)用

2.1 學(xué)習(xí)數(shù)據(jù)庫建立

機(jī)器學(xué)習(xí)訓(xùn)練數(shù)據(jù)庫是本文建立網(wǎng)絡(luò)的重要組成部分，為了驗(yàn)證模型的預(yù)測(cè)能力，本文對(duì)具有不同滲透率的天然巖心的CT掃描結(jié)果進(jìn)行提取分割，并采用OpenFOAM開源框架建立三維孔隙網(wǎng)絡(luò)模型，計(jì)算每個(gè)樣本的孔滲參數(shù)，形成訓(xùn)練數(shù)據(jù)庫。

實(shí)驗(yàn)巖心選擇江蘇油田3塊砂巖油藏天然巖心，基本參數(shù)見表1。對(duì)3塊巖心進(jìn)行Micro-CT掃描，CT掃描分辨率為2.98 μm，掃描后尺寸原始尺寸為1920×1920×1536像素，考慮到巖心掃描邊緣效應(yīng)，每塊巖心掃描結(jié)果將截取1200×1200×1200像素作為基礎(chǔ)數(shù)據(jù)，圖像二值化算法主要采用孔隙度校正方法保證圖像分割的正確性，二值化結(jié)果見圖6。

表1 實(shí)驗(yàn)巖心基本參數(shù)Table 1 Parameters of experimental core

考慮到計(jì)算機(jī)運(yùn)算能力，對(duì)每個(gè)基礎(chǔ)數(shù)據(jù)進(jìn)行下向采樣完成分割數(shù)據(jù)，每組基礎(chǔ)數(shù)據(jù)的采樣數(shù)為：

圖6 巖心二維切片F(xiàn)ig. 6 Two-dimensional cross-sectional images of cores

式中，Iori為原始數(shù)據(jù)像素尺寸,Isub為采樣子塊像素尺寸，Istride為采樣間隔像素，考慮到計(jì)算資源，Isub取200，Istride取0，每塊巖心原始數(shù)據(jù)能夠獲取218個(gè)子塊數(shù)據(jù)3塊巖心共生成654個(gè)子塊數(shù)據(jù)，訓(xùn)練數(shù)據(jù)流程見圖7。

巖石的絕對(duì)滲透率是巖石的固有屬性，其大小僅僅取決于巖石的孔隙結(jié)構(gòu)，是用來表征油藏的非均質(zhì)性和各向異性的重要參數(shù)，利用達(dá)西定律(公式3)可計(jì)算巖石單相絕對(duì)滲透率。

圖7 訓(xùn)練數(shù)據(jù)生成流程圖Fig. 7 Work flow of training data generation

式中，v表示達(dá)西速度；k為絕對(duì)滲透率；ΔP為驅(qū)替壓差；μ為流體黏度；L為巖樣的長(zhǎng)度。

當(dāng)考慮到所通過巖石的流體為不可壓縮流體，密度為定常數(shù)；流體為牛頓流體，其動(dòng)態(tài)黏度為常數(shù)；流動(dòng)過程為穩(wěn)定流，流體流速不會(huì)隨時(shí)間變化；流速較小，流動(dòng)過程為層流?？梢詫avier-Stakes方程簡(jiǎn)化為：

式中，??表示散度算子，?2表示梯度算子，v為介質(zhì)中流相的流體速度，?為拉普拉斯算符，p表示介質(zhì)中流相流體的壓力。

對(duì)于每個(gè)子塊，其數(shù)據(jù)維度為200×200×200的二值化三維數(shù)據(jù)體，通過開源軟件OpenFoam建立三維單相流模型，通過求解Navier-Stakes方程，能夠?qū)Ψ€(wěn)定狀態(tài)的壓力及流速進(jìn)行求取，模型壓差設(shè)置為1 Pa，模擬時(shí)間總時(shí)間為1 s，時(shí)間步長(zhǎng)為0.1 s，計(jì)算結(jié)果如圖8、圖9所示。

圖8 三維單相流模型壓力分布(1 s)Fig. 8 Pressure distribution of three-dimensional singlephase flow model (1 s)

圖9 三維單相流模型流速分布(1 s)Fig. 9 Volocity distribution of three-dimensional singlephase flow model (1 s)

根據(jù)模型求得的穩(wěn)定狀態(tài)下的截面流速，通過式(3)獲取到每個(gè)數(shù)字巖心子塊的滲透率，654組樣本計(jì)算結(jié)果如圖10所示。樣本有效孔隙度分布范圍為0.05～0.35，平均孔隙度為0.21，滲透率分布范圍為5～4042 mD，平均滲透率為616 mD，有效孔隙度與滲透率在對(duì)數(shù)坐標(biāo)上呈明顯的直線關(guān)系，因此在設(shè)計(jì)訓(xùn)練樣本時(shí)，選擇有效孔隙度和滲透率的對(duì)數(shù)值作為預(yù)測(cè)結(jié)果。

2.2 卷積神經(jīng)網(wǎng)絡(luò)設(shè)計(jì)

模型每個(gè)樣本的輸入數(shù)據(jù)為三維數(shù)字巖心數(shù)據(jù)體，維度為200×200×200，輸出維度為2×1，分別代表孔隙度的大小和滲透率的對(duì)數(shù)，根據(jù)訓(xùn)練數(shù)據(jù)輸入輸出維度，3D CNN模型的網(wǎng)絡(luò)結(jié)構(gòu)設(shè)計(jì)參數(shù)如表2所示。

模型主要由卷積層、池化層、激活層和全連接層構(gòu)成。模型卷積核設(shè)置均為5×5×5，激活函數(shù)均采用ReLU激活函數(shù)，誤差函數(shù)選擇平方和誤差函數(shù)。每層對(duì)應(yīng)的特征圖個(gè)數(shù)分別為4，4，8，8，16，16，32，32，64，64，13 824，2。對(duì)于全連接層，運(yùn)用一個(gè)13 824維的全連接層將分布式特征數(shù)據(jù)集集合到一個(gè)13 824維的樣本空間中，全連接層進(jìn)行特征融合，最后的二維結(jié)果分別代表預(yù)測(cè)的有效孔隙度和滲透率對(duì)數(shù)值。將654組三維CT數(shù)據(jù)分隨機(jī)挑選50組作為測(cè)試數(shù)據(jù)，其余604組數(shù)據(jù)作為訓(xùn)練數(shù)據(jù)，迭代次數(shù)設(shè)置為訓(xùn)練代數(shù)設(shè)置為400。

圖10 訓(xùn)練樣本孔隙度滲透率分布Fig. 10 Porosity and permeability distribution of training data

表2 3D CNN模型參數(shù)設(shè)置Table 2 The parameter of 3D CNN

2.3 損失函數(shù)

損失函數(shù)是模型對(duì)數(shù)據(jù)擬合程度的反映，由于本次研究目的是預(yù)測(cè)，所以選擇最小平方誤差判別(Minimum Squared-Error，MSE)作為損失函數(shù)。

式中，yi和f(xi)分別表示第i個(gè)樣本真實(shí)值和預(yù)測(cè)值，m為樣本個(gè)數(shù)。訓(xùn)練的目標(biāo)是使MSE變小，因此訓(xùn)練過程就成為了一個(gè)無約束優(yōu)化過程。

2.4 學(xué)習(xí)率

學(xué)習(xí)率是卷積神經(jīng)網(wǎng)絡(luò)中的一個(gè)重要的超參數(shù)，合適的學(xué)習(xí)率是訓(xùn)練出好模型的關(guān)鍵要素之一，本次模型學(xué)習(xí)率分別設(shè)置了(0.001,0.003,0.005,0.01,0.03)以研究不同學(xué)習(xí)率對(duì)模型訓(xùn)練結(jié)果的影響。

2.5 結(jié)果分析

圖11展示了不同學(xué)習(xí)率下相同參數(shù)3D CNN模型訓(xùn)練過程中訓(xùn)練誤差的變化情況，當(dāng)學(xué)習(xí)率α=0.03時(shí)，由于學(xué)習(xí)率過大，訓(xùn)練誤差振蕩嚴(yán)重，訓(xùn)練一定代數(shù)誤差曲線開始上揚(yáng)；α=0.001時(shí)，由于學(xué)習(xí)率過小，模型收斂速度慢，在達(dá)到迭代次數(shù)后模型仍未達(dá)到較好的訓(xùn)練結(jié)果； 當(dāng)α=0.005和0.01時(shí)，模型快速收斂，在較小的迭代次數(shù)后即達(dá)到水平，后期訓(xùn)練誤差基本保持不變；α=0.003時(shí)，模型訓(xùn)練誤差下降較為平穩(wěn)，在訓(xùn)練后期訓(xùn)練逐漸平穩(wěn)且收斂速度較高，因此本文選定α=0.003訓(xùn)練3D CNN預(yù)測(cè)模型。

圖11 不同學(xué)習(xí)率3D CNN模型訓(xùn)練誤差變化Fig. 11 Training error of 3D CNN model with different learning rates

圖12 測(cè)試數(shù)據(jù)累計(jì)誤差分布Fig. 12 accumulative error of test data

圖13 測(cè)試數(shù)據(jù)孔隙度滲透率分布Fig. 13 Porosity and permeability distribution of test data

通過測(cè)試集50組數(shù)據(jù)對(duì)模型的泛化能力進(jìn)行檢測(cè)，從圖12中可以看出，模型孔隙度預(yù)測(cè)誤差小于5%的占95.74%，滲透率預(yù)測(cè)誤差小于10%的占91.49%，并且在孔滲關(guān)系上與真實(shí)數(shù)據(jù)有著較強(qiáng)的符合關(guān)系(如圖13)，說明基于此模型具有較強(qiáng)的泛化能力。

對(duì)于建立的預(yù)測(cè)模型，從訓(xùn)練數(shù)據(jù)中選取10個(gè)樣本與OpenFOAM計(jì)算速度進(jìn)行對(duì)比。OpenFOAM每次均需要進(jìn)行建模及數(shù)模運(yùn)算，數(shù)值模擬方法平均用時(shí)在為1.58 h，而通過訓(xùn)練好的機(jī)器學(xué)習(xí)模型預(yù)測(cè)用時(shí)僅為0.03 s，為數(shù)字巖心滲透率的快速預(yù)測(cè)提供了有效的解決方案。

3 結(jié)論

(1)根據(jù)實(shí)際天然巖心CT掃描結(jié)果運(yùn)用開源流體力學(xué)計(jì)算模塊OpenFOAM生成了數(shù)字巖心孔滲預(yù)測(cè)數(shù)據(jù)庫，通過三維卷積神經(jīng)網(wǎng)絡(luò)建立了數(shù)字巖心孔隙快速預(yù)測(cè)模型，并取得了較好的訓(xùn)練結(jié)果，說明機(jī)器學(xué)習(xí)在數(shù)字巖心特征提取上具有較大潛力。

(2)分析了學(xué)習(xí)率對(duì)模型訓(xùn)練效果的影響，學(xué)習(xí)率過大訓(xùn)練誤差會(huì)發(fā)生大幅振蕩且后期曲線上揚(yáng)，學(xué)習(xí)率過小則收斂速度較慢，本次模型的最佳學(xué)習(xí)率為0.003，能夠達(dá)到較快的收斂速度且能有效避免過擬合。

(3)分析預(yù)測(cè)結(jié)果與真實(shí)孔滲，預(yù)測(cè)誤差小于10%的概率結(jié)果占90%以上，說明模型具有較強(qiáng)的泛化能力，能夠?qū)?shù)字巖心進(jìn)行快速孔滲預(yù)測(cè)，建立的預(yù)測(cè)模型與傳統(tǒng)數(shù)值模擬方法在預(yù)測(cè)速度上具有數(shù)量級(jí)的提升，實(shí)現(xiàn)了數(shù)字巖心孔滲高效率、高精度預(yù)測(cè)，能夠有效降低生產(chǎn)成本，大幅度提高工作效率。

4 展望

目前，前人通過機(jī)器學(xué)習(xí)在圖像、視頻及語音領(lǐng)域進(jìn)行了大量的研究，但對(duì)基于CT掃描結(jié)果建模的數(shù)字巖心孔滲參數(shù)預(yù)測(cè)還尚未有報(bào)道。本文提出了采用三維卷積神經(jīng)網(wǎng)絡(luò)進(jìn)行的數(shù)字巖心孔滲參數(shù)預(yù)測(cè)方法，該方法通過三維卷積能夠?qū)?shù)字巖心進(jìn)行自動(dòng)特征提取，并完成快速預(yù)測(cè)，訓(xùn)練好的模型預(yù)測(cè)耗時(shí)小于1 s，展示了三維卷積神經(jīng)網(wǎng)絡(luò)對(duì)于數(shù)字巖心特征提取的潛力。

然而，機(jī)器學(xué)習(xí)對(duì)于數(shù)字巖心更為復(fù)雜的特征提取是否具有同樣的效果還需要進(jìn)一步驗(yàn)證，需要在以下方面進(jìn)行深入研究：

(1)數(shù)字巖心兩相滲流特征。目前采用的訓(xùn)練樣本是采用單相流模擬結(jié)果進(jìn)行預(yù)測(cè)，可進(jìn)一步選擇兩相流模擬生成訓(xùn)練樣本，探索機(jī)器學(xué)習(xí)對(duì)數(shù)字巖心相滲的預(yù)測(cè)。

(2)對(duì)于孔喉的非均質(zhì)性，可以通過對(duì)水驅(qū)后的數(shù)字巖心油水飽和度結(jié)果進(jìn)行學(xué)習(xí)，預(yù)測(cè)數(shù)字巖心不同的孔喉動(dòng)用的難易程度，為儲(chǔ)層評(píng)價(jià)提供新思路。