基于WGAN的多孔材料三維重建

張傲克, 錢(qián)宇航, 齊 紅

(吉林大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院, 長(zhǎng)春 130012)

0 引 言

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

1 Wasserstein GAN

1.1 GAN

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

(1)

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

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

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

1.2 Wasserstein距離

對(duì)式(1)進(jìn)一步分析可得

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

(2)

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

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

(3)

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

(4)

1.3 Lipschitz連續(xù)

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

2 評(píng)估標(biāo)準(zhǔn)

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

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

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

(5)

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

(6)

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

(7)

2.2 Minkowski泛函

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

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

χ=V-E+F-O

(12)

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

2.3 滲透率

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

(13)

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

φeff=Vflow/V

(14)

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

μΔv=p

(15)

(16)

2.4 后期圖像處理

根據(jù)Nyquist準(zhǔn)則可知, 對(duì)一個(gè)無(wú)噪聲低通信道, 最高碼元傳輸速率B滿(mǎn)足關(guān)系

B=2W

(17)

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

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

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

3 基于WGAN的三維重建

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

判別器中將三維卷積神經(jīng)網(wǎng)絡(luò)與三維批處理層作為一個(gè)模塊, 批處理層用于加快神經(jīng)網(wǎng)絡(luò)的收斂速度同時(shí)提高穩(wěn)定性。除最后一個(gè)模塊外三維卷積神經(jīng)網(wǎng)絡(luò)中卷積核的大小為(4,4,4), 步長(zhǎng)為(2,2,2), 并以(1,1,1)的大小進(jìn)行填充, 以保留邊界信息。同時(shí), 其余模塊中均有線(xiàn)性整流函數(shù)作為激活函數(shù), 以引入更多的非線(xiàn)性因素增強(qiáng)模型的擬合能力。在最后一個(gè)模塊中, 不使用激活函數(shù)與批處理層, 并且步長(zhǎng)為(1,1,1),不進(jìn)行填充, 卷積核的大小仍為(4,4,4)。對(duì)卷積神經(jīng)網(wǎng)絡(luò)和批處理層, 其超參數(shù)均采用默認(rèn)值。最終將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ù), 邊長(zhǎng)表示三維圖像的邊長(zhǎng)。生成器除最后一個(gè)模塊外將三維反卷積神經(jīng)網(wǎng)絡(luò)與三維批處理層作為一個(gè)模塊, 三維反卷積神經(jīng)網(wǎng)絡(luò)中卷積核的大小均為(4,4,4), 步長(zhǎng)為(2,2,2), 以(1,1,1)的大小進(jìn)行填充, 并且每個(gè)模塊中均使用線(xiàn)性整流函數(shù)。其中在第1個(gè)模塊中, 僅有一層三維反卷積神經(jīng)網(wǎng)絡(luò), 卷積核大小仍為(4,4,4), 但步長(zhǎng)為(1,1,1), 并且不進(jìn)行任何填充。最后, 通過(guò)雙曲正切函數(shù)處理輸出結(jié)果。反卷積神經(jīng)網(wǎng)絡(luò)和批處理層均采用默認(rèn)超參數(shù)值。生成器的示意圖以及輸入輸出通道數(shù)如圖4所示。由于輸入數(shù)據(jù)維數(shù)為4。因此, 將每層網(wǎng)絡(luò)的寬度表示其通道數(shù), 邊長(zhǎng)表示三維圖像的邊長(zhǎng)。

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

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

4 基本流程

4.1 Berea砂巖

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

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

4.2 Beadpack

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

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

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

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

(18)

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

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

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

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

圖9 Berea砂巖生成圖像的徑向平均兩點(diǎn)相關(guān)函數(shù)值與3個(gè)笛卡爾坐標(biāo)方向的兩點(diǎn)相關(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中, 筆者采用箱線(xiàn)圖對(duì)原始圖像與合成圖像進(jìn)行分析。箱線(xiàn)圖的上下兩條線(xiàn)表示最大值和最小值, 而長(zhǎng)方形的上下兩邊分別表示上四分位數(shù)和下四分位數(shù), 即數(shù)據(jù)經(jīng)過(guò)排序后處于25%和75%的位置。方格中間的虛線(xiàn)表示數(shù)據(jù)的中位數(shù)。

通過(guò)結(jié)果可看出, 由該模所生成的圖像具有與原始圖像類(lèi)似的拓?fù)湫再|(zhì), 即具有相似的連通性。圖11給出了基于WGAN的三維重建圖形的滲透率。從圖11可看出, 生成圖像的有效孔隙率與滲透率的關(guān)系表現(xiàn)出的相關(guān)性與原始圖形類(lèi)似, 并且坐標(biāo)點(diǎ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

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

圖12 Beadpack生成圖像徑向平均兩點(diǎn)相關(guān)函數(shù)值與3個(gè)笛卡爾坐標(biāo)方向的兩點(diǎn)相關(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中的中值濾波、 雙邊濾波和線(xiàn)性濾波進(jìn)行比較, 最終選用中值濾波對(duì)圖像進(jìn)行處理。從圖13中可看出, 圖像的清晰度明顯提高, 巖石表面更加連續(xù)、 光滑。

圖13 原始重建圖像和經(jīng)過(guò)平滑處理后的結(jié)果Fig.13 The original reconstructed image and the result of smoothing

5 結(jié) 語(yǔ)

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

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

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