胡 進, 查文舒, 劉鎮(zhèn)領(lǐng), 李道倫
(1.合肥工業(yè)大學(xué) 數(shù)學(xué)學(xué)院,安徽 合肥 230601; 2.渤海鉆探油氣井測試分公司, 河北 廊坊 065007)
重構(gòu)數(shù)字巖心的方法[1-4]主要分為物理實驗方法和數(shù)值重建方法兩大類。物理實驗方法借助高精度儀器獲取巖心的平面圖像,然后對平面圖像進行重建即可得到數(shù)字巖心。物理重建的常用方法有:序列成像法[5]、核磁共振成像法[6]、X射線計算機層析成像(X-CT)掃描法[7]和聚焦離子束電子顯微鏡(FIB-SEM)掃描法[8]等。這些方法能夠?qū)r心樣本的真實孔隙結(jié)構(gòu)進行定性描述, 但視野小,而且花費代價大。為便于分析與獲得更大尺度孔隙結(jié)構(gòu),基于數(shù)字巖心重建與重構(gòu)方法得到了重視。
數(shù)值重建方法是利用得到的電鏡掃描圖片, 結(jié)合不同的數(shù)學(xué)統(tǒng)計方法對巖樣進行重構(gòu)。這類方法統(tǒng)計參數(shù)的取得方式比較靈活, 花費的成本比較低,不僅可以大范圍地進行重構(gòu),而且可以滿足重構(gòu)對精度的要求。基于數(shù)值的重構(gòu)方法主要有過程法[9]、順序指示模擬技術(shù)[10]、模擬退火法[11]、多點地質(zhì)統(tǒng)計法[12]、馬爾科夫鏈蒙特卡羅方法[13]等。
多點地質(zhì)統(tǒng)計法一直是數(shù)值重建研究的熱點,它是相對于兩點地質(zhì)統(tǒng)計法[14]而言的。傳統(tǒng)的兩點地質(zhì)統(tǒng)計方法利用變差函數(shù)來研究地質(zhì)變量空間的相關(guān)性,主要把握空間上兩點之間的相關(guān)性,難以表征復(fù)雜的空間結(jié)構(gòu)和再現(xiàn)復(fù)雜目標的幾何形態(tài)。針對這種方法的不足,多點地質(zhì)統(tǒng)計學(xué)方法應(yīng)運而生。在多點地質(zhì)統(tǒng)計學(xué)中, 應(yīng)用訓(xùn)練圖像[15]代替變差函數(shù)表達地質(zhì)變量的空間結(jié)構(gòu)性, 這使得多點地質(zhì)統(tǒng)計學(xué)建模方法與計算機圖像紋理合成領(lǐng)域展示出目的的一致性,即都是隨機生成與訓(xùn)練圖像擁有相似特征的圖像[16]。
多點地質(zhì)統(tǒng)計學(xué)方法與計算機圖形學(xué)算法的融合逐步得到了學(xué)者重視。文獻[17]首次將計算機圖形學(xué)算法中的圖像縫合算法用于多點地質(zhì)統(tǒng)計的條件模式模擬,并命名為CIQ算法。CIQ方法與多點地質(zhì)統(tǒng)計學(xué)模式模擬方法直接用新模擬區(qū)域的值更新重疊區(qū)域值的方法不同,它是將重疊區(qū)域沿著重合區(qū)域誤差最小處進行縫合??p合方法的使用緩解了原先直接取代方法造成的重疊區(qū)域邊界突出的缺陷[18],得到了更好的模擬結(jié)果。
模擬結(jié)果的對比如圖1所示。這種方法適用于規(guī)則矩形的重合區(qū)域。
圖1 重疊區(qū)域直接覆蓋和最小誤差拼接的模式模擬效果示意圖
文獻[19]將計算機圖切算法應(yīng)用于地質(zhì)巖相反演模擬。通過使用馬爾科夫鏈融合后驗信息,圖形切割技術(shù)能夠從訓(xùn)練圖像中找到合適的模式加入到模擬結(jié)構(gòu)中。
CIQ算法和計算機圖切算法都成功應(yīng)用在地質(zhì)模擬中,對于特征信息量較小的河道圖像,重構(gòu)效果好。但是,對于本文中有著大量特征信息的巖心圖像,在用上述算法對其重構(gòu)時,會出現(xiàn)部分特征的缺失和重復(fù)等不足。本文在文獻[17,19]基礎(chǔ)上,針對大樣本的頁巖巖心數(shù)字圖像,提出了一種基于多樣子圖的圖像紋理重構(gòu)算法。這種方法從訓(xùn)練圖像中挑選出特征子圖,可以避免合成圖中的重要信息缺失;同時設(shè)置比例,也可避免合成圖中的大面積的重復(fù);另外每次都是從挑選的特征子圖中選擇匹配塊,避免了全局搜索,大大提高了合成速度。
頁巖氣儲層孔隙結(jié)構(gòu)特征對頁巖氣勘探開發(fā)具有重要意義[20]。通過掃描電子顯微鏡可以對這些微裂隙的微觀結(jié)構(gòu)進行更直觀地觀察,獲得孔隙的大小、形態(tài)及有機質(zhì)分布等,從而為研究頁巖氣的成藏機理、運移規(guī)律及開發(fā)提供基礎(chǔ)數(shù)據(jù)。本文用頁巖巖心的電鏡掃描數(shù)字圖像作為原始訓(xùn)練圖像,如圖2所示。為去除噪聲,對圖2進行三值化處理[21],通過設(shè)定合適的閾值,將圖像上像素點的灰度值變?yōu)?、128或255,其中0表示有機質(zhì),255表示白色礦石結(jié)果,如圖3所示。
圖2 原始訓(xùn)練圖像
為了避免合成結(jié)果中出現(xiàn)大面積的特征信息重復(fù)或缺失,本文方法不直接從圖3中挑選匹配塊,而是首先人工選取m幅的特征子圖(本文中m取11),如圖4所示。這些特征子圖包含著原始頁巖巖心圖像中的礦石、裂縫、有機質(zhì)等重要信息。然后給特征子圖預(yù)設(shè)定在合成中的比例,依次為a1,a2,…,am。通過設(shè)置不同的ai值(i=1,2,…,m),可以控制各組分在合成圖像中的比例,達到比例可調(diào)控的目的。
圖4 11幅特征子圖
本文算法的合成過程與文獻[22-23]類似, 每次從訓(xùn)練圖像中挑選ε個最相似的匹配塊隨機選擇一個進行匹配,然后用Graph Cut算法決定哪部分將輸出到合成圖中, 這樣不斷地擴充合成區(qū)域直至完成整幅圖像。
不同的是,本文中的訓(xùn)練圖像是從圖3中挑選出的多幅特征子圖,即有多個訓(xùn)練圖像,因此在匹配的過程中第1步要解決的問題就是訓(xùn)練圖像的選擇。
(1)
bm表示在合成過程中每幅子圖所占當前合成區(qū)域的比例與其預(yù)設(shè)比例的差值。在合成過程中盡量要求每幅子圖的Nm/N接近它的預(yù)設(shè)比例am,因此每次選擇最小的bm值對應(yīng)的特征子圖作為下一個訓(xùn)練圖像,這樣就能保證在最后的合成圖像中每幅特征子圖所占比例與預(yù)設(shè)比例大致相同。
重疊區(qū)域的拼接也是影響重構(gòu)結(jié)果的關(guān)鍵因素。本文沿用計算機圖切技術(shù)中的Graph Cut 能量最小化算法完成重疊區(qū)域的拼接,它最早是由 文獻[24]提出的。
本文通過復(fù)制訓(xùn)練圖像中的矩形像素塊到輸出圖像中合成紋理,再通過圖切算法在矩形紋理塊中選出最優(yōu)的不規(guī)則區(qū)域用于合成。
在Graph Cut算法執(zhí)行過程中,最佳路徑由許多相關(guān)的像素點組合而成,像素點的生成是從每一對像素點中選擇而來,最簡單的匹配標準是比較相鄰像素點的顏色誤差。設(shè)t、u為重疊區(qū)域中的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)
這樣的邊值對應(yīng)著重疊區(qū)域的模式差異值大小,使用基于最小能量的圖像剪切技術(shù)可以將2個模式以最小誤差的方式拼接起來,過程如圖5所示。圖5中的左邊為重疊區(qū)域示意圖,右邊紅線代表接縫。
圖5 2個模式以最小誤差方式的拼接過程
(1) 根據(jù)實驗需要,對原始的數(shù)字巖心圖像做三值化處理,提取圖中的礦石、有機質(zhì)、裂縫等重要信息。
(2) 從處理后的三值化圖像中選取m幅特征子圖作為訓(xùn)練圖像。
(3) 輸入m幅訓(xùn)練圖像、初始參數(shù)(匹配塊大小p、重疊區(qū)域大小O和候選匹配塊數(shù)量ε等)。
(4) 計算bm的值,從對應(yīng)最小的bm值的訓(xùn)練圖像中隨機選取一個匹配塊開始重構(gòu)。
(5) 迭代直到整個圖像重構(gòu)完成。步驟如下:① 根據(jù)重疊區(qū)域的相似性大小,在對應(yīng)bm值最小的訓(xùn)練圖像中搜索找到下一個匹配塊;② 從ε個最相似的匹配塊中隨機選擇一個輸入到待合成區(qū)域中;③ 在待合成區(qū)域中使用Graph Cut技術(shù)將待合成圖像與匹配塊的重疊區(qū)域進行誤差最小化拼接。
三值化后的頁巖巖心數(shù)字圖像(圖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。初始參數(shù)p、O和ε的設(shè)定對于本文的重構(gòu)結(jié)果也具有一定的影響[25]。從經(jīng)驗上講,重構(gòu)過程中匹配塊大小為原始訓(xùn)練圖像的1/4~1/6之間為最好,故本文匹配塊大小p=200×200。對于重疊區(qū)域,如果使用較小的重疊區(qū)域,無法保持空間結(jié)構(gòu)的連續(xù)性,但是使用過大的重疊區(qū)域,會導(dǎo)致重構(gòu)時間的大大增加,故重疊區(qū)域大小設(shè)置為匹配塊的1/6為宜(O=35×200)。而對于候選匹配塊數(shù)量ε,取值太小容易出現(xiàn)大面積的訓(xùn)練圖像復(fù)制,增大候選匹配塊數(shù)目減少了對訓(xùn)練圖像的復(fù)制,增加了匹配的多樣性,故本文中ε=10。
本文算法實現(xiàn)的部分重構(gòu)結(jié)果如圖6所示,大小為1 000×1 000。CIQ算法在相同的初始參數(shù)條件下實現(xiàn)的部分重構(gòu)結(jié)果如圖7所示,大小為1 000×1 000。
圖6 本文算法重構(gòu)結(jié)果
圖7 CIQ方法重構(gòu)結(jié)果
從實驗結(jié)果上看,CIQ方法重構(gòu)出的圖像出現(xiàn)了一定面積的重復(fù)(圖7中紅色標記區(qū)域),而且缺失了很多重要的特征,而本文重構(gòu)的結(jié)果基本保持了原始訓(xùn)練圖像具有的特征,還因為引入了比例控制,各特征組分比例也基本相同,同時每次搜索特征子圖隨機選擇匹配塊,搜索范圍變小,整個重構(gòu)所需時間更少。
本文算法引進了控制比例的方法,下面來驗證算法的可調(diào)控性。
重構(gòu)圖像由36塊紋理塊合成,算法中的比例參數(shù)為am,以第6幅特征子圖為例,在本文算法中a6設(shè)置為0.025,因此在重構(gòu)結(jié)果中該特征只出現(xiàn)了一次。當將a6設(shè)置成0.05時,實驗結(jié)果如圖8所示。
圖8 當a6=0.05時的重構(gòu)結(jié)果
從圖8可以看出,當a6設(shè)置成原來的2倍時,該特征子圖在重構(gòu)結(jié)果中出現(xiàn)了2次。接下來將a6設(shè)置成0.075,重構(gòu)結(jié)果如圖9所示。
從圖9可以看出,當a6設(shè)置成原來的3倍時,該部分特征在重構(gòu)結(jié)果中出現(xiàn)了3次,且重構(gòu)后的圖像仍然保持著原圖的其他特征。由此可以看出本文算法在基于控制比例的條件下,可以重構(gòu)出較好的結(jié)果。
圖9 當a6=0.075時的重構(gòu)結(jié)果
本文提出了基于圖像紋理合成的特征可調(diào)控的數(shù)字巖心重構(gòu)方法,具有以下特點:
(1) 對大樣本的頁巖巖心圖像重構(gòu)效果好,能基本保持原始訓(xùn)練圖像的特征,各組分比例也大致相同。
(2) 快速?,F(xiàn)有的實時或近實時圖像重構(gòu)算法需要從全局圖像搜索匹配塊,而本文算法在其特征子圖中選取匹配塊,避免了在原始訓(xùn)練圖像中全局搜索,提高了匹配的速度。
(3) 較少的參數(shù)設(shè)置。本文在匹配的時候設(shè)置了一個比例挑選條件,無需設(shè)置很多的參數(shù),就能夠?qū)崿F(xiàn)對特征比例的控制。