流場壓縮感知滲透率計算

郭 龍，姚樹新，鄭法威

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

多孔材料(例如吸附劑、鋰電池、巖石、建筑材料等)的評價檢測過程中均需要對其進行計算機斷層掃描(Computed Tomography，CT)，從而得到三維多孔介質模型，然后對模型進行滲透率分析。這種模型滲透率分析需要計算流體動力學(Computational Fluid Dynamics，CFD)技術。流體模擬目前是計算機算法的一個重要研究分支[1]。基于工程物理中的流體力學模型，已有眾多計算研究對CFD算法進行了優化。滲透率計算的準確程度通常與流體的細節相關，例如在地質學上要計算出更精確的數字巖心滲透率，即需要更多的網格采樣點。網格點的增多使計算開銷增大，導致模擬速度變慢，同時也引起了內存需求增加。

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

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

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

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

1 滲透率計算步驟

1.1 使用LBM求解NS方程組

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

(1)

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

(2)

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

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

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

(3)

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

(4)

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

(5)

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

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

1.2 壓縮感知

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

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

‖r‖0=k

(6)

通常流場在x、y、z方向上的分量均不是稀疏的，則需要尋找一種正交稀疏變換方法Ψ將r轉換為變換域下的向量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中恢復出r。需要求解的精確流場切片即為式(8)解集中最稀疏的解。于是問題轉換為求解一個l0范數最小的優化問題

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

(9)

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

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

(10)

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

2 速度場稀疏重建

2.1 稀疏變換方法

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

(a)

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

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

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

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

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

2.2 采樣基

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

(11)

(12)

2.3 流場重構算法

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

(a)

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

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

表4 恢復重構時間對比

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

3 結束語

本文在多孔介質的滲透率計算LBM模擬中引入了壓縮感知，利用切片上1/4的流場采樣點，以較低的MSE恢復了高精度的流場信息。基于壓縮感知理論，構建了滲透率計算中流體模擬壓縮感知采樣方法與稀疏恢復方法。

根據計算網格的特性，本文使用均勻采樣矩陣作為壓縮感知的采樣基。對比針對l1范數最小的4種重構算法，且對比了曲波變換、小波變換的多種壓縮基的稀疏度。最后，利用實際巖心的256×256切片模擬實際的滲透率計算。從模擬結果中可以看出，壓縮感知在滲透率計算中的應用能夠將低網格數的速度場切片計算結果恢復到高精度的速度場切片結果，從而提高計算效率。