孟胤全，蔣建國，吳吉春

南京大學(xué) 地球科學(xué)與工程學(xué)院，表生地球化學(xué)教育部重點(diǎn)實(shí)驗(yàn)室，南京 210023

確定多孔介質(zhì)的滲透率對于地下水運(yùn)動、油氣開采、二氧化碳地下封存以及核廢料泄漏等問題的研究具有重要意義（Neuman, 2005；Tsang et al.,2015；Zhou et al., 2019）。達(dá)西壓力梯度法（Miguel and Serrenho, 2007）測定多孔介質(zhì)滲透率的周期較長且受環(huán)境干擾較大，因此無法大量應(yīng)用。隨著計(jì)算流體力學(xué)的快速發(fā)展，孔隙尺度的多孔介質(zhì)流體流動數(shù)值模擬已經(jīng)廣泛應(yīng)用于計(jì)算滲透率。然而，復(fù)雜的邊界條件與非線性的Navier-Stokes 方程使得數(shù)值模擬計(jì)算成本較高，從而限制了多孔介質(zhì)的計(jì)算尺寸。根據(jù)達(dá)西定律，多孔介質(zhì)的滲透率由孔隙空間結(jié)構(gòu)決定，兩者之間映射關(guān)系可轉(zhuǎn)化為機(jī)器學(xué)習(xí)有監(jiān)督問題。前人選擇使用隨機(jī)森林模型、梯度提升模型、人工神經(jīng)網(wǎng)絡(luò)（ANN）和卷積神經(jīng)網(wǎng)絡(luò)（CNNs）等等機(jī)器學(xué)習(xí)回歸模型預(yù)測多孔介質(zhì)的滲透率（Araya and Ghezzehei, 2019 ；Rabbani et al.,2020 ；Srisutthiyakorn et al., 2016）。

作為計(jì)算機(jī)視覺領(lǐng)域最為經(jīng)典的模型, CNNs已被廣泛用于預(yù)測多孔介質(zhì)數(shù)字圖像所對應(yīng)的滲透率（Kamrava et al., 2019；Tang et al., 2022；Tian et al., 2020；Wu et al., 2018）。然而，孔隙空間結(jié)構(gòu)數(shù)據(jù)量、CNNs 模型的復(fù)雜度與GPU 內(nèi)存之間的矛盾限制了多孔介質(zhì)的研究尺寸與模型的預(yù)測性能（Kashefi et al., 2021）。例如，尺寸為256×256×256立方體素的三維多孔介質(zhì)數(shù)字圖像具有超過1000萬個體素，直接輸入至復(fù)雜度高的深層3D CNN 時會給GPU 帶來巨大的負(fù)擔(dān)。然而降低CNNs 模型的復(fù)雜度會導(dǎo)致欠擬合從而降低預(yù)測性能，增加GPU 的內(nèi)存或數(shù)量則會大幅增加計(jì)算成本。因此前人考慮應(yīng)用多孔介質(zhì)領(lǐng)域的專家知識對孔隙空間結(jié)構(gòu)進(jìn)行手動特征提取，而后再輸入至機(jī)器學(xué)習(xí)模型中。Kashefi 等人將孔隙—固相邊界的位置信息提取為點(diǎn)云數(shù)據(jù)并輸入至點(diǎn)云神經(jīng)網(wǎng)絡(luò)，大幅度降低了模型的數(shù)據(jù)輸入量并取得了很好的預(yù)測結(jié)果。此外，平均曲率積分、多點(diǎn)相關(guān)函數(shù)等孔隙空間結(jié)構(gòu)特征也往往作為人工神經(jīng)網(wǎng)絡(luò)（ANN）或支持向量機(jī)回歸（SVR）或XgBoost 等機(jī)器學(xué)習(xí)模型的輸入（Rabbani et al., 2019 ；R?ding et al., 2020 ；Tian et al., 2020）。然而提取多孔介質(zhì)孔隙結(jié)構(gòu)特征參數(shù)的過程往往復(fù)雜而較難實(shí)現(xiàn)，且參數(shù)一旦冗余并且經(jīng)過組合后輸入至機(jī)器學(xué)習(xí)模型中，則該模型應(yīng)用于滲透率值范圍窄且樣本量少的多孔介質(zhì)數(shù)據(jù)集時，會產(chǎn)生明顯的過擬合（Tian et al., 2020）。因此，本研究考慮一種易于理解與實(shí)現(xiàn)的孔隙結(jié)構(gòu)特征提取方法，能夠通過機(jī)器學(xué)習(xí)模型有效地建立起孔隙結(jié)構(gòu)特征參數(shù)與滲透率的映射關(guān)系。

多孔介質(zhì)的切片在空間上是連續(xù)的，可視為一種空間序列（Zhang et al., 2022）。對于任意一個多孔介質(zhì)，本研究提取其所有切片的單一孔隙結(jié)構(gòu)特征以形成一個多維向量。該向量中的元素數(shù)即多孔介質(zhì)的切片數(shù)。切片的孔隙結(jié)構(gòu)特征保留了空間的連續(xù)性，所組成的向量因此可視為序列數(shù)據(jù)。我們將其作為經(jīng)典機(jī)器學(xué)習(xí)模型與長短期記憶神經(jīng)網(wǎng)絡(luò)（LSTM）模型的輸入。LSTM 模型在自然語言處理（NLP）、語音翻譯等序列建模問題上取得了許多成功（Hochreiter and Schmidhuber, 1997），本文采用該模型來處理切片孔隙結(jié)構(gòu)特征的連續(xù)性。

本文首先介紹多孔介質(zhì)的生成和獲取方法、切片的孔隙結(jié)構(gòu)特征提取方法。其次介紹 3 種常用的經(jīng)典機(jī)器學(xué)習(xí)模型與長短期記憶神經(jīng)網(wǎng)絡(luò)模型的原理。再次，介紹模型的應(yīng)用方法并以兩個案例展示模型在測試集上的預(yù)測性能。而后討論本研究方法的優(yōu)勢。最后，對文章內(nèi)容進(jìn)行總結(jié)。

1 數(shù)據(jù)與方法

1.1 數(shù)據(jù)準(zhǔn)備

應(yīng)用有監(jiān)督機(jī)器學(xué)習(xí)模型預(yù)測多孔介質(zhì)滲透率的前提是獲取足夠多的樣本，包括三維的多孔介質(zhì)數(shù)據(jù)與對應(yīng)的滲透率值。在實(shí)踐過程中采用成像技術(shù)獲取大批量真實(shí)孔隙結(jié)構(gòu)斷層掃描數(shù)據(jù)受時間、經(jīng)濟(jì)成本的制約（Song et al., 2019）。而為了提高效率，利用數(shù)值模擬來生成多孔介質(zhì)這一方法被廣泛采納（Graczyk et al., 2020；Volkhonskiy et al.,2022）。本文采用沉降法生成由球狀顆粒填充的多孔介質(zhì)（Vold, 1960；Pilotti, 1998），如圖1 所示。直徑可變的球顆粒從頂部依次釋放，在下落過程中如果其中一個顆粒撞擊到另一個顆粒，則該顆粒會沿著下部顆粒的表面下落。這一下落過程一直持續(xù)到顆粒勢能達(dá)到局部最小值，此時該顆粒與其他三個球顆粒接觸或撞擊到底部。此沉降法生成的多孔介質(zhì)保留了真實(shí)孔隙結(jié)構(gòu)的拓?fù)涮卣鳎紫抖韧ǔ４笥?.4。為降低孔隙度，我們將部分球顆粒的直徑擴(kuò)大10%，使得孔隙度降低至0.38 以下，同時配位數(shù)保持不變。我們采用周期性邊界條件，將球顆粒沿x和y軸周期性放置。

圖1 生成多孔介質(zhì)的示意圖Fig. 1 The schematic diagram of the generation of porous media

多孔介質(zhì)的滲透率k由達(dá)西定律計(jì)算得出，即u=k ?P/L，其中u為平均流速，?P是壓力差，L是多孔介質(zhì)在流動方向上的長度。因此，當(dāng)計(jì)算出水流的平均流速u時即可確定滲透率值k。Navier-Stokes 方程描述具有復(fù)雜幾何邊界的多孔介質(zhì)中的水流，而格子玻爾茲曼方法（LBM）采用簡便的反彈方案處理無滑移邊界條件，是求解Navier-Stokes方程的有效手段。本文采用LBM 的D3Q19 模型首先計(jì)算出多孔介質(zhì)中z軸方向水流的平均流速，然后再根據(jù)達(dá)西定律計(jì)算出多孔介質(zhì)在該方向上的滲透率（Qian et al., 1992；Wolf-Gladrow, 2004）。

應(yīng)用上述方法本研究生成了1000 個由球顆粒填充的多孔介質(zhì)（圖2a）并計(jì)算出了相應(yīng)的滲透率值，我們將其作為研究案例1。此外，我們還從公開可用的DeePore 多孔介質(zhì)數(shù)據(jù)集（Rabbani et al., 2020）中隨機(jī)選取了5000 個多孔介質(zhì)樣本（圖2b）, 其滲透率通過孔隙網(wǎng)絡(luò)模型（PNM）計(jì)算得出，我們將其作為研究案例2。

圖2 本研究使用的兩類多孔介質(zhì)Fig. 2 Two types of porous media used in this study

1.2 特征提取方法

本研究認(rèn)為，三維多孔介質(zhì)由空間中連續(xù)的二維切片構(gòu)成，從而形成一種空間序列。切片的孔隙結(jié)構(gòu)特征在空間上仍然具有連續(xù)性，因此可視為一組序列數(shù)據(jù)。本研究將三維多孔介質(zhì)表示為由二維孔隙結(jié)構(gòu)特征參數(shù)所組成的多維向量，由于向量的元素具有連續(xù)性，我們稱此類向量為特征序列。特征序列的長度即向量的維數(shù)、切片的數(shù)量。

長期以來，含有孔隙度、比表面積這兩個全局物理參數(shù)的半經(jīng)驗(yàn)公式Kozeny-Carman 方程（Carman，1939）被廣泛應(yīng)用于預(yù)測多孔介質(zhì)的滲透率。本研究提取多孔介質(zhì)切片的孔隙度與比表面積，分別組成孔隙度序列與比表面積序列來表征三維多孔介質(zhì)（圖3）。切片的孔隙度和比表面積可以表示為數(shù)字圖像像素的函數(shù)?？紫抖仁强紫断袼財?shù)與總像素數(shù)之比。二維孔隙結(jié)構(gòu)的比表面積定義為孔隙總周長與其總面積之比（Rabbani et al., 2014；Yu et al., 2009），本研究將比表面積表示為孔隙—顆粒邊界的像素數(shù)與孔隙像素數(shù)之比。此外，本研究還考慮了切片孔隙結(jié)構(gòu)的連通性并以圖像歐拉數(shù)來表示（Ohser et al., 2002 ；Shamsi et al., 2021 ；Vogel et al., 2010）。作為圖像拓?fù)湫再|(zhì)的量度，歐拉數(shù)由孔隙部分的連通分支與其中固相顆粒之差確定。對切片圖像進(jìn)行旋轉(zhuǎn)、拉伸等操作并不會改變歐拉數(shù)，且歐拉數(shù)越小，孔隙的連通性越好。我們借助scikit-image 包，以4-鄰域（1-連接）的模式提取每一切片的圖像歐拉數(shù)。

圖3 提取多孔介質(zhì)切片特征的示意圖Fig. 3 Schematic diagram of extracting slice features of a porous medium

本小節(jié)以4 個樣本為例，說明特征序列具有表征多孔介質(zhì)孔隙空間結(jié)構(gòu)的能力。我們采用沉降法生成兩個孔隙度為0.38 的多孔介質(zhì)作為樣本1 與樣本2，兩者分別由直徑為0.3 mm、0.6 mm的球顆粒填充，z軸方向上的尺寸為400。我們從DeePore 多孔介質(zhì)數(shù)據(jù)集中選取兩個孔隙度為0.10、z軸方向上的尺寸為256 的巖石樣本作為樣本3 和樣本4。我們將這些樣本沿著z軸方向表示為切片序列（圖4），然后由特征提取方法得到每個樣本的孔隙度序列、比表面積序列和歐拉數(shù)序列（圖5）。對于固相均為球顆粒且孔隙度相同的樣本1 與樣本2，兩者孔徑不同這一信息隱含在各自特征序列的震蕩中。對于非均質(zhì)性更強(qiáng)的樣本3 與樣本4，即便孔隙度相同，由不同地質(zhì)條件所形成的孔隙空間結(jié)構(gòu)以及由此產(chǎn)生的特征序列都能夠?qū)烧呙黠@區(qū)分開。因此，我們通過特征提取方法建立的序列能夠表示任意的多孔介質(zhì)并且不會存在重復(fù)。我們計(jì)算每個樣本的三種特征序列兩兩之間的相關(guān)度，并以Pearson 相關(guān)系數(shù)表示，結(jié)果如圖6 所示。我們發(fā)現(xiàn)，對于樣本1 與樣本2 而言，兩者的孔隙度序列與比表面積序列的Pearson 相關(guān)系數(shù)都接近-1，說明由球顆粒填充的多孔介質(zhì)其孔隙度序列與比表面積序列線性相關(guān)性高，且為負(fù)相關(guān)。而孔隙度序列與歐拉數(shù)序列雖然也為負(fù)相關(guān)，但線性相關(guān)性明顯減弱。兩個樣本的比表面積序列與歐拉數(shù)序列的Pearson 相關(guān)系數(shù)分別為0.77 與0.66。對于樣本3 與樣本4 而言，孔隙度序列與比表面積序列之間仍然有很高的負(fù)線性相關(guān)性，但是歐拉數(shù)序列與孔隙度序列或比表面積序列的線性相關(guān)性則顯著減弱。在本研究中，代表多孔介質(zhì)孔隙空間結(jié)構(gòu)的特征序列可以作為經(jīng)典機(jī)器學(xué)習(xí)模型或長短期記憶神經(jīng)網(wǎng)絡(luò)模型的輸入。三種特征序列也可以合并為一個多維數(shù)組而后輸入至長短期記憶神經(jīng)網(wǎng)絡(luò)中。

圖4 沿z軸方向不同位置的多孔介質(zhì)切片F(xiàn)ig. 4 Slices of the porous media at different positions along the z-axis direction

圖5 多孔介質(zhì)的特征序列Fig. 5 Feature sequences of the porous media

圖6 特征序列的Pearson 相關(guān)系數(shù)熱圖Fig. 6 Pearson correlation coefficient heat map of feature sequences

1.3 機(jī)器學(xué)習(xí)方法

本研究采用經(jīng)典機(jī)器學(xué)習(xí)模型k近鄰（KNN）、隨機(jī)森林（Random Forest）、支持向量機(jī)回歸（SVR），長短期記憶神經(jīng)網(wǎng)絡(luò)（LSTM）來預(yù)測多孔介質(zhì)的滲透率。算法介紹如下。

k近鄰算法會尋找最接近輸入樣本的其他k個訓(xùn)練樣本并以它們的平均標(biāo)簽值作為模型的預(yù)測結(jié)果（Cover and Hart，1967）。k近鄰模型的復(fù)雜度由最近鄰數(shù)量k控制，k值越小則表示模型越復(fù)雜。隨機(jī)森林是一種相對簡單的bagging 集成算法（Breiman, 2001），旨在訓(xùn)練多個隨機(jī)且獨(dú)立的決策樹從而提供更好的預(yù)測結(jié)果。支持向量機(jī)回歸算法使得所有樣本數(shù)據(jù)接近超平面并保持總偏差最?。–ortes and Vapnik, 1995），模型采用結(jié)構(gòu)風(fēng)險最小化原則來學(xué)習(xí)數(shù)據(jù)特征，可以很好地處理小樣本集的非線性或高維回歸問題。

長短期記憶神經(jīng)網(wǎng)絡(luò)LSTM 是循環(huán)神經(jīng)網(wǎng)絡(luò)（RNN）的變體（Hochreiter et al., 2001），通過門控機(jī)制來控制網(wǎng)絡(luò)中信息的流動和記憶，可以有效解決反向傳播過程中的梯度消失或梯度爆炸，往往適用于時間序列建模問題。LSTM 的網(wǎng)絡(luò)結(jié)構(gòu)和工作流程如圖7 所示，網(wǎng)絡(luò)中引入長期狀態(tài)ct來存儲特定時刻的歷史信息，并將該信息非線性地輸出到隱藏層的外部狀態(tài)ht。長期狀態(tài)ct的計(jì)算公式如下：

圖7 長短期記憶神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)示意圖Fig. 7 Schematic diagram of the structure of long short-term memory neural network

其中ft∈[0,1]D、it∈[0,1]D與ot∈[0,1]D分別是控制信息傳輸路徑的遺忘門、輸入門與輸出門。⊙表示向量對應(yīng)元素的乘積，ct-1是上一時刻的長期狀態(tài)，c~t是通過雙曲正切函數(shù)得到的候選狀態(tài)，其定義為：

遺忘門、輸入門與輸出門的計(jì)算公式如下：

其中W*、U*與b*均為神經(jīng)網(wǎng)絡(luò)的可學(xué)習(xí)參數(shù)，* ∈{c,i,f,o} 。

在水文地質(zhì)領(lǐng)域，LSTM 模型廣泛應(yīng)用于降雨徑流模擬預(yù)測、洪水預(yù)測以及地下水位預(yù)測等水文模型的預(yù)測預(yù)報問題（Ni et al., 2020 ；Xiang et al.,2020 ；Zhang et al., 2018）。本研究將LSTM 模型的時間記憶和預(yù)測能力應(yīng)用于空間序列問題。我們通過PyTorch 平臺開發(fā)并訓(xùn)練輸入了不同特征序列的LSTM 模型。經(jīng)典機(jī)器學(xué)習(xí)模型與LSTM 模型的預(yù)測性能通過決定系數(shù)R2和均方根誤差RMSE 兩個指標(biāo)來評價，定義如下：

其中n是多孔介質(zhì)的樣本量，與分別表示樣本i的真實(shí)滲透率與預(yù)測滲透率，表示所有的平均值。

2 模型結(jié)果與討論

2.1 模型應(yīng)用方法

首先，令X1，X2，X3，…，XN表示多孔介質(zhì)樣本集。向量Xn=[x1x2…x400]表示案例1 中的第n個樣本，Xn=[x1x2…x256]則表示案例2 中的第n個樣本。向量Xn即通過上文特征提取方法取得的孔隙結(jié)構(gòu)特征序列，向量的元素數(shù)等于多孔介質(zhì)的切片數(shù)，同時也是LSTM 模型的時間（空間）步長。令Yn（n=1，2，3，…，N）表示第n個樣本的滲透率。其次，將多孔介質(zhì)的特征序列X=[X1X2X3…XN]與滲透率標(biāo)簽值y=[Y1Y2Y3…YN]按照7∶3的比例劃分為訓(xùn)練集（Xtrain，ytrain）與測試集（Xtest，ytest）。再次，應(yīng)用訓(xùn)練集（Xtrain，ytrain）對經(jīng)典機(jī)器學(xué)習(xí)模型進(jìn)行交叉驗(yàn)證與網(wǎng)格化參數(shù)搜索以確定相對最佳的超參數(shù)。對于LSTM 模型，應(yīng)當(dāng)確定數(shù)據(jù)的輸入結(jié)構(gòu)、神經(jīng)網(wǎng)絡(luò)的結(jié)構(gòu)、損失函數(shù)以及優(yōu)化器等等。經(jīng)過數(shù)次測試可以確定適當(dāng)?shù)呐看笮。╞atch size）、epoch、學(xué)習(xí)率以及其他超參數(shù)。LSTM 模型在每個時刻（空間中的位置）輸入向量Xn的一個元素。代表同一多孔介質(zhì)的不同Xn向量可合并為一個多維數(shù)組，此時LSTM 模型在每個時刻（空間中的位置）的輸入是來自不同 向量的元素組合，即切片孔隙結(jié)構(gòu)特征參數(shù)的組合。這種處理方式在案例2 中進(jìn)行了展示。而后，應(yīng)用搜索到的超參數(shù)訓(xùn)練模型。訓(xùn)練集輸入至LSTM 網(wǎng)絡(luò)中，網(wǎng)絡(luò)的權(quán)重在訓(xùn)練過程中由反向傳播算法與優(yōu)化器不斷調(diào)整。最后，將測試集（Xtest，ytest）代入訓(xùn)練完成的模型中，并以決定系數(shù)R2和均方根誤差RMSE 定量描述模型的預(yù)測性能。

2.2 案例1：由球顆粒填充的多孔介質(zhì)

本案例樣本集包含了1000 個利用沉降法生成的多孔介質(zhì)，相應(yīng)的滲透率通過LBM 計(jì)算得出。多孔介質(zhì)的尺寸為128×128×400 立方體素，體素尺寸為20 μm。滲透率均勻分布在33 Darcy 與800 Darcy 之間，平均值為416 Darcy，標(biāo)準(zhǔn)差為220 Darcy。三維多孔介質(zhì)的孔隙度范圍為0.27~0.38，平均值為0.34，標(biāo)準(zhǔn)差為0.03。本研究定義三維多孔介質(zhì)的比表面積為孔隙—顆粒邊界處的體素數(shù)與孔隙體素數(shù)之比，其范圍為0.25~0.62。應(yīng)用6-鄰域（1-連接）模式，我們還計(jì)算了多孔介質(zhì)的三維歐拉數(shù)以表示其連通性。本案例多孔介質(zhì)的孔隙空間結(jié)構(gòu)參數(shù)如表4 所示。應(yīng)用特征提取方法，含有6553600 個體素的多孔介質(zhì)被壓縮成維數(shù)為400 的向量，也即特征序列，數(shù)據(jù)壓縮率為99.994%。我們將孔隙度序列、比表面積序列和歐拉數(shù)序列分別輸入至經(jīng)典機(jī)器學(xué)習(xí)模型與LSTM 模型中。為促進(jìn)模型收斂，訓(xùn)練過程中我們對滲透率值取log。經(jīng)典機(jī)器學(xué)習(xí)模型通過Python scikit-learn 進(jìn)行訓(xùn)練與測試。模型的超參數(shù)通過交叉驗(yàn)證與網(wǎng)格搜索得到，表現(xiàn)最佳預(yù)測性能的超參數(shù)組合如表1、表2 和表3 所示。

表1 k近鄰模型的超參數(shù)Table 1 Hyperparameters of k-nearest neighbor model

表2 隨機(jī)森林模型的超參數(shù)Table 2 Hyperparameters of random forest model

表3 支持向量機(jī)回歸模型的超參數(shù)Table 3 Hyperparameters of support vector machine regression model

表4 案例1多孔介質(zhì)的屬性Table 4 Properties of porous media in case 1

各個模型的預(yù)測結(jié)果如圖8 所示。在經(jīng)典機(jī)器學(xué)習(xí)模型中，我們比較輸入不同特征序列的模型各自的預(yù)測結(jié)果。輸入比表面積序列的k近鄰、隨機(jī)森林以及支持向量機(jī)回歸均取得了各自最佳的預(yù)測結(jié)果，R2分?jǐn)?shù)分別為0.981、0.979 和0.972，RMSE 損失分別為0.041、0.043 和0.050，而輸入孔隙度序列的k近鄰和隨機(jī)森林則取得了各自最差的預(yù)測結(jié)果，R2分?jǐn)?shù)分別為0.974 和0.948。輸入歐拉數(shù)序列的支持向量機(jī)回歸同樣表現(xiàn)不佳，R2分?jǐn)?shù)為0.965，RMSE 損失為0.056。

圖8 案例1中模型的預(yù)測滲透率（logD）和真實(shí)滲透率（logD）的對比Fig. 8 Comparisons of predicted permeability (logD) and true permeability (logD) of models in case 1

與上述經(jīng)典機(jī)器學(xué)習(xí)模型相比，長短期記憶神經(jīng)網(wǎng)絡(luò)LSTM 模型的預(yù)測能力有明顯提升。我們將3 種特征序列分別輸入至LSTM 模型中。由PyTorch平臺搭建的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)設(shè)置如下。當(dāng)神經(jīng)網(wǎng)絡(luò)的超參數(shù)batch first 設(shè)置為True且輸入尺寸設(shè)置為1 時，輸入張量的大小為batch size×sequence length×1?？紤]到樣本量較少，故選擇較小的batch size 更為合適，經(jīng)嘗試后batch size 設(shè)置為16。且本案例的多孔介質(zhì)其切片數(shù)為400，因此輸入張量的大小為16×400×1。本模型將4 個LSTM 層堆疊在一起，參數(shù)num_layers 設(shè)置為4，并將hidden_size，即隱藏狀態(tài)h的特征數(shù)量設(shè)置為2048。LSTM模型將其學(xué)習(xí)到的孔隙結(jié)構(gòu)序列知識保存在外部狀態(tài)ht中。如圖7 所示，我們將序列中最后一個位置的ht，即hL，輸入至兩層全連接神經(jīng)網(wǎng)絡(luò)，網(wǎng)絡(luò)節(jié)點(diǎn)依次為2048 和1。輸入孔隙度序列、比表面積序列與歐拉數(shù)序列的LSTM 模型，初始學(xué)習(xí)率分別設(shè)置為1e-5、1e-4 以及1e-3。訓(xùn)練過程中我們選擇Adam 優(yōu)化器并將 epoch 設(shè)置為500，當(dāng)epoch為50、150、250、350 時，學(xué)習(xí)率分別固定為初始學(xué)習(xí)率的0.1、0.05、0.025 和0.01 倍。網(wǎng)絡(luò)的其余超參數(shù)采用默認(rèn)值。本研究使用一塊NVIDIA Tesla v100 16GB GPU 來訓(xùn)練模型，訓(xùn)練過程持續(xù)4 小時。如圖8 所示，輸入孔隙度序列、比表面積序列和歐拉數(shù)序列的LSTM 模型，其R2分?jǐn)?shù)分別為0.977、0.985 和0.985，RMSE 損失分別為0.045、0.036 與0.036，預(yù)測結(jié)果優(yōu)于經(jīng)典機(jī)器學(xué)習(xí)模型。

上述結(jié)果說明了以下兩點(diǎn)：（1）相比于輸入其他兩種特征序列，每個輸入比表面積序列的模型均可提供較好的預(yù)測結(jié)果，這表明比表面積序列能夠更好地表征多孔介質(zhì)的孔隙空間結(jié)構(gòu)；（2）具備序列處理能力的長短期記憶神經(jīng)網(wǎng)絡(luò)LSTM，可以更有效地建立特征序列與滲透率之間的映射關(guān)系，取得比經(jīng)典機(jī)器學(xué)習(xí)模型更高的預(yù)測分?jǐn)?shù)。

2.3 案例2：來自DeePore數(shù)據(jù)集的多孔介質(zhì)

Rabbani 等(2020) 開發(fā)了一套名為DeePore 的深度學(xué)習(xí)工作流程，能夠基于顯微斷層掃描圖像快速估計(jì)出多孔介質(zhì)的滲透率。DeePore 多孔介質(zhì)數(shù)據(jù)集基于60 個真實(shí)巖石的顯微斷層結(jié)構(gòu)，通過改變孔隙度與圖像增強(qiáng)等方法將樣本量擴(kuò)充至17700。他們假定多孔介質(zhì)為各向同性，并采用孔隙網(wǎng)絡(luò)模型（PNM）計(jì)算出x、y和z方向上滲透率的算術(shù)平均值，其單位為像素的平方，像素尺寸為5 μm。 此外，他們提取多孔介質(zhì)的三個垂直中平面并分別為它們計(jì)算固相與孔隙之間的距離變換，將結(jié)果堆疊后再用作2D CNN 的初始特征圖，該模型的R2分?jǐn)?shù)為0.9。

本案例從DeePore 多孔介質(zhì)數(shù)據(jù)集中隨機(jī)選擇5000 個尺寸為256×256×256 立方體素的樣本。滲透率的范圍為0.001~134 Darcy，平均值為10.6 Darcy，標(biāo)準(zhǔn)差為15.10 Darcy；孔隙度的數(shù)值范圍為0.1~0.45，平均值為0.27，標(biāo)準(zhǔn)差為0.10。其余的孔隙空間結(jié)構(gòu)參數(shù)如表5 所示。

表5 案例2多孔介質(zhì)的屬性Table 5 Properties of porous media in case 2

本案例依然使用交叉驗(yàn)證與網(wǎng)格搜索方法來獲取經(jīng)典機(jī)器學(xué)習(xí)模型的超參數(shù)。與案例1 的不同之處在于LSTM 模型的輸入，我們此次合并三種特征序列來共同表征多孔介質(zhì)的孔隙空間結(jié)構(gòu)，其中歐拉數(shù)序列進(jìn)行了標(biāo)準(zhǔn)化。每一個具有16777216 個體素的DeePore 多孔介質(zhì)樣本被壓縮成維數(shù)為256 的向量，或是尺寸為256×3 的數(shù)組。模型的預(yù)測結(jié)果如圖9 所示。在三種經(jīng)典機(jī)器學(xué)習(xí)模型中，輸入比表面積序列的k近鄰取得了最高的R2分?jǐn)?shù)0.963 與最低的RMSE 損失0.119。隨機(jī)森林和支持向量機(jī)回歸在輸入比表面積序列時同樣取得了相對較好的預(yù)測結(jié)果。我們由此認(rèn)為，比表面積序列更適合表征多孔介質(zhì)的孔隙空間結(jié)構(gòu)。LSTM 模型的網(wǎng)絡(luò)層數(shù)與案例1 相同，而hidden_size 增加至4096。初始學(xué)習(xí)率與batch size 分別設(shè)置為1e-4 與32。全連接層的結(jié)構(gòu)以及其他相關(guān)超參數(shù)的設(shè)置與案例1 保持了一致。我們使用一塊NVIDIA TITAN RTX 24GB GPU來訓(xùn)練模型，訓(xùn)練過程持續(xù)約26 小時。經(jīng)過500 個epoch 后，LSTM 模型的R2達(dá)到0.984，而RMSE 為0.077。與經(jīng)典機(jī)器學(xué)習(xí)模型中結(jié)果最好的k 近鄰相比，LSTM 模型的R2提升了2.18%，RMSE 降低了35.29%。此外，LSTM 模型還可以合并更多的多孔介質(zhì)孔隙結(jié)構(gòu)序列數(shù)據(jù)使得模型的魯棒性進(jìn)一步提升，而這是經(jīng)典機(jī)器學(xué)習(xí)模型無法做到的。

圖9 案例2中模型的預(yù)測滲透率（px2）和真實(shí)滲透率（px2）的對比Fig. 9 Comparisons of predicted permeability (px2) and true permeability (px2) of models in case 2

2.4 討論

本研究提出，由多孔介質(zhì)切片的孔隙度、比表面積和歐拉數(shù)所組成的特征序列可表征孔隙空間的結(jié)構(gòu)。將這三種特征序列輸入至經(jīng)典機(jī)器學(xué)習(xí)模型中，模型可以在數(shù)秒鐘內(nèi)完成訓(xùn)練并給出可靠的預(yù)測結(jié)果。在案例1 中，輸入比表面積序列的模型表現(xiàn)最好，輸入歐拉數(shù)序列的模型次之，而輸入孔隙度序列的模型表現(xiàn)最差；在案例2 中，輸入比表面積序列的模型表現(xiàn)最佳，其次是輸入孔隙度序列的模型，而輸入歐拉數(shù)序列的模型表現(xiàn)最差。由此可見，比表面積序列相對于另外兩種特征序列，能夠更好地表征多孔介質(zhì)的滲透能力，而歐拉數(shù)序列雖然在案例1 中表征能力相對于孔隙度序列要好一些，但對于孔隙空間結(jié)構(gòu)更加復(fù)雜、隨機(jī)的案例2 樣本，歐拉數(shù)序列便無法有效地表征多孔介質(zhì)的滲透能力。深度學(xué)習(xí)長短期記憶神經(jīng)網(wǎng)絡(luò)LSTM模型在兩個案例中的預(yù)測性能優(yōu)于經(jīng)典機(jī)器學(xué)習(xí)模型，訓(xùn)練時間分別為4 小時和26 小時。我們將輸入特征序列的LSTM 模型與直接輸入數(shù)字圖像的3D CNN 模型在同一多孔介質(zhì)數(shù)據(jù)集上進(jìn)行比較，兩者的預(yù)測結(jié)果基本一致，但前者的訓(xùn)練效率更高：具有4 個卷積層、卷積核尺寸為4×4×4 且特征圖輸出尺寸為64×8×8 的3D CNN 模型，在案例2中取得的R2為0.987，RMSE 為0.068。我們將3D CNN 模型的epoch 設(shè)置為150，訓(xùn)練時間持續(xù)66小時，平均每個epoch 用時26.4 分鐘，而LSTM 模型平均每epoch 僅需3.12 分鐘。

兩個案例的多孔介質(zhì)樣本在形狀尺寸、孔隙結(jié)構(gòu)和滲透率計(jì)算方法上存在顯著差異，因此我們采用遷移學(xué)習(xí)方法來驗(yàn)證LSTM 模型的泛化性能。我們用案例1 較少的樣本集對案例2 中訓(xùn)練完成的LSTM 模型進(jìn)行重新訓(xùn)練：應(yīng)用300 個訓(xùn)練樣本預(yù)測出其余700 個樣本的滲透率值。滲透率的單位為保持統(tǒng)一，已由Darcy 轉(zhuǎn)化為像素的平方。訓(xùn)練時長為4.65 小時，R2分?jǐn)?shù)為0.959，RMSE 損失為0.109?？紤]到案例2 由PNM 計(jì)算得到的滲透率與真實(shí)值之間存在5%~30%的簡化誤差（Rabbani et al., 2020），LSTM 模型的泛化能力是相當(dāng)好的。我們認(rèn)為，輸入特征序列的長短期記憶神經(jīng)網(wǎng)絡(luò)LSTM 模型在前期應(yīng)用大量樣本訓(xùn)練完成后，對于從未見過的、陌生的多孔介質(zhì)樣本集，可以借助樣本集的少量訓(xùn)練樣本為其余樣本提供準(zhǔn)確的滲透率預(yù)測結(jié)果。此外，深度學(xué)習(xí)模型的預(yù)測性能與訓(xùn)練樣本量呈正比是當(dāng)前公認(rèn)的事實(shí)，我們可以預(yù)見，如果增加訓(xùn)練樣本量使LSTM 模型得到充分的再訓(xùn)練，預(yù)測準(zhǔn)確度會有明顯提升。

3 結(jié)論

本文提出了一種直觀的方法來提取多孔介質(zhì)切片的孔隙結(jié)構(gòu)特征并將它們組成向量。向量中的元素保留了切片的空間連續(xù)性，故我們稱此類向量為特征序列。本研究提取多孔介質(zhì)每一張切 片的孔隙度、比表面積和歐拉數(shù)，分別構(gòu)建孔隙度序列、比表面積序列與歐拉數(shù)序列來表征多孔介質(zhì)的三維孔隙空間，而后應(yīng)用3 個經(jīng)典機(jī)器學(xué)習(xí)模型與深度學(xué)習(xí)長短期記憶神經(jīng)網(wǎng)絡(luò)LSTM 模型來預(yù)測多孔介質(zhì)的滲透率。我們設(shè)置了兩種孔隙結(jié)構(gòu)差異明顯的多孔介質(zhì)數(shù)據(jù)集作為研究案例。研究結(jié)果表明，輸入比表面積序列的模型相較于輸入另外兩種特征序列的模型能夠取得相對更好的預(yù)測結(jié)果，說明比表面積序列更能代表多孔介質(zhì)的孔隙空間。此外，具備序列信息處理能力的LSTM 模型，其預(yù)測結(jié)果相比于經(jīng)典機(jī)器學(xué)習(xí)模型更好。在案例1 中，輸入比表面積序列的LSTM 模型其R2分?jǐn)?shù)為0.985，RMSE損失為0.036。在案例2 中，輸入特征序列數(shù)組的LSTM 模型其R2分?jǐn)?shù)為0.984，RMSE 損失為0.077。用案例2 樣本集訓(xùn)練完成的LSTM 模型，當(dāng)應(yīng)用于案例1 時，僅需要少量樣本對模型進(jìn)行再訓(xùn)練，同樣可以給出良好的預(yù)測結(jié)果。

就目前而言，大多數(shù)預(yù)測多孔介質(zhì)滲透率的機(jī)器學(xué)習(xí)框架都以CNNs 為主，但三維多孔介質(zhì)的數(shù)據(jù)量、神經(jīng)網(wǎng)絡(luò)的復(fù)雜度以及GPU 的內(nèi)存限制使得該方法計(jì)算成本高昂。本研究提出的特征提取方法能夠有效降低多孔介質(zhì)的數(shù)據(jù)量，使機(jī)器學(xué)習(xí)模型在保持較高預(yù)測準(zhǔn)確度的同時，大幅降低計(jì)算成本。此外，長短期記憶神經(jīng)網(wǎng)絡(luò)LSTM 模型能夠從特征序列中學(xué)習(xí)切片的空間連續(xù)性，這對于今后研究多孔介質(zhì)的三維重建尤為重要。