張傲克, 錢宇航, 齊 紅

(吉林大學 計算機科學與技術學院, 長春 130012)

0 引 言

多孔材料在工程應用中具有重要價值[1-2]。根據“結構決定性質, 性質決定用途”原則, 如何采用適當的方法探究其內在聯(lián)系是人們關注的重點。目前, 通過計算機斷層掃描技術[3]可直接獲取3D圖像, 但該方法存在高成本、 分辨率低等缺點。 已有的三維重建方法, 包括基于優(yōu)化的方法[4-5]、 直接取樣法[6-7]和多點統(tǒng)計法[8-9]等, 對三維重建問題計算時間復雜度、 空間復雜度和準確率仍然是衡量模型性能的重要標準, 而上述方法在其不同方面均存在一定的缺陷。但通過這些算法可知, 提取一個圖像的形態(tài)學特征, 并利用該潛在空間中的向量進行三維重建是解決該問題的有效手段之一。筆者以Berea砂巖和Beadpack為數據集, 將WGAN(Wasserstein Generative Adversarial Networks)改造為用于處理三維重建問題的深度學習模型, 通過對相關評價標準的測試得知該深度學習模型具有較高準確性與較低時間復雜性等特點, 并具有較好的泛化能力。

1 Wasserstein GAN

1.1 GAN

在傳統(tǒng)的機器學習算法中, 通常僅考慮現(xiàn)有的數據模式, 而無法生成新的事物。Goodfellow等[10]提出的生成對抗網絡是目前非常流行的無監(jiān)督學習方法之一, 該方法通過學習數據的隱式分布生成新的事物, 并且具有極好的效果[11]。在該深度學習模型中使用一個生成器與一個判別器通過相互博弈實現(xiàn)分布近似。其中, 生成器通過給定輸入生成需要的結果, 而判別器則用于判斷所輸入內容是否為真, 生成器和判別器進行以下最小-最大(min-max)博弈

(1)

圖1 GAN訓練流程Fig.1 The training process of GAN

通過式(1)可知, 在訓練過程中首先固定生成器, 并通過訓練使判別器達到最優(yōu), 即讓判別器輸入生成數據時判別為假, 輸入真實數據時判別為真。再固定判別器, 訓練生成器使其達到最優(yōu), 即讓判別器在輸入生成數據時判別為真。重復上述過程, 進行迭代訓練達到或接近納什均衡點, 此時GAN處于收斂狀態(tài)。

圖1為GAN的訓練流程。但GAN存在收斂速度慢、 模式崩潰和過度泛化等問題[12], 而Wasserstein GAN則可在一定程度上解決上述問題。

1.2 Wasserstein距離

對式(1)進一步分析可得

V(G,D*)=2J(pr(x)‖pg(x))-2log 2

(2)

而J散度又存在飽和區(qū)間, 同時根據“當Pr與P-g的支撐集是高維空間中的低維流形時,pr與pg重疊部分測度為0的幾率為1”[13]結論, 可知GAN在初始化狀態(tài)下就可能出現(xiàn)嚴重問題。因此, 需要對J散度進行替代。

最優(yōu)輸運起源于Monge問題, 并由Kantorovich通過對偶理論加以發(fā)展, 進而提出了Kantorovich-Rubinstein距離等重要概念[14]。最優(yōu)輸運的相關理論在流體力學、 生物醫(yī)學和人工智能等領域中有重要應用價值。其中, Wasserstein距離[15](簡稱W距離)是最優(yōu)輸運理論中的重要概念, 其代表諸多輸運方法中的最小代價。W距離為

(3)

通過對W距離分析可知, 其具有交換性、 非負性和指標性, 因此解決了J散度的相關問題。并且對生成圖像的概率分布, 可通過求解線性規(guī)劃的對偶問題解決。最終損失函數可更新為[13]

(4)

1.3 Lipschitz連續(xù)

通過式(4)可知,fω(x)需要滿足k-Lipschitz連續(xù)。為滿足這個條件, WGAN中使用了一種稱為權重裁剪的策略。由于, 在WGAN中Lipschitz連續(xù)并不關注于梯度的具體大小, 而是要保證梯度在合理范圍內。因此, 通過將fω(x)中的ω限制在[-c,c], 進而使函數的輸出值控制在一定范圍內, 通過這種方式保證輸出值的變化在一個有限的范圍內。

2 評估標準

2.1 兩點相關函數

兩點相關函數計算了隨機距離為‖r‖的兩點落在同一個相中的概率[15]。其公式為

S2(r)=P(x∈P,x+r∈P),x,r∈n

(5)

其中S2(0)是孔隙率φ, 孔隙率描述了多孔材料存儲流體的能力, 并且當r→∞時,S2趨近于φ2。通常由于其各向異性, 需沿x,y,z3個笛卡爾坐標方向計算S2(r)的值, 以及其徑向平均值。同時可通過函數圖像獲取更多信息。由文獻[16]可知, 在一個各向同性的多孔材料中, 比表面積Sv與S2滿足

(6)

其中比表面積的大小與物體的吸附能力有關。并通過式(5)可計算出顆粒相與多孔相中的弦長。其計算公式為[17]

(7)

2.2 Minkowski泛函

采用Minkowski泛函對多孔材料進行數據分析是通過形態(tài)學特征刻畫多孔材料特性的基本方法之一。Minkowski 泛函起源于有關凸集的理論, 積分幾何中的Hadwiger理論表明可用d+1維參數描述d維的嵌入空間。因此, 對三維圖像可采用4個Minkowski泛函刻畫其形態(tài)學特征[17]。

式(8)是孔隙率, 表示孔隙所占體積與重建對象的總體體積之比。式(9)用于計算比表面積。式(10)是物體表面的平均曲率, 其中r1和r2是面積元dS的最小和最大曲率半徑。式(11)與Euler特性密切相關, 它表示了總體曲率。但在實際計算過程中考慮Euler數, 式(11)可表示如下

χ=V-E+F-O

(12)

其中V是頂點數,E是邊數,F是面數,O是對象的數量。Euler數表示了介質的連通性等性質。但由于計算工具受限, 在實際數據分析過程中沒有計算表面的平均曲率。其余參數通過image J的MorphoLib J進行計算。

2.3 滲透率

滲透率[17]表示流體穿過多孔材料的能力, 高滲透性材料意味著流體更容易通過相應介質。滲透率可看做Darcy定律在多孔材料中的應用。Darcy定律表述如下

(13)

其中U是平均流體速度,p0是壓力梯度,k是介質的滲透率,μ是流體的動態(tài)粘度。在計算中, 需要考慮具有一定流量的動態(tài)連通部分, 因此需要計算有效孔隙率

φeff=Vflow/V

(14)

滲透率有很多實際應用, 如從多孔巖石中提取石油, 以計算相關數據。為計算原始三維圖像以及生成的三維圖像中的滲透率, 可采用描述粘性不可壓縮流體動量守恒的運動方程, 即Navier-Stokes方程, 關系式為

μΔv=p

(15)

(16)

2.4 后期圖像處理

根據Nyquist準則可知, 對一個無噪聲低通信道, 最高碼元傳輸速率B滿足關系

B=2W

(17)

圖2 CT圖像的局部放大結果Fig.2 The local magnification result of CT image

其中W是帶寬。在圖形學中, 經常遇到鋸齒、 摩爾紋和車輪效應等走樣現(xiàn)象。走樣是指低頻采樣造成的信息失真。由于CT掃描過程存在走樣現(xiàn)象(見圖2), 因此考慮相應方法對圖像進行后期處理。由于已經完成采樣, 產生了頻譜的混疊, 所以只能通過濾波等操作處理相應的圖像, 使生成的巖石表面更加平滑、 清晰。

采用opencv[19]中的相應函數對圖像進行處理, 包括中值濾波、 雙邊濾波和線性濾波。不同的濾波器對圖像平滑的方式不同, 例如中值濾波是一種非線性濾波器, 它根據鄰域窗口中像素值的中值確定該點處的像素值。通過對每種方法進行嘗試, 最終根據實際情況選擇合適的處理結果作為最終的生成圖像。

3 基于WGAN的三維重建

3.1 網絡架構

判別器中將三維卷積神經網絡與三維批處理層作為一個模塊, 批處理層用于加快神經網絡的收斂速度同時提高穩(wěn)定性。除最后一個模塊外三維卷積神經網絡中卷積核的大小為(4,4,4), 步長為(2,2,2), 并以(1,1,1)的大小進行填充, 以保留邊界信息。同時, 其余模塊中均有線性整流函數作為激活函數, 以引入更多的非線性因素增強模型的擬合能力。在最后一個模塊中, 不使用激活函數與批處理層, 并且步長為(1,1,1),不進行填充, 卷積核的大小仍為(4,4,4)。對卷積神經網絡和批處理層, 其超參數均采用默認值。最終將1×64×64像素的圖像轉換為1×1×1的距離值。判別器的示意圖以及相應的輸入輸出通道數如圖3所示。

圖3 判別器示意圖Fig.3 The schematic diagram of the discriminator

由于輸入數據維數為4。因此, 用每層網絡的寬度表示其通道數, 邊長表示三維圖像的邊長。生成器除最后一個模塊外將三維反卷積神經網絡與三維批處理層作為一個模塊, 三維反卷積神經網絡中卷積核的大小均為(4,4,4), 步長為(2,2,2), 以(1,1,1)的大小進行填充, 并且每個模塊中均使用線性整流函數。其中在第1個模塊中, 僅有一層三維反卷積神經網絡, 卷積核大小仍為(4,4,4), 但步長為(1,1,1), 并且不進行任何填充。最后, 通過雙曲正切函數處理輸出結果。反卷積神經網絡和批處理層均采用默認超參數值。生成器的示意圖以及輸入輸出通道數如圖4所示。由于輸入數據維數為4。因此, 將每層網絡的寬度表示其通道數, 邊長表示三維圖像的邊長。

圖4 生成器示意圖Fig.4 The schematic diagram of the generator

根據Arjovsky等[13]提出的權重裁剪簡單地實現(xiàn)Lipschitz連續(xù), 并將相應的參數取為±0.01。但在實際訓練中, 應將權重值取為盡可能大的值如±1 000 000, 否則將出現(xiàn)訓練失敗等情況。同時, 自適應學習率算法應采用非動量算法, 一般使用RMSProp, 該方法已被證明具有良好的效果[20]并得到廣泛應用。

4 基本流程

4.1 Berea砂巖

圖5所示為Berea砂巖, 是一種淺灰色至淺黃色, 顆粒大小為細粒到中粒的砂巖。Berea砂巖通常用于建筑材料, 同時也是石油和天然氣的重要來源。

采用的CT掃描圖像來自于文獻[18], 該圖像的原始大小為400×400×400像素, 以32為步長, 并將原始圖像分割成大小為64×64×64像素的子圖像。將得到的三維圖像作為真實數據在訓練中使用。

4.2 Beadpack

Beadpack是一種由相同大小的顆粒無序堆積而成的人造巖心。筆者采用的Beadpack原始圖像大小為500×500×500像素, 并將原始圖像以步長32進行分割。在訓練過程中, 使用子圖像作為真實圖像的數據集。

圖5 Berea砂巖的三維圖像 圖6 Beadpack三維圖像 Fig.5 Three dimensional image of Berea sandstone Fig.6 Beadpack 3D image

首先通過使用滿足高斯正態(tài)分布并且通道數為512的隨機噪聲, 固定預先經過初始化的生成器, 并生成相應數據。然后將生成數據與真實數據分別輸入判別器, 并計算相應的損失函數值, 通過誤差逆?zhèn)鞑ニ惴ǜ屡袆e器中的相應權重。參數值的更新公式如下

ωk+1=ωk-αf(ωk)

(18)

完成該輪對判別器的訓練。固定判別器, 對判別器輸入生成數據并計算相應的損失函數, 通過誤差逆?zhèn)鞑ニ惴ǜ律善鞯南鄳獏?見圖7)。

圖7 訓練流程示意圖Fig.7 The schematic diagram of the training process

圖8為生成圖像與迭代次數關系。由圖8可知, 在訓練初期, 無法觀察到具有孔隙結構的圖像, 在立方體表面只有模糊的斑點。但在迭代5 000次時, 已經能看見相對清晰的孔隙。而在10 000次～20 000次之間, 訓練效果趨于穩(wěn)定。但當訓練次數繼續(xù)增加達到40 000次左右時, 孔隙數量有一定的減少。

圖8 生成圖像與迭代次數關系Fig.8 Generate a graph of the relationship between the image and the number of iterations

4.3 數據分析

通過式(5)可計算出重建圖像的徑向平均兩點相關函數值與沿著x,y,z3個方向的歸一化兩點相關函數值(見圖9)。

圖9 Berea砂巖生成圖像的徑向平均兩點相關函數值與3個笛卡爾坐標方向的兩點相關函數值Fig.9 The radial average two-point correlation function value of the image generated by Berea sandstone and the two-point correlation function value of three Cartesian coordinate directions

由圖9可知, 基于WGAN的模型在孔隙率、 比表面積和函數擬合等方面均具有較好的效果。而在圖10中, 筆者采用箱線圖對原始圖像與合成圖像進行分析。箱線圖的上下兩條線表示最大值和最小值, 而長方形的上下兩邊分別表示上四分位數和下四分位數, 即數據經過排序后處于25%和75%的位置。方格中間的虛線表示數據的中位數。

通過結果可看出, 由該模所生成的圖像具有與原始圖像類似的拓撲性質, 即具有相似的連通性。圖11給出了基于WGAN的三維重建圖形的滲透率。從圖11可看出, 生成圖像的有效孔隙率與滲透率的關系表現(xiàn)出的相關性與原始圖形類似, 并且坐標點的分布也比較接近。

圖10 原始圖像與生成圖像的Minkowski functional比較結果Fig.10 Comparison results of Minkowski functional between original image and generated image

圖11 原始圖像與生成圖像的滲透率比較Fig.11 Comparison of permeability between original image and generated image

為驗證基于WGAN的三維重建模型能廣泛適用于不同數據集, 采用Beadpack進行測試(見圖12)。通過對其歸一化兩點相關函數的對比圖像可知, 其在笛卡爾坐標方向和徑向方向均具有較好的擬合效果。從而進一步表明了該模型的泛化能力。

圖12 Beadpack生成圖像徑向平均兩點相關函數值與3個笛卡爾坐標方向的兩點相關函數值Fig.12 Beadpack generates the radial radial average two-point correlation function value of the image and the two-point correlation function value in three Cartesian coordinate directions

4.4 后期處理

使用opencv中的中值濾波、 雙邊濾波和線性濾波進行比較, 最終選用中值濾波對圖像進行處理。從圖13中可看出, 圖像的清晰度明顯提高, 巖石表面更加連續(xù)、 光滑。

圖13 原始重建圖像和經過平滑處理后的結果Fig.13 The original reconstructed image and the result of smoothing

5 結 語

基于WGAN的三維重建模型具有較低時間復雜度。基于隨機數的三維重建模型在NVDIA 3060上完成對berea砂巖和Beadpack圖像的84 000次迭代僅需3.5～4 h左右。并且在1 h左右即可觀察到明顯效果。當訓練完成后, 便可在秒級時間單位下完成三維重建, 同時通過訓練得到的生成模型可存儲并重復利用, 而傳統(tǒng)方法則需要數小時完成三維重建, 并且這種通過顯式分布進行學習的方法的時間復雜度隨著圖像的體素的增多而變大, 因此基于深度學習的模型在效率上具有明顯優(yōu)勢。同時根據相應的性能評估標準可知, 生成圖像與真實圖像有較好的擬合效果。

通過對實驗結果進行分析可知, 該模型仍存在以下不足: 1) 該模型可能會出現(xiàn)模式崩塌等常見問題, 從而導致生成圖像的單一性; 2) 由幾何的相應知識可知, 在整個過程中使用單張二維圖像進行訓練會造成信息的缺失, 從而導致生成圖像不能無限接近于真實圖像。因此, 通過采用更多數目的非平行圖像可能會有更好的效果。

基于WGAN的三維重建模型相關參數并未進行精細調控, 例如在訓練過程中可實時關注損失函數的變化情況, 通過人工改變學習率的方法提升訓練效果。同時文中所示的權重截斷超參數的范圍為有效范圍, 但實際范圍可以進一步縮小。網絡中的參數與網絡層數量可以得到進一步地調整與縮減, 在保證高性能的同時, 能提高訓練速度, 并減小空間復雜度。另外, 批尺寸大小與圖像的大小也是值得考慮的因素, 通過找到相應參數的最優(yōu)取值, 能加快訓練速度, 提高內存利用率, 并優(yōu)化訓練效果。并且, 傳統(tǒng)方法與深度學習方法相結合也是值得考慮的方向之一?？梢钥紤]采用模擬退火算法對重建結果進行再重建, 由于通過深度學習方法可在低時間復雜度的情況下實現(xiàn)較高準確率的三維重建, 因此可以結合高效率的傳統(tǒng)三維重建算法進行優(yōu)化, 以實現(xiàn)準確率與效率的權衡。