郭 龍，姚樹新，鄭法威

(1.中海油信息科技有限公司，廣東 深圳 518000；2.中國科學(xué)院 深圳先進技術(shù)研究院，廣東 深圳 518000)

多孔材料(例如吸附劑、鋰電池、巖石、建筑材料等)的評價檢測過程中均需要對其進行計算機斷層掃描(Computed Tomography，CT)，從而得到三維多孔介質(zhì)模型，然后對模型進行滲透率分析。這種模型滲透率分析需要計算流體動力學(xué)(Computational Fluid Dynamics，CFD)技術(shù)。流體模擬目前是計算機算法的一個重要研究分支[1]?；诠こ涛锢碇械牧黧w力學(xué)模型，已有眾多計算研究對CFD算法進行了優(yōu)化。滲透率計算的準(zhǔn)確程度通常與流體的細節(jié)相關(guān)，例如在地質(zhì)學(xué)上要計算出更精確的數(shù)字巖心滲透率，即需要更多的網(wǎng)格采樣點。網(wǎng)格點的增多使計算開銷增大，導(dǎo)致模擬速度變慢，同時也引起了內(nèi)存需求增加。

鉆井時，通過巖屑迅速CT掃描工具，可獲得數(shù)字巖心并求取在鉆地層的滲透率。但鉆井速度較快，要得到實時的滲透率信息，仍需要長時間的流體模擬與三維渲染，這使得數(shù)字巖心掃描方法的速度優(yōu)勢已不存在。因此眾多研究者們致力于探尋各種方法，在加快流體模擬速度的同時保證模擬的細節(jié)。通過細節(jié)質(zhì)量確保最終計算結(jié)果準(zhǔn)確與實時的計算要求。

流體模擬的兩大類基本方法為：歐拉網(wǎng)格法與拉格朗日粒子法。其中，歐拉網(wǎng)格法在空間中設(shè)置固定的網(wǎng)格采樣點，通過求解流體方程得到短時間內(nèi)網(wǎng)格上速度、壓強、密度[2]。文獻[3]用格子玻爾茲曼方法(Lattice Boltzmann Method，LBM)得到了Nnvier-Stokes(NS)方程的數(shù)值解，并用該方法進行三維多孔介質(zhì)內(nèi)滲透率的計算。拉格朗日粒子法處理流體為小顆粒組成的粒子系統(tǒng)，使每個顆粒與周圍顆粒互相作用，以求解短時間內(nèi)粒子的速度、加速度、位置等變化[4]。光滑粒子動力學(xué)(Smoothed Particle Hydrodynamics，SPH)[5-6]是一種經(jīng)典的粒子法。在滲透率計算中，粒子法的時間復(fù)雜度隨粒子數(shù)以幾何倍數(shù)增加，且精度較低；而歐拉網(wǎng)格法的數(shù)值精度高且便于實現(xiàn)，是較為常用的方法。

歐拉網(wǎng)格的LBM求解器可分為：碰撞、遷移、計算宏觀量、邊界處理共4個步驟。網(wǎng)格等同于流場的采樣點，也將壓縮感知中的采樣點與流場格子相對應(yīng)。滲透率計算需要統(tǒng)計每個格點的瞬時流量，為減小誤差，通常采用如下3種方法在計算量與精度之間平衡：(1)自適應(yīng)網(wǎng)格；(2)額外細節(jié)；(3)模型縮減。

歐拉方法通常使用自適應(yīng)網(wǎng)格節(jié)省計算資源。文獻[7]使用了八叉樹網(wǎng)格優(yōu)化采樣格點。文獻[8]用粗網(wǎng)格投影，先計算粗網(wǎng)格，再使內(nèi)部細網(wǎng)格滿足無散條件。歐拉網(wǎng)格也可以在粗網(wǎng)格模擬的基礎(chǔ)上，添加湍流細節(jié)，例如文獻[9]在粗網(wǎng)格上的工作。歐拉網(wǎng)格的模型縮減法通過數(shù)據(jù)降維，降低計算的復(fù)雜度。文獻[10]在2006年首次使用主成分分析了將滿足不可壓縮條件的高精度速度場進行降維，獲得代表流場特征的正交基，然后在低維模型中求解方程。結(jié)果表明，其模擬效率大幅提高，但精度也相應(yīng)地顯著降低。文獻[11]在其基礎(chǔ)上將模擬場景進行分塊降維再重新拼接，獲得了在超大尺度空間實時模擬的效率。以上優(yōu)化方法均可減少歐拉方法的計算量，但基本都是以精度大幅下降為代價。傳統(tǒng)網(wǎng)格取樣受限于經(jīng)典的Nyquist采樣理論，在計算時無法將采樣數(shù)量降至定理指出的最小采樣數(shù)。因此，需要借助壓縮感知(Compressive Sensing，CS)突破Nyquist采樣理論的限制，從而降低粗網(wǎng)格計算量，然后再恢復(fù)到高精度計算結(jié)果。由文獻[12～14]建立起的壓縮感知理論，顛覆了Nyquist經(jīng)典采樣定理。其建立了數(shù)學(xué)理論框架，在滿足一定的限制條件下，可以通過一些重構(gòu)優(yōu)化算法從少量的采樣數(shù)據(jù)中較好地還原出原始數(shù)據(jù)。

1 滲透率計算步驟

1.1 使用LBM求解NS方程組

Navier-Stokes方程組是一組描述粘性不可壓縮流體運動的方程。用于生成速度場的NS方程組為

(1)

式中，ρ是流體的密度；u是流體的速度；t是時間。網(wǎng)格中每一個采樣點上均有一個帶方向的速度矢量，組成了一個速度場。三維空間的每個速度矢量有3個分量，分別展開成ux、uy、uz標(biāo)量場，后續(xù)壓縮感知可處理的數(shù)據(jù)即是這種標(biāo)量場切片。g是外力對格點內(nèi)密度造成的加速度，v是流體的粘滯系數(shù)。在統(tǒng)計物理學(xué)中，式(1)等價于Boltzmann方程的粒子速度分布描述函數(shù)。文獻[15]指出了格子玻爾茲曼方程是玻爾茲曼方程的有限差分求解形式，即LBM可以通過差分格式求解NS方程。Boltzmann-BGK方程對粒子速度分布函數(shù)f時空變化進行了描述[16]

(2)

式中，f=f(x，c，t)是粒子速度分布函數(shù)；c為格子上的離散速度；g為外力及表面力產(chǎn)生的加速度；f對c的偏導(dǎo)數(shù)是速度空間中的單方向?qū)?shù)；τ為上述單一松弛時間模型中的特性松弛時間。在方程右側(cè)是BGK碰撞模型，其中feq是無量綱的離散速度空間的局部平衡態(tài)分布函數(shù)。選擇合適的f是構(gòu)造格子Boltzmann模型的關(guān)鍵，而f的形式與構(gòu)造離散速度模型相關(guān)。文獻[15]提出的DnQm系列模型是目前被廣泛應(yīng)用的模型，其中n表示維數(shù)，m是離散速度c的方向個數(shù)。D2Q7、D2Q9、D3Q15模型，如圖1所示。

圖1 DnQm系列模型的離散速度方向

這類模型的f可以寫成以下形式

(3)

式中，cs是格子的聲速，也是一個常數(shù)；ρ是當(dāng)前格子的流體密度；u是瞬時的流體速度；ci是本格子第i個方向的離散速度；wa是權(quán)系數(shù)。通過替換等式中的平衡速度ueq來實現(xiàn)外力(或壓力梯度)作用。

(4)

其中，F(xiàn)為外力項，其在速度u的坐標(biāo)軸x、y、z方向上分解后，可以以加速度g=τF/ρ的形式與速度分量相加得到平衡速度。由外力作用后的f(t+Δt)，以式(5)的形式傳播迭代。從而完成一次“碰撞、遷移、計算宏觀量、邊界處理”計算循環(huán)。

(5)

式中，Δt為一次迭代所消耗的格子時間，其需要經(jīng)過格子尺度轉(zhuǎn)化為真實世界時間，具體模擬方法已較為成熟，本文不再贅述。在同樣的真實世界尺度下，不同精度的網(wǎng)格計算滲透率具有差異。如圖2所示，在100×41×41的管道模型中，格子數(shù)量降低至10×5×5。則本該是近似拋物面包圍面積的區(qū)域被線性插值包圍面積所替代，從而造成流場計算誤差，影響滲透率的計算。

根據(jù)滲透率的定義，速度場在某個方向上的投影切片均為該切片上的滲透率。本文中的流場是指流體速度場。如何在較小網(wǎng)格精度下得到高精度的速度場，從而得到高精度的滲透率計算結(jié)果，成為了流場壓縮感知的重要問題。

1.2 壓縮感知

壓縮感知理論[14]可推導(dǎo)出：只要流場截面中包含眾多零或通過變換得到一個稀疏域，即可將采樣頻率降至低于Nyquist采樣定理要求的頻率，然后通過重構(gòu)從稀疏流場中恢復(fù)高精度流場。壓縮感在大量減少流場采樣的同時，保證將滲透率較好地計算還原出來。

定義一個n維信號向量r∈Rn，當(dāng)r中只有k個非零采樣點而其余均為0時，若k?n，則稱r是稀疏的。向量r的稀疏度為k/n，而r的l0范數(shù)含義是其中非零元素的個數(shù)。

‖r‖0=k

(6)

通常流場在x、y、z方向上的分量均不是稀疏的，則需要尋找一種正交稀疏變換方法Ψ將r轉(zhuǎn)換為變換域下的向量s，使得s是稀疏的。常用的稀疏變換方法為離散余弦、小波變換等。變換如式(7)所示。

r=ψs

(7)

式中，Ψ為壓縮基，是一個n×n的正交矩陣。再定義一個m維向量b∈Rm，其中m?n。b是由r通過矩陣M下采樣得到的，M為采樣基。m×n的感知矩陣A是采樣基與壓縮基的乘積，即整個壓縮感知可用下式表達

b=Mr=Mψs=As

(8)

其目的是從流場截面b中恢復(fù)出r。需要求解的精確流場切片即為式(8)解集中最稀疏的解。于是問題轉(zhuǎn)換為求解一個l0范數(shù)最小的優(yōu)化問題

min(‖s‖0) s.t.b=As

(9)

但求解l0范數(shù)最小問題是NP難問題。因此，需要放寬優(yōu)化條件為求l1范數(shù)，其含義是向量中各元素的絕對值之和，表達式為

min(‖s‖1) s.t.b=As

(10)

式(10)是壓縮感知的基追蹤形式。根據(jù)式(7)，若希望恢復(fù)的是一個N×N的二維高精度流場，則轉(zhuǎn)換成向量r后是一個N2×1的向量，稀疏向量也是一個N2×1的向量s。此時，壓縮基Ψ的大小為N2×N2。該矩陣相對較大，計算矩陣乘法時的速度也較慢。因此使用的壓縮矩陣是直接與N×N的流場矩陣相乘，對其進行直接變換，不需要將流場數(shù)據(jù)轉(zhuǎn)換為一維向量。

2 速度場稀疏重建

2.1 稀疏變換方法

考慮到流場在某方向上的分量切片與普通灰度圖像的相似性，可以用灰度圖像的壓縮感知理論來處理流場信息。對于較平滑的斑點圖像，主要有兩種稀疏變換方法：離散小波變換和離散曲波變換。二維離散小波變換每層變換均會分解為4個大小相同的系數(shù)矩陣，其長寬均是被分解數(shù)據(jù)的1/2。若邊長出現(xiàn)非2的冪次數(shù)值，則在分解前補0。圖3是對一個真實巖心圖像計算所得流場切片的Haar小波變換示意圖。離散小波變換后的系數(shù)矩陣數(shù)據(jù)量大幅減少，只需用大約10 %的數(shù)據(jù)即可高概率的還原至原速度場切片。

(a)

壓縮感知中稀疏度越高，則可在相同數(shù)據(jù)量下恢復(fù)出精度更高的流場信息。本文對比了Haar、Daubechies、Symlets Biorthogonal共3種小波在多層分解條件下的稀疏度，如表1所示。

表1 流場切片在多種子波壓縮感知中的稀疏度

表1中的數(shù)據(jù)顯示，使用Symlets Biorthogonal小波比Haar與Daubechies變換后的矩陣稀疏度更小?？紤]到小波基較多，限于篇幅，本文使僅對比幾類常用的小波基。下文中所對比使用的小波變換默認采用Symlets Biorthogonal小波基。

曲波變換(Curvelet)是基于傅里葉變換與小波變換的高度各向異性的一種改進，常用于恢復(fù)圖形邊緣和抑制周邊噪聲。曲波變換各向異性的優(yōu)勢在于拉伸、平移的基礎(chǔ)上同時引入了一個旋轉(zhuǎn)變化。本文用MATLAB軟件的CurveLab-2.0程序包實現(xiàn)了Wrapping形式的曲波變換，并將曲波變換后得到的系數(shù)矩陣稀疏度與小波變換進行比較，如表2所示。通常其在進行滲透率計算時，總是向特定方向上施加一個壓力梯度以驅(qū)動流體流動，流場在壓力梯度方向的投影絕對值明顯大于其他方向，造成了各向異性。而二維曲波變換在各向異性圖像中生成的稀疏矩陣系數(shù)更少，稀疏程度更高，更適合作為流場切片壓縮感知的壓縮基。

表2 小波變換和曲波變換稀疏度對比

2.2 采樣基

通常LBM模擬框架中，使用的是規(guī)則網(wǎng)格?？s減網(wǎng)格相當(dāng)于將高精度速度場r變成低精度流場b的過程。一般低精度流場也是規(guī)則網(wǎng)格，若隨機采樣，則采樣點是不規(guī)則網(wǎng)格難以匹配，因此下文使用均勻下采樣。由于孔隙分布是隨機的，這種采樣滿足隨機采樣標(biāo)準(zhǔn)。構(gòu)造的采樣矩陣，將網(wǎng)格精度縮減為原來的0.5倍。對于一個N×N的二維高精度流場，其某一方向上的切片是一個N2×1的向量，得到的低精度流場大小為(N/2)×(N/2)，采樣信號b是一個N2/4×1的向量。此時采樣矩陣M的大小為N2/4×N2，其為一個較大的矩陣。這里將采樣矩陣拆分成兩個矩陣，在采樣時不需要將速度場切片轉(zhuǎn)換為一維向量，而是與兩個N×(N/2)的矩陣相乘。該采樣矩陣的構(gòu)造原理如下：當(dāng)網(wǎng)格數(shù)下降至0.5倍原始網(wǎng)格時，需要將2×2的4個數(shù)據(jù)合成1個位于中心的數(shù)據(jù)，每個數(shù)據(jù)雙線性插值的系數(shù)均為1/2；采樣矩陣先左乘流場矩陣得到(N/2)×N的中間結(jié)果，相當(dāng)于對流場切面矩陣的縱向進行單次線性插值，然后用中間結(jié)果右乘采樣矩陣的轉(zhuǎn)置，相當(dāng)于再進行了一次橫向的線性插值。

(11)

(12)

2.3 流場重構(gòu)算法

重構(gòu)算法的主要目的是求解式(10)，即稀疏優(yōu)化問題。為了解決稀疏優(yōu)化問題，本文對比了幾個高引用率的算法：梯度投影稀疏重建(Gradient Projection for Sparse Reconstruction，GPSR)[17]、譜投影梯度L1最小化(Spectral Projected-Gradient forL1minimization，SPGL1)[18]、快速迭代收縮閾值算法(Fast Iterative Shrinkage Thresholding Algorithm，F(xiàn)ISTA)[19]與凸集投影(Projection onto Convex Set，POCS)[20]。在均勻取樣的條件下，本文比較了這些取樣和稀疏優(yōu)化的算法。首先進行高精度網(wǎng)格的LBM模擬得到高精度速度場，隨機選擇其中一個X分量的256×256格子的切片1。然后進行0.5倍精度的網(wǎng)格模擬，得到低精度速度場。如圖4所示，在低精度速度場對應(yīng)切片圖4(a)，用壓縮感知稀疏重建方法恢復(fù)速度場得到切片圖4(b)。對比切片1與切片2得到兩者的均方誤差(Mean-Square Error，MSE)，即可得出不同稀疏恢復(fù)方法在流場壓縮感知中的適應(yīng)性，如表3所示。

(a)

表3 小波變換和曲波變換稀疏度對比

從表3的MSE列可以看出，4種稀疏優(yōu)化方法求解而恢復(fù)出的高精度流場切片與直接用高精度網(wǎng)格計算的結(jié)果較相似。原高精度網(wǎng)格計算的滲透率為16.1 md，低精度網(wǎng)格計算的滲透率為15.4 md。4種方法恢復(fù)出的切面滲透率與原始高精度流場計算的滲透率最大差0.2 md，從低精度網(wǎng)格恢復(fù)的滲透率與高精度網(wǎng)格模擬的結(jié)果較為接近。

表4 恢復(fù)重構(gòu)時間對比

作者使用Intel I7-8700K CPU @ 4.7 GHz與雙通道16 GB@3.8 GHz內(nèi)存多線程運行MATLAB，并得到了表4的運行時間數(shù)據(jù)。切片大小設(shè)置為256×256網(wǎng)格，計算流場共進行了2 000次迭代達到平衡態(tài)，用時275 ms。若使用1 024×1 024網(wǎng)格計算，則用時2 140 ms。使用高精度的流場計算時間減去低精度計算時間，再減去重建恢復(fù)時間即可得到算法節(jié)省的計算時間。在4種稀疏優(yōu)化方法中，GPSR方法和POCS恢復(fù)重建速度較快。

3 結(jié)束語

本文在多孔介質(zhì)的滲透率計算LBM模擬中引入了壓縮感知，利用切片上1/4的流場采樣點，以較低的MSE恢復(fù)了高精度的流場信息?；趬嚎s感知理論，構(gòu)建了滲透率計算中流體模擬壓縮感知采樣方法與稀疏恢復(fù)方法。

根據(jù)計算網(wǎng)格的特性，本文使用均勻采樣矩陣作為壓縮感知的采樣基。對比針對l1范數(shù)最小的4種重構(gòu)算法，且對比了曲波變換、小波變換的多種壓縮基的稀疏度。最后，利用實際巖心的256×256切片模擬實際的滲透率計算。從模擬結(jié)果中可以看出，壓縮感知在滲透率計算中的應(yīng)用能夠?qū)⒌途W(wǎng)格數(shù)的速度場切片計算結(jié)果恢復(fù)到高精度的速度場切片結(jié)果，從而提高計算效率。