張傲克, 錢宇航, 齊 紅

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

0 引 言

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

1 Wasserstein GAN

1.1 GAN

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

(1)

圖1 GAN訓(xùn)練流程Fig.1 The training process of GAN

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

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

1.2 Wasserstein距離

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

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

(2)

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

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

(3)

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

(4)

1.3 Lipschitz連續(xù)

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

2 評估標準

2.1 兩點相關(guān)函數(shù)

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

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

(5)

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

(6)

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

(7)

2.2 Minkowski泛函

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

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

χ=V-E+F-O

(12)

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

2.3 滲透率

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

(13)

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

φeff=Vflow/V

(14)

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

μΔv=p

(15)

(16)

2.4 后期圖像處理

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

B=2W

(17)

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

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

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

3 基于WGAN的三維重建

3.1 網(wǎng)絡(luò)架構(gòu)

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

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

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

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

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

4 基本流程

4.1 Berea砂巖

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

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

4.2 Beadpack

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

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

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

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

(18)

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

圖7 訓(xùn)練流程示意圖Fig.7 The schematic diagram of the training process

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

圖8 生成圖像與迭代次數(shù)關(guān)系Fig.8 Generate a graph of the relationship between the image and the number of iterations

4.3 數(shù)據(jù)分析

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

圖9 Berea砂巖生成圖像的徑向平均兩點相關(guān)函數(shù)值與3個笛卡爾坐標方向的兩點相關(guān)函數(shù)值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的模型在孔隙率、 比表面積和函數(shù)擬合等方面均具有較好的效果。而在圖10中, 筆者采用箱線圖對原始圖像與合成圖像進行分析。箱線圖的上下兩條線表示最大值和最小值, 而長方形的上下兩邊分別表示上四分位數(shù)和下四分位數(shù), 即數(shù)據(jù)經(jīng)過排序后處于25%和75%的位置。方格中間的虛線表示數(shù)據(jù)的中位數(shù)。

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

圖10 原始圖像與生成圖像的Minkowski functional比較結(jié)果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的三維重建模型能廣泛適用于不同數(shù)據(jù)集, 采用Beadpack進行測試(見圖12)。通過對其歸一化兩點相關(guān)函數(shù)的對比圖像可知, 其在笛卡爾坐標方向和徑向方向均具有較好的擬合效果。從而進一步表明了該模型的泛化能力。

圖12 Beadpack生成圖像徑向平均兩點相關(guān)函數(shù)值與3個笛卡爾坐標方向的兩點相關(guān)函數(shù)值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 原始重建圖像和經(jīng)過平滑處理后的結(jié)果Fig.13 The original reconstructed image and the result of smoothing

5 結(jié) 語

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

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

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