胡 進, 查文舒, 劉鎮領, 李道倫
(1.合肥工業大學 數學學院,安徽 合肥 230601; 2.渤海鉆探油氣井測試分公司, 河北 廊坊 065007)
重構數字巖心的方法[1-4]主要分為物理實驗方法和數值重建方法兩大類。物理實驗方法借助高精度儀器獲取巖心的平面圖像,然后對平面圖像進行重建即可得到數字巖心。物理重建的常用方法有:序列成像法[5]、核磁共振成像法[6]、X射線計算機層析成像(X-CT)掃描法[7]和聚焦離子束電子顯微鏡(FIB-SEM)掃描法[8]等。這些方法能夠對巖心樣本的真實孔隙結構進行定性描述, 但視野小,而且花費代價大。為便于分析與獲得更大尺度孔隙結構,基于數字巖心重建與重構方法得到了重視。
數值重建方法是利用得到的電鏡掃描圖片, 結合不同的數學統計方法對巖樣進行重構。這類方法統計參數的取得方式比較靈活, 花費的成本比較低,不僅可以大范圍地進行重構,而且可以滿足重構對精度的要求。基于數值的重構方法主要有過程法[9]、順序指示模擬技術[10]、模擬退火法[11]、多點地質統計法[12]、馬爾科夫鏈蒙特卡羅方法[13]等。
多點地質統計法一直是數值重建研究的熱點,它是相對于兩點地質統計法[14]而言的。傳統的兩點地質統計方法利用變差函數來研究地質變量空間的相關性,主要把握空間上兩點之間的相關性,難以表征復雜的空間結構和再現復雜目標的幾何形態。針對這種方法的不足,多點地質統計學方法應運而生。在多點地質統計學中, 應用訓練圖像[15]代替變差函數表達地質變量的空間結構性, 這使得多點地質統計學建模方法與計算機圖像紋理合成領域展示出目的的一致性,即都是隨機生成與訓練圖像擁有相似特征的圖像[16]。
多點地質統計學方法與計算機圖形學算法的融合逐步得到了學者重視。文獻[17]首次將計算機圖形學算法中的圖像縫合算法用于多點地質統計的條件模式模擬,并命名為CIQ算法。CIQ方法與多點地質統計學模式模擬方法直接用新模擬區域的值更新重疊區域值的方法不同,它是將重疊區域沿著重合區域誤差最小處進行縫合。縫合方法的使用緩解了原先直接取代方法造成的重疊區域邊界突出的缺陷[18],得到了更好的模擬結果。
模擬結果的對比如圖1所示。這種方法適用于規則矩形的重合區域。

圖1 重疊區域直接覆蓋和最小誤差拼接的模式模擬效果示意圖
文獻[19]將計算機圖切算法應用于地質巖相反演模擬。通過使用馬爾科夫鏈融合后驗信息,圖形切割技術能夠從訓練圖像中找到合適的模式加入到模擬結構中。
CIQ算法和計算機圖切算法都成功應用在地質模擬中,對于特征信息量較小的河道圖像,重構效果好。但是,對于本文中有著大量特征信息的巖心圖像,在用上述算法對其重構時,會出現部分特征的缺失和重復等不足。本文在文獻[17,19]基礎上,針對大樣本的頁巖巖心數字圖像,提出了一種基于多樣子圖的圖像紋理重構算法。這種方法從訓練圖像中挑選出特征子圖,可以避免合成圖中的重要信息缺失;同時設置比例,也可避免合成圖中的大面積的重復;另外每次都是從挑選的特征子圖中選擇匹配塊,避免了全局搜索,大大提高了合成速度。
頁巖氣儲層孔隙結構特征對頁巖氣勘探開發具有重要意義[20]。通過掃描電子顯微鏡可以對這些微裂隙的微觀結構進行更直觀地觀察,獲得孔隙的大小、形態及有機質分布等,從而為研究頁巖氣的成藏機理、運移規律及開發提供基礎數據。本文用頁巖巖心的電鏡掃描數字圖像作為原始訓練圖像,如圖2所示。為去除噪聲,對圖2進行三值化處理[21],通過設定合適的閾值,將圖像上像素點的灰度值變為0、128或255,其中0表示有機質,255表示白色礦石結果,如圖3所示。

圖2 原始訓練圖像
為了避免合成結果中出現大面積的特征信息重復或缺失,本文方法不直接從圖3中挑選匹配塊,而是首先人工選取m幅的特征子圖(本文中m取11),如圖4所示。這些特征子圖包含著原始頁巖巖心圖像中的礦石、裂縫、有機質等重要信息。然后給特征子圖預設定在合成中的比例,依次為a1,a2,…,am。通過設置不同的ai值(i=1,2,…,m),可以控制各組分在合成圖像中的比例,達到比例可調控的目的。

圖4 11幅特征子圖
本文算法的合成過程與文獻[22-23]類似, 每次從訓練圖像中挑選ε個最相似的匹配塊隨機選擇一個進行匹配,然后用Graph Cut算法決定哪部分將輸出到合成圖中, 這樣不斷地擴充合成區域直至完成整幅圖像。
不同的是,本文中的訓練圖像是從圖3中挑選出的多幅特征子圖,即有多個訓練圖像,因此在匹配的過程中第1步要解決的問題就是訓練圖像的選擇。

(1)
bm表示在合成過程中每幅子圖所占當前合成區域的比例與其預設比例的差值。在合成過程中盡量要求每幅子圖的Nm/N接近它的預設比例am,因此每次選擇最小的bm值對應的特征子圖作為下一個訓練圖像,這樣就能保證在最后的合成圖像中每幅特征子圖所占比例與預設比例大致相同。
重疊區域的拼接也是影響重構結果的關鍵因素。本文沿用計算機圖切技術中的Graph Cut 能量最小化算法完成重疊區域的拼接,它最早是由 文獻[24]提出的。
本文通過復制訓練圖像中的矩形像素塊到輸出圖像中合成紋理,再通過圖切算法在矩形紋理塊中選出最優的不規則區域用于合成。
在Graph Cut算法執行過程中,最佳路徑由許多相關的像素點組合而成,像素點的生成是從每一對像素點中選擇而來,最簡單的匹配標準是比較相鄰像素點的顏色誤差。設t、u為重疊區域中的2個相鄰像素點,B(t)、C(t)分別為像素點t在舊紋理塊和新紋理塊上的像素顏色值,B(u)、C(u)分別為像素點u在舊紋理塊和新紋理塊上的像素顏色值。連接t、u2個頂點的邊的值定義為:
M(t,u,B,C)=|B(t)-C(t)|+
|B(u)-C(u)|
(2)
這樣的邊值對應著重疊區域的模式差異值大小,使用基于最小能量的圖像剪切技術可以將2個模式以最小誤差的方式拼接起來,過程如圖5所示。圖5中的左邊為重疊區域示意圖,右邊紅線代表接縫。

圖5 2個模式以最小誤差方式的拼接過程
(1) 根據實驗需要,對原始的數字巖心圖像做三值化處理,提取圖中的礦石、有機質、裂縫等重要信息。
(2) 從處理后的三值化圖像中選取m幅特征子圖作為訓練圖像。
(3) 輸入m幅訓練圖像、初始參數(匹配塊大小p、重疊區域大小O和候選匹配塊數量ε等)。
(4) 計算bm的值,從對應最小的bm值的訓練圖像中隨機選取一個匹配塊開始重構。
(5) 迭代直到整個圖像重構完成。步驟如下:① 根據重疊區域的相似性大小,在對應bm值最小的訓練圖像中搜索找到下一個匹配塊;② 從ε個最相似的匹配塊中隨機選擇一個輸入到待合成區域中;③ 在待合成區域中使用Graph Cut技術將待合成圖像與匹配塊的重疊區域進行誤差最小化拼接。
三值化后的頁巖巖心數字圖像(圖3)大小為1 279×887,特征子圖為圖像(圖4)大小依次為209×206、210×202、212×210、221×211、222×211、220×209、200×205、222×212、618×497、562×494、586×468。初始參數p、O和ε的設定對于本文的重構結果也具有一定的影響[25]。從經驗上講,重構過程中匹配塊大小為原始訓練圖像的1/4~1/6之間為最好,故本文匹配塊大小p=200×200。對于重疊區域,如果使用較小的重疊區域,無法保持空間結構的連續性,但是使用過大的重疊區域,會導致重構時間的大大增加,故重疊區域大小設置為匹配塊的1/6為宜(O=35×200)。而對于候選匹配塊數量ε,取值太小容易出現大面積的訓練圖像復制,增大候選匹配塊數目減少了對訓練圖像的復制,增加了匹配的多樣性,故本文中ε=10。
本文算法實現的部分重構結果如圖6所示,大小為1 000×1 000。CIQ算法在相同的初始參數條件下實現的部分重構結果如圖7所示,大小為1 000×1 000。

圖6 本文算法重構結果

圖7 CIQ方法重構結果
從實驗結果上看,CIQ方法重構出的圖像出現了一定面積的重復(圖7中紅色標記區域),而且缺失了很多重要的特征,而本文重構的結果基本保持了原始訓練圖像具有的特征,還因為引入了比例控制,各特征組分比例也基本相同,同時每次搜索特征子圖隨機選擇匹配塊,搜索范圍變小,整個重構所需時間更少。
本文算法引進了控制比例的方法,下面來驗證算法的可調控性。
重構圖像由36塊紋理塊合成,算法中的比例參數為am,以第6幅特征子圖為例,在本文算法中a6設置為0.025,因此在重構結果中該特征只出現了一次。當將a6設置成0.05時,實驗結果如圖8所示。

圖8 當a6=0.05時的重構結果
從圖8可以看出,當a6設置成原來的2倍時,該特征子圖在重構結果中出現了2次。接下來將a6設置成0.075,重構結果如圖9所示。
從圖9可以看出,當a6設置成原來的3倍時,該部分特征在重構結果中出現了3次,且重構后的圖像仍然保持著原圖的其他特征。由此可以看出本文算法在基于控制比例的條件下,可以重構出較好的結果。

圖9 當a6=0.075時的重構結果
本文提出了基于圖像紋理合成的特征可調控的數字巖心重構方法,具有以下特點:
(1) 對大樣本的頁巖巖心圖像重構效果好,能基本保持原始訓練圖像的特征,各組分比例也大致相同。
(2) 快速。現有的實時或近實時圖像重構算法需要從全局圖像搜索匹配塊,而本文算法在其特征子圖中選取匹配塊,避免了在原始訓練圖像中全局搜索,提高了匹配的速度。
(3) 較少的參數設置。本文在匹配的時候設置了一個比例挑選條件,無需設置很多的參數,就能夠實現對特征比例的控制。