屠紅艷，張 挺*，夏鵬飛，杜 奕

（1.上海電力大學計算機科學與技術(shù)學院，上海 200090；2.上海第二工業(yè)大學工學部，上海 201209）

0 引言

空間信息是指包含空間特征的數(shù)據(jù)集合［1-2］?？臻g信息的特征表征是地質(zhì)勘探、空間科學、石油工程和生物學等領域的重要研究課題之一［3-4］。直接獲取大范圍空間信息的成本較高，通過數(shù)值方法重建是獲取空間信息的主要手段之一。目前數(shù)值重建方法主要分為“確定”性重建方法和“不確定”性重建方法［5］。“不確定”性重建方法的不確定性主要體現(xiàn)在重建結(jié)果的隨機性。目前，作為典型的“不確定”性重建方法，多點信息統(tǒng)計法（Multiple-Point Statistics，MPS）通過提取訓練圖像（Training Image，TI）的內(nèi)在特征完成空間數(shù)據(jù)重建［6］，但在每次重建時都需要重新掃描TI而導致整體重建速度較慢。

近年來，深度學習已經(jīng)廣泛應用于與特征提取相關(guān)的研究領域，在深度生成模型方面，無監(jiān)督學習模型取得了突破性進展［7］，其中，變分自編碼器（Variational Auto-Encoder，VAE）被視為無監(jiān)督學習領域的重要工具之一，在深度生成模型領域得到了越來越多的應用［8-9］。

本文結(jié)合信息論將VAE 應用于空間信息的不確定性重建研究，將結(jié)合了信息論的VAE 方法——信息變分自編碼器（Information Variational Auto-Encoder，IVAE）與原始VAE 方法作對比，在重建效果上有一定改進；與一些傳統(tǒng)的重建方法相比，在速度上更具有優(yōu)勢，可以大大縮短訓練時間。實驗驗證了本方法的有效性。

1 本方法核心思想

1.1 變分自編碼器模型

VAE 是一種基于變分貝葉斯網(wǎng)絡結(jié)構(gòu)的深度生成模型，它使用兩種神經(jīng)網(wǎng)絡建立了兩種概率密度分布模型：一種用于對原始輸入數(shù)據(jù)進行變分推理，生成隱式變量的變分概率分布，稱為編碼器（Encoder）；另一種是根據(jù)生成的隱式變量變分概率分布，生成訓練數(shù)據(jù)的近似概率分布來重建數(shù)據(jù)，稱為解碼器（Decoder）。VAE結(jié)構(gòu)如圖1所示。

圖1 VAE結(jié)構(gòu)Fig.1 Structure of VAE

圖1中，“+”和“*”分別表示元素相加和相乘，X'表示重建結(jié)果。在編碼階段，由于給定的數(shù)據(jù)集X通常是高維的，具有復雜的依賴關(guān)系，所以VAE 首先通過編碼器學習訓練數(shù)據(jù)特征，由于誤差在傳播過程中需要經(jīng)過一個不連續(xù)的、沒有梯度的采樣層，而隨機梯度下降（Stochastic Gradient Descent，SGD）雖然可以處理隨機輸入，但不能處理隨機操作，所以需要“重新參數(shù)化”，即通過引入輔助參數(shù)ε來實現(xiàn)。ε從標準正態(tài)分布N（0，1）中采樣得到，使得中間隱變量z和高斯分布的標準差和均值σ、μ之間的關(guān)系可用式（1）表示，從而可以使用SGD 方法進行優(yōu)化。如圖1所示，σ和μ由編碼器深度神經(jīng)網(wǎng)絡計算得到，z的采樣可以通過以下方式完成［10］：

VAE 模型通過引入?yún)?shù)φ和θ分別控制編碼器和解碼器的神經(jīng)網(wǎng)絡優(yōu)化。由于真實數(shù)據(jù)X的隱變量z的分布是不可知的，為此在編碼器網(wǎng)絡中引入用φ參數(shù)化的神經(jīng)網(wǎng)絡編碼器模型qφ(z|x)來取代無法確定的真實后驗分布Pθ(z|x)，并采用KL（Kullback-Leibler）散度［11］來度量編碼器模型qφ(z|x)和真實后驗分布Pθ(z|x)的相似性，訓練優(yōu)化約束參數(shù)φ和θ，目標是最小化KL 散度，即最大化證據(jù)下界函數(shù)L，具體的相關(guān)含義見參考文獻［12］。證據(jù)下界函數(shù)定義如下：

其中：DKL表示KL散度；lb表示以2為底的對數(shù)；Pθ(x'|z)是用θ參數(shù)化的神經(jīng)網(wǎng)絡解碼器；Pθ(z)表示隱變量z的概率分布。最后將中間隱變量z輸入到解碼器模型Pθ(z)Pθ(x'|z)得到最終的空間數(shù)據(jù)重建結(jié)果X'。

1.2 香農(nóng)信息熵與費雪信息量

從上節(jié)VAE 模型的證據(jù)下界的定義可以發(fā)現(xiàn)，原始VAE模型只考慮KL 散度作為正則化懲罰項，所以很難在隱變量z的表示學習和重建結(jié)果的似然最大化之間取得平衡［13］。因為散度誤差只控制VAE 編碼網(wǎng)絡的表示方式，對解碼網(wǎng)絡沒有控制，所以提出將隨機向量X的香農(nóng)信息熵和費雪信息量相結(jié)合應用于VAE 的編碼階段和解碼階段。設X的香農(nóng)信息熵為H（X），費雪信息量為J（X）。在香農(nóng)信息熵的使用中通常將之轉(zhuǎn)化為熵權(quán)形式［14］：

其中exp 表示以e 為底的指數(shù)函數(shù)。當X表示隨機向量時，費雪信息量J(X)表示一個費雪信息量矩陣，那么N（X）與J(X)的跡之間滿足關(guān)系：

其中：K取決于分布形式和隨機向量的維數(shù)，K≥1 是當且僅當X為一個高斯隨機向量時成立，這個性質(zhì)的證明見參考文獻［15］，tr表示矩陣的跡。由式（4）可以發(fā)現(xiàn)信息量J（X）與香農(nóng)熵權(quán)N（X）之間存在一個平衡關(guān)系，可以通過香農(nóng)熵權(quán)和費雪信息間的關(guān)系來平衡似然估計和真實數(shù)據(jù)與隱變量z之間的依賴關(guān)系。

1.3 IVAE主要思想

在本節(jié)中利用信息量來平衡VAE 模型的編碼階段和解碼階段，并提出新的規(guī)則懲罰項。IVAE 模型不僅考慮KL 散度作為懲罰項，而且在編碼階段和解碼階段分別考慮其費雪信息作為懲罰項，即在費雪信息量的約束下最大化證據(jù)下界，并重新構(gòu)造證據(jù)下界L（θ，φ；X）如下：

與原證據(jù)下界式（2）相比，式（5）增加了兩個分別調(diào)節(jié)編碼器和解碼器的懲罰項，λz調(diào)節(jié)編碼器網(wǎng)絡，λx'調(diào)節(jié)解碼器網(wǎng)絡。Fx'和Fz均為正常數(shù)，分別表示解碼器網(wǎng)絡和編碼器網(wǎng)絡中期望的費雪信息值。Fx'和Fz越大，表明更傾向于用θ和φ參數(shù)化的模型進行分布估計；反之則表明弱化了分布建模，增強了香農(nóng)熵權(quán)的影響，而費雪信息的估計可以根據(jù)它的定義直接計算。

下面根據(jù)編碼器網(wǎng)絡和解碼器網(wǎng)絡分別討論模型的優(yōu)化目標。先以編碼器網(wǎng)絡為例，在編碼器網(wǎng)絡中考慮費雪信息量，則設式（5）中λx'為0，編碼器網(wǎng)絡的優(yōu)化目標Le（θ，φ；X）包含了散度誤差項和信息誤差項：

后驗分布qφ(z|x)為正態(tài)分布，Pθ(z)為標準正態(tài)分布，故式（6）中的KL散度計算為：

式（6）中的費雪信息量根據(jù)定義［12］得到：

將式（7）和（8）代入式（6），得到優(yōu)化目標為：

同理，在解碼器網(wǎng)絡中將費雪信息量與證據(jù)下界函數(shù)結(jié)合，則設式（5）中λz為0，解碼器網(wǎng)絡的優(yōu)化目標Ld包含了重建項和信息誤差項：

對比IVAE 與原始VAE，可見IVAE 的證據(jù)下界公式（5）包含一個重建項和兩個分別考慮了信息量和KL 散度的懲罰項，而VAE只考慮KL散度正則化懲罰項，IVAE通過結(jié)合費雪信息量和香農(nóng)信息熵來平衡似然估計和數(shù)據(jù)與隱變量之間的依賴關(guān)系，并提出新的規(guī)則懲罰項，使得VAE 模型結(jié)合考慮KL散度和空間數(shù)據(jù)的信息量，從而提高空間數(shù)據(jù)重建質(zhì)量。

2 空間信息重建方法的流程

2.1 算法步驟

VAE 本質(zhì)上是在傳統(tǒng)自編碼器的基礎上，將神經(jīng)網(wǎng)絡計算結(jié)果的平均值加上“高斯噪聲”，用方差神經(jīng)網(wǎng)絡來動態(tài)調(diào)整噪聲強度。引入編碼噪聲和KL 正則化使編碼到隱空間中的數(shù)據(jù)特征分布不斷地迭代優(yōu)化以接近標準高斯分布。本文將IVAE模型應用于空間信息重建，分別在編碼器網(wǎng)絡和解碼器網(wǎng)絡中將信息誤差作為約束項。具體步驟如下：

步驟1 設計訓練網(wǎng)絡結(jié)構(gòu)，初始化所有參數(shù)，應用SGD更新參數(shù)，選擇學習率。

步驟2 使用設計好的深度神經(jīng)網(wǎng)絡對訓練數(shù)據(jù)進行迭代擬合，得到訓練數(shù)據(jù)的概率分布和信息量，根據(jù)式（6）對編碼器網(wǎng)絡進行優(yōu)化。

步驟3 從N（0，1）中進行隨機采樣，將隨機采樣結(jié)果和編碼器網(wǎng)絡得到的σ、μ代入式（1）計算中間結(jié)果z，當訓練誤差達到要求時保存網(wǎng)絡參數(shù)。

步驟4 將中間結(jié)果z作為解碼器網(wǎng)絡的輸入，根據(jù)式（10）對解碼器網(wǎng)絡進行迭代優(yōu)化，對z進行逐層解碼，解碼器的解碼結(jié)果即為重建結(jié)果，保存網(wǎng)絡參數(shù)。

2.2 技術(shù)路線

首先通過數(shù)字圖像處理技術(shù)對采集的圖像的色彩空間進行分析，確定分割閾值，對分割后的圖像進行形態(tài)學處理，最終獲得孔隙幾何形態(tài)結(jié)構(gòu)較好的訓練數(shù)據(jù)。然后對訓練數(shù)據(jù)應用本文所提方法，先將數(shù)據(jù)輸入編碼器進行編碼，然后通過采樣得到中間結(jié)果，最后將中間結(jié)果輸入到解碼器中進行解碼得到最終結(jié)果。技術(shù)路線如圖2所示。

圖2 技術(shù)路線Fig.2 Technology roadmap

3 實驗與結(jié)果分析

本實驗擬采用真實巖石數(shù)據(jù)作為空間信息重建的數(shù)據(jù)源。由于真實巖石孔隙數(shù)據(jù)難以用數(shù)學公式或者某種定量的語句描述，故巖石的孔隙分布具有很強的不確定性特征，適合采用不確定性重建方法重建孔隙結(jié)構(gòu)。實驗樣品為直徑3 mm 的圓柱形砂巖，通過采用同步輻射射線掃描砂巖樣品獲得該砂巖樣本的三維空間數(shù)據(jù)，分辨率為每體素10 μm。

3.1 訓練圖像

從上述真實砂巖中截取80×80×80體素的體數(shù)據(jù)作為TI，孔隙度為0.170 5，外表面、剖面圖（X=40，Y=40，Z=40）和孔隙如圖3 所示。圖3 的砂巖中僅包含兩種狀態(tài)值：孔隙和骨架，其中藍色表示孔隙，灰色表示骨架。輸入數(shù)據(jù)集為80×80×80的三維矩陣，矩陣中的元素有兩種狀態(tài)值｛0，1｝，1 表示孔隙，0表示骨架。訓練數(shù)據(jù)共512 000（80×80×80）條。

圖3 訓練圖像Fig.3 Training image

3.2 重建結(jié)果比較

在下面的實驗中，分別采用IVAE、VAE 和其他幾種經(jīng)典的不確定性空間信息重建方法SNESIM（Single Normal Equation SIMulation）［16］、DS（Direct Sampling）［17］重建砂巖圖像并進行對比分析。重建圖像如圖4所示。

圖4 各方法的重建圖像對比Fig.4 Comparison of reconstructed images by different methods

可見四種重建方法的結(jié)果均與TI（圖3）具有相似的結(jié)構(gòu)和長連通的孔隙空間，下面將量化比較各重建結(jié)果。

3.2.1 孔隙度比較

孔隙度是指巖樣中所有孔隙空間體積之和與該巖樣體積的比值，孔隙度?的定義如下：

其中：VP表示孔隙體積，V表示巖石體積。多次重建結(jié)果的孔隙度均值在一定程度上反映了重建方法的重建質(zhì)量，孔隙度方差則表明重建方法的穩(wěn)定性。本文使用的訓練數(shù)據(jù)的孔隙度為0.170 5，分別用IVAE、VAE、SNESIM 和DS 進行30 次砂巖重建。如表1 所示，IVAE 重建的三維圖像的孔隙度更接近訓練數(shù)據(jù)，而且方差最小，說明重建質(zhì)量好而且波動性小。

表1 各方法30次重建結(jié)果的孔隙度均值與方差Tab.1 Porosity mean and variance of 30 reconstruction results by using each method

3.2.2 變差函數(shù)比較

變差函數(shù)通常作為空間數(shù)據(jù)的空間連續(xù)性的評價工具，定義如下：

其中：Var 表示方差，h表示空間距離，x和x+h分別表示空間兩點位置，Z（x）表示空間點的狀態(tài)值。

本文分別比較了TI 和IVAE、VAE、SNESIM 以及DS 的重建圖像在X、Y、Z方向的變差函數(shù)，如圖5 所示，IVAE 重建的空間數(shù)據(jù)與TI在3 個方向上（尤其是Z方向）的變差函數(shù)更為接近，表明IVAE 重建圖像的孔隙結(jié)構(gòu)特征與TI 更相近，說明IVAE重建圖像的質(zhì)量最好。

圖5 TI與SNESIM、VAE、DS、IVAE重建圖像的變差函數(shù)Fig.5 Variograms of TI and reconstructed images by using SNESIM，VAE，DS and IVAE

3.2.3 孔隙分布比較

通過對TI 和SNESIM、VAE、DS 以及IVAE 的重建圖像的孔隙進行分析，可以得到重建圖像內(nèi)部孔隙的數(shù)量、大小和孔隙直徑的分布情況。孔徑近似定義為：

其中V為孔隙的體積。表2表示TI和30次重建圖像的平均孔隙數(shù)目，表3 表示TI 和重建圖像中的孔徑。如表2 和表3 所示，IVAE與TI的孔隙數(shù)和孔徑更為接近。

表2 TI和每種方法的30張重建圖像的平均孔隙數(shù)Tab.2 Average numbers of pores in TI and 30 reconstructed images by using each method

表3 TI和每種重建方法重建圖像中的孔徑Tab.3 Pore diameters in TI and reconstructed images by using each method

圖6 為TI 和每種方法的重建圖像的孔隙分布情況，由于孔徑分布圖像不能為TI 和重建結(jié)果之間的差異度提供定量測量，因此，在孔徑分布中引入差異度函數(shù)DD，定義為：

圖6 TI和各方法的重建圖像孔隙直徑分布Fig.6 Distribution of pore diameters of TI and reconstructed images by using each method

其中：reSNESIM表示重建圖像方法；Nd和nd分別表示孔徑為d的訓練圖像與重建圖像的孔隙數(shù)，DD越小表示差異度越小，從表4可以看出IVAE 與TI的孔徑分布差異度最小，表明了本方法更具有優(yōu)越性。

表4 四種方法與訓練圖像的孔徑分布差異度Tab.4 Differences of pore diameter distribution between four methods and TI

3.2.4 CPU性能與內(nèi)存比較

本文的實驗環(huán)境：CPU 型號為Inter Core i5，8 GB 內(nèi)存，GPU 型號為Nvidia GeForce GTX970（4 GB）。為了比較每種方法的平均性能，表5 記錄了IVAE、VAE、SNESIM 和DS 方法30次重建時的CPU/GPU 平均利用率和重建時間。表5的重建時間列出了各方法第1次重建時間和剩余29次重建的平均時間。因為IVAE和VAE第1次重建時間包含了模型的訓練時間，經(jīng)過第1次重建時的訓練，IVAE和VAE模型的參數(shù)被保存，之后的每次重建時間將被縮短。而如SNESIM和DS傳統(tǒng)空間數(shù)據(jù)重建方法在每次的訓練過程中都需要重新掃描數(shù)據(jù)，所以將IVAE 和VAE 與傳統(tǒng)空間數(shù)據(jù)重建方法分成第1 次和其余29次時間比較。如表5 所示，IVAE 和VAE 在第1 次重建時間和其他29次重建時間上均具有較大優(yōu)勢，并且均可利用GPU 進行重建計算，而其他經(jīng)典的重建方法只能依靠CPU完成重建，所以在重建時間和CPU負荷方面IVAE和VAE更有優(yōu)勢。

表5 四種方法在30次重建中的平均內(nèi)存消耗、CPU/GPU利用率和運行時間Tab.5 Average memory consumption，CPU/GPU utilization and running time of each method in 30 reconstructions of each method

包括IVAE 的深度學習方法比傳統(tǒng)方法具有的另外一個明顯優(yōu)勢在于：每次進行新的空間數(shù)據(jù)重建時，SNESIM 和DS等傳統(tǒng)重建方法都要重新掃描訓練數(shù)據(jù)來建立一個新的模式庫，但IVAE并不需要重新學習。因為傳統(tǒng)重建方法一般只在內(nèi)存中存儲訓練模式，導致每當重建過程結(jié)束，內(nèi)存中的訓練模式數(shù)據(jù)就被清除；但是IVAE模型經(jīng)過第一輪的訓練后將訓練模型的參數(shù)永久地存儲在硬盤，可以用于未來的重建過程，當需要重建與之前訓練模型同類的空間信息時，可直接將空間信息的結(jié)構(gòu)特征參數(shù)用于重建，大大減少了重建的總時間。

4 結(jié)語

傳統(tǒng)的不確定性空間信息重建方法在每次重建時都需要掃描訓練數(shù)據(jù)以建立訓練模型，因此會花費較多訓練時間，并且會占用大量的CPU 和內(nèi)存資源。本文將費雪信息量和變分自編碼器結(jié)合應用于空間信息不確定性重建。通過編碼器提取輸入數(shù)據(jù)的特征，再建模得到輸入數(shù)據(jù)的概率分布并在編碼器中根據(jù)信息量重新定義證據(jù)下界函數(shù)，然后通過無偏采樣得到與輸入數(shù)據(jù)具有相似分布的數(shù)據(jù)作為解碼器的輸入數(shù)據(jù)，最后的輸出數(shù)據(jù)具有與輸入數(shù)據(jù)相似的結(jié)構(gòu)特征。實驗表明本方法得到的重建結(jié)果具有更高的精度和更優(yōu)的結(jié)構(gòu)特征，在內(nèi)存/CPU 占用率和時間消耗方面也優(yōu)于傳統(tǒng)不確定性空間信息重建方法。