張愉玲 , 邢會林 , , 李三忠 , , 逄 碩 , , 劉駿標 , 張熔鑫

(1.深海圈層與地球系統(tǒng)前沿科學中心, 海底科學與探測技術教育部重點實驗室, 中國海洋大學 海洋高等研究院/海洋地球科學學院, 山東 青島 266100; 2.中國海洋大學 海底科學與工程計算國際中心, 山東 青島266100; 3.青島海洋科學與技術國家實驗室 海洋礦產資源評價與探測技術功能實驗室, 山東 青島 266100)

0 引 言

利用CT 掃描等對巖石內部結構進行成像是數(shù)字巖心、礦物乃至底層構造識別的基礎。隨著巖石成像技術的快速發(fā)展, 數(shù)字圖像被廣泛應用。在巖石圖像處理中, 裂隙或者薄構造等重要結構相對于巖體的體積要小得多, 需要高精度與高分辨率的成像數(shù)據。如果對采集的全部數(shù)據進行存儲處理, 不僅占用大量的存儲空間, 而且面對如此龐大的數(shù)據存儲量, 可視化及分析處理也將遇到前所未有的挑戰(zhàn)(翟明國等, 2018)。另外, 圖像特征提取過程中往往會丟失部分重要的特征, 且可能在后續(xù)去噪過程中完全丟失。因此, 需要理清巖石中的主要特征, 并采取有效的圖像處理技術對巖石圖像特征進行精確提取,從而在保留主要特征的前提下降低數(shù)據存儲量。

近年來, 圖像處理技術快速發(fā)展, 蟻群算法(Dorigo et al., 1999; Dorigo and Stüzle, 2003)、神經網絡技術(Zhang et al., 2018)、遺傳算法(Rodrigues et al.,2017)、免疫算法(Yu et al., 2016)等被廣泛運用到計算機圖像處理中。其中, 蟻群算法最早由 Alberto Colorni 在1991 年提出(Colorni and Dorigo et al.,1991), 它作為一種具有正反饋性、魯棒性和并行性的高效率仿生算法, 適用于復雜的圖像數(shù)據處理中。已有學者將蟻群算法運用到圖像特征提取(Ayd?n and U?ur, 2011)、圖像分割(Ma et al., 2009;Arnay et al., 2017)等復雜圖像處理問題。但蟻群算法搜索時間長, 容易停滯而陷入局部最優(yōu), 且修復裂隙特征容易產生多余特征。為了降低數(shù)據復雜度并提高效率, Ester et al. (1996)提出了DBSCAN 聚類算法, 該方法能夠自動分簇, 不需要人工輸入簇數(shù)就可將一組龐大的數(shù)據分為多個簇。此外, 在圖像識別中,Canny 邊緣檢測算子(Canny, 1986)被提出, 并用來準確識別圖像的邊界特征。以上這些方法已經廣泛的運用于人臉圖像識別(孫珊珊, 2016)、醫(yī)學(Yoursefi et al., 2012)、遙感(趙芳等, 2020)和航空航天(陳俠等,2019)等領域, 在地學中, 蟻群算法主要用于斷層識別(葉濤等, 2018), 但未有應用于巖石學的圖像處理中。巖石中的裂隙所占的體積雖小, 但其對于巖石的宏觀力學性質及滲流特征有著重要影響。比如裂隙是含油氣儲層主要的儲集空間和運移通道, 準確提取巖石裂縫特征有助于評估儲層的連通性和滲透性, 進而提高油氣的采收率。本文將綜合利用DBSCAN 聚類算法、螞蟻算法和Canny 算子三種算法, 對巖石精細裂隙特征進行提取、修復和連接, 確保在圖像特征提取過程中保留巖石精細結構的主要幾何形狀特征, 為下一步相關分析計算奠定基礎。

下面簡要介紹蟻群算法和Canny 算子的基本理論, 提出蟻群算法、DBSCAN 聚類算法和Canny 算子的混合算法的框架; 隨后基于PANDAS 軟件平臺(Xing, 2017)對巖石圖像特征提取, 并運用混合算法對提取特征進行修復, 對比單獨使用蟻群算法和混合算法的應用效果; 最終, 獲得較完整的裂隙特征數(shù)據, 并對計算結果進行討論。

1 基于蟻群和Canny 邊緣檢測算子混合算法的圖像處理方案

1.1 蟻群算法原理

20 世紀50 年代, 許多學者開始相繼研究各種仿生算法以尋求最優(yōu)化問題, 其中, 蟻群算法以獨特的優(yōu)勢得到了許多學者的關注。在蟻群算法中, 最經典且與現(xiàn)實生活緊密相關的是旅行商問題(traveling salesmen problem, TSP)。它主要是描述如何找到最短路徑、快速遍歷所有城市。將常見的TSP 數(shù)學問題進行模型的建立, 可以更好理解蟻群算法。

一只螞蟻隨機從n座城市中的某一座城市出發(fā),對所有的城市只訪問一遍絕不重復, 直到遍歷完所有城市螞蟻將停止搜索。設n座城市的集合式中為C={C1,C2, …,Cn}, 兩兩城市之間的距離集合為L={li j|ci,cj?C}, 則G=(C,L)構成有向圖可以直觀地表現(xiàn)出螞蟻遍歷城市的場景, 將現(xiàn)實生活中的問題轉化為對最短Hamilton 圈的求解。城市坐標為(xi,yi)和兩個城市之間的距離滿足下列公式:

設t時刻位于第i座城市上螞蟻的數(shù)量為bi(t),則系統(tǒng)中螞蟻的總數(shù)量為在t時刻從城市i到j的信息素大小為τij, 信息素增量為Δτij(t)。當一只螞蟻從一個城市轉移到另一個城市, 信息素揮發(fā)系數(shù)ε將隨時間不斷的局部更新,Δτij(t)=1/(n×Lnn),Lnn為最鄰近算法求得螞蟻遍歷所有城市的距離。更新規(guī)則滿足(Dorigo, 1997):

當所有螞蟻都完成轉移, 信息素揮發(fā)系數(shù)ρ將隨時間不斷的全局更新, Δτij(t)=1/L,L為螞蟻經過所有城市的最短長度, 更新的規(guī)則滿足:

在t時刻從城市i到j的啟發(fā)函數(shù)為ηij(t)。并且啟發(fā)函數(shù)與距離dij的關系為:

螞蟻轉移規(guī)則為:q為隨機數(shù),q0是0～1 之間的值, 當q≤q0時, 則從城市i轉移到城市j:

當q>q0時, 則選擇下一個城市的j的概率:

allowedk是螞蟻可以選擇下一個沒有經過的城市的集合, 它與螞蟻所經過的城市形成的禁忌表tabuk互為補集。一般q>q0時隨機訪問下一個城市,為了方便應用, 此處更改為訪問allowedk的一個城市。期望啟發(fā)式因子β決定了信息素在路徑中的影響程度, 局部殘留的信息素揮發(fā)系數(shù)ε決定螞蟻選擇其他路徑的概率, 全局殘留的信息素揮發(fā)系數(shù)ρ決定了全局搜索能力及收斂速度。

1.2 Canny 邊緣檢測算子

Canny 雙閾值算子是一種圖像處理技術中常用的利用雙閾值和梯度方向的邊緣檢測算法。它先用高斯函數(shù)對初始圖像進行平滑, 將平滑后的圖像計算梯度大小和方向, 根據梯度差值來尋找圖像邊緣。

設二維高斯函數(shù)為:

將F(x,y)利用高斯函數(shù)進行卷積平滑處理后得到圖像:

采用2×2 鄰域內一階偏導數(shù)有限差分計算I(x,y)的梯度幅值和梯度方向, 將點在x方向和y方向的一階偏導數(shù)分別記為Kx(x,y),Ky(x,y):

之后對梯度幅值進行非極大值抑制, 若M(x,y)比θ(x,y)方向的兩個相鄰像素點的幅值小, 則歸零。若像素點梯度的幅值M(x,y)大于或等于θ(x,y)方向的兩個相鄰像素點的梯度幅值, 則將該點判定為可能為邊緣點。再對該點進行判定, 設置高、低2 個閾值, 大于高閾值的認為是強邊緣點進行保留, 介于兩者之間的認為是弱邊緣點, 若采用8 連鄰域像素檢查周圍有強邊緣點, 則保留該點。

1.3 混合算法處理流程

為了將巖石圖片特征提取缺失的細線狀裂隙進行修復, 本文將會進行以下幾個操作: 首先將巖石片狀圖片轉化為含有色彩信息的像素數(shù)據, 并提取出主要的裂隙特征。利用DBSCAN 聚類算法將裂隙特征點劃分為多個簇, 每個簇利用蟻群算法進行裂隙連接。然后, 將多余的裂隙利用Canny 算子得到圖像求交集去除多余特征點。詳細的步驟如下:

步驟1: 根據巖石片狀圖像像素點灰度均分為多個區(qū)域, 每個區(qū)域賦予一種顏色信息并做標記,提取出裂隙顏色信息。

步驟2: 根據細線狀結構斷裂處與周圍點像素特征, 提取了此次需要的重要像素點, 降低數(shù)據量,提高算法的運行速率。

步驟3: 初始化蟻群參數(shù)β、ε、ρ, 并且根據城市坐標計算兩座城市之間的距離, 對信息素進行初始化。設置迭代次數(shù)Tmax。

步驟4: 設置DBSCAN 聚類算法的鄰域閾值Eps, 點數(shù)閾值Minpts。

步驟5:X1k,X2k表示兩個數(shù)據點k維坐標, 表示兩個點之間距離。當一個點與周圍的點滿足且周圍點數(shù)大于等于Minpts時, 該點為核心點, 以核心點逐漸向外訪問,分為簇C(i)并標記該點到訪問集合v。直到未訪問集合訪問完為止。

步驟6: 將螞蟻隨機分布到n座城市上來進行搜尋, 并進入迭代次數(shù)循環(huán),N=N+1。

步驟7: 螞蟻數(shù)目也開始循環(huán),k=k+1(0

步驟8: 螞蟻根據概率選擇滿足的下一個目的地, 轉移到該目標城市后將城市加入到禁忌表, 直到禁忌表滿結束搜索。

步驟9: 如果k

步驟10: 根據公式(2)對信息素進行滿足蒸發(fā)機制的全局更新, 并再一次進入到步驟6 開始下一次迭代, 直到達到最大迭代次數(shù)。最終找到全局最優(yōu)路徑。

步驟11: 設置灰度閾值提取裂隙圖像, 設置Canny算子雙閾值提取裂隙圖像, 將兩幅圖與裂隙修復圖像求解交集去除多余特征。

2 基于混合算法的巖石圖像分析

2.1 巖石圖像數(shù)據的選取

應用提出的混合算法, 以圖1a 巖石圖像為例進行分析。利用PANDAS 軟件將圖1a 所示的巖石不同特征提取出來, 包含礦物顆粒邊界及其裂隙(為描述方便, 本文以下統(tǒng)稱為裂隙)、巖石基質、巖石礦物, 分別如圖1b～e 所示。由圖1b 細線狀為主的裂隙特征圖像可以看出, 在特征提取的過程中, 雖然保留了主要部分及其特征, 但巖石的裂隙圖像在提取過程中部分缺失, 如圖1b 中紅色框線標注所示。

2.2 蟻群參數(shù)的選取

較好的蟻群算法的參數(shù)可以使算法獲得全局最優(yōu)路徑和較短的運行時間。為驗證本次實驗的有效性,應選取合適的蟻群參數(shù)。根據前人對bayg29, att48等常用數(shù)據庫的研究可以判斷出合適的蟻群參數(shù)(甘屹和李勝, 2011), 對于不同裂隙數(shù), 為獲得全局最優(yōu)路徑, 蟻群算法參數(shù)范圍相同, 但具體的取值不同。根據全局最優(yōu)路徑大小, 選取最合適蟻群參數(shù)需要大量的實驗, 沒有規(guī)律和理論作為參考。以巖石裂隙圖像1b 為例, 可通過大量的實驗分析參數(shù)對算法的影響, 來選取適合的參數(shù)。

該數(shù)值實驗算法主要涉及參數(shù)β、ε和ρ, 對其不同組合進行測試研究。隨著螞蟻數(shù)增加, 搜索能力增強, 但搜索速度呈指數(shù)增加。本文螞蟻數(shù)設置為150, 迭代次數(shù)為200, 圖像均勻分割為三份,Eps設置為15, 缺省參數(shù)設置為q0=0.9。為了防止信息素積累或揮發(fā)過快, 在公式(2)和(3)中,ε和ρ是介于0～1的參數(shù)。在公式(5)和(6)中,β越大計算量越大。本文選取0<β<10, 0<ε<1, 0<ρ<1 進行實驗。β以1 開始步長2 遞增、ε和ρ以0.1 開始步長0.2 遞增, 調整蟻群參數(shù), 判斷對全局最優(yōu)路徑和運行時間的影響。

數(shù)值實驗運行結果顯示, 參數(shù)ρ固定時,β和ε兩個值的變化影響全局最優(yōu)路徑和運行時間(圖2)。ρ=0.1時, 7≤β≤9, 0.1≤ε≤0.3 時全局最優(yōu)路徑小(圖2a)。ρ=0.1 時, 隨著ε的增大, 算法運行時間不變, 運行時間隨著β的增大逐漸增加, 當β≥3 時, 運行時間趨于穩(wěn)定(圖2b)。

而參數(shù)ε固定時,β和ρ兩個的值變化影響全局最優(yōu)路徑和運行時間(圖3)。ε=0.1 時, 7≤β≤9, 0.1≤ρ≤0.7 時全局最優(yōu)路徑小(圖3a)。ε=0.1 時, 隨著ρ的增大, 算法運行時間不變, 運行時間隨著β的增大逐漸增加, 當β≥3 時, 運行時間趨于穩(wěn)定(圖3b)。

參數(shù)β固定時,ρ和ε兩個的值變化也影響全局最優(yōu)路徑和運行時間(圖4)。β=9 時, 0.1≤ρ≤0.7, 0.1≤ε≤0.5 時全局最優(yōu)路徑小(圖4a)。β=9 時, 隨著ρ和ε的增大, 算法運行時間不變(圖4b)。

綜合以上研究, 參數(shù)β、ε和ρ對全局最優(yōu)路徑都有一定的影響, 總的來說, 7≤β≤9, 0.1≤ε≤0.3, 0.1≤ρ≤0.7 范圍內全局最優(yōu)路徑較小。參數(shù)ε和ρ對運行時間影響非常小,β對運行時間影響大, 但β≥3后運行時間趨于穩(wěn)定。并且可以看出, 當全局最優(yōu)路徑較小時, 不可避免要消耗一定的時間。

為了防止不同裂隙之間相互干擾, 應選取較短的全局最優(yōu)路徑連接巖石裂隙, 并基于合適的蟻群算法參數(shù)進行實驗。針對圖1b 部分特征缺失和產生多余特征的情況進行研究, 根據大量實驗得到圖2～4, 本文基于β=9,ε=0.1,ρ=0.3 蟻群算法參數(shù), 從巖石圖像裂隙連接和修復兩個方面, 對實驗結果進行分析和討論。

圖1 圖像特征提取Fig.1 Image feature extraction

2.3 巖石圖像裂隙連接和修復

圖2 ρ=0.1 時, β 和ε 的變化對全局最優(yōu)路徑的影響(a)和β 和ε 的變化對運行時間的影響(b)Fig.2 Th e influence of β and ε on the global optimal path (a) and the influence of β and ε on the running time (b)(ρ=0.1)

圖3 ε=0.1 時, β 和ρ 的變化對全局最優(yōu)路徑的影響(a)及β 和ρ 的變化對運行時間的影響(b)Fig.3 Th e influence of β and ρ on the global optimal path (a) and the influence of β and ρ on the running time (b)(ε=0.1)

圖4 β=9 時, ρ 和ε 的變化對全局最優(yōu)路徑的影響(a)及ρ 和ε 的變化對運行時間的影響(b)Fig.4 The influence of ρ and ε on the global optimal path (a) and the influence of ρ and ε on the running time (b)(β=9)

為了合理地判斷實驗效果, 本文通過定性和定量兩種方式進行對比。定性即通過對比蟻群算法、混合算法連接和修復的圖像與人工連接和修復的圖像。定量則通過利用誤差來驗證蟻群算法和混合算法圖像修復效果, 誤差公式為:

其中, 期望圖像是人工手動修復圖像, 將圖像修復像素數(shù)與理想修復像素數(shù)建立誤差公式, 誤差絕對值越小則修復效果越好。當誤差為負且絕對值較大時, 裂隙大部分未修復; 當誤差為正且絕對值較大時, 產生許多多余的特征。

本實驗選取圖1b 進行修復, 由于在實驗過程中DBSCAN 算法的鄰域閾值Eps對修復效果影響最大,將比較不同Eps對圖像的影響(點數(shù)閾值(Minpts)設置為固定參數(shù)2 不變)。對比蟻群算法和混合算法實驗效果差異, 實驗結果如圖5、6 所示。

圖5 不同鄰域閾值Eps 對巖石圖形修復的影響(紅色裂隙為修復產生)Fig.5 The influence of different neighborhood threshold Eps on rock image restoration

圖6 基于巖石圖像的蟻群算法和混合算法修復誤差對比Fig.6 Comparison of r ock image repair errors of ant colony algorithm and hybrid algorithm methods

蟻群算法結果顯示, 隨著Eps的增加, 誤差絕對值先變小后變大(圖6)。Eps=5 時, 誤差為負, 許多斷裂特征未連接(圖5b、6), 說明Eps過小會導致大部分斷裂特征無法連接。Eps=8 時, 誤差為0, 部分斷裂無法連接, 且產生多余特征(圖5c、6), 說明誤差為0 只是因為期望圖像和修復圖像點數(shù)相同,但不是理想情況, 并且Eps較小會導致部分斷裂無法連接, 且產生部分多余特征。Eps≥16 時, 誤差為正且大于100%, 斷裂特征基本連接但也有多余特征產生(圖5d、e、6), 說明Eps適中或偏大時雖然斷裂特征基本連接但產生許多多余特征。結果證明,蟻群算法雖然可以修復裂隙, 但產生的修復圖像干擾過多, 與理想裂隙修復圖像相似度低。

混合算法結果顯示, 隨著Eps的增加, 誤差絕對值也先變小后變大。Eps=5 和Eps=8 時, 誤差為負且絕對值較大, 許多斷裂特征未連接(圖5f、g, 6);說明使用混合算法時,Eps小會導致部分斷裂特征無法連接, 且Eps越小連接的特征越少。Eps=16時, 誤差為0, 斷裂特征基本連接但也有多余特征產生(圖5h、6), 說明誤差為0 只是因為期望圖像和修復圖像點數(shù)相同, 但不是理想情況,Eps適中時連接效果好,產生極少的多余特征, 且混合算法中誤差已經能夠較準的反映斷裂連接狀況。Eps=30 時, 誤差為正且小于50%, 斷裂特征大部分連接但也有多余特征產生(圖5i、6), 說明混合算法產生的多余特征少, 但Eps過大也會導致部分特征不能修復。結果證明, 混合算法修復裂隙且產生干擾少, 與理想裂隙圖像相似度高。

結合圖5、6 對比蟻群算法和混合算法可知, 混合算法降低了蟻群算法的誤差, 且蟻群算法產生的誤差越大, 混合算法能夠降低誤差的差值越大。這說明Canny 算子去除多余特征, 優(yōu)化了蟻群算法連接的圖像。同時, 相比于蟻群算法, 混合算法修復圖像更接近理想裂隙圖像。

為了驗證混合算法的適用性, 本文選取了裂隙更復雜的巖石圖片進行實驗。由于蟻群參數(shù)的范圍具有適用性, 在范圍內, 這里依然選取β=9,ε=0.1,ρ=0.3 蟻群算法參數(shù), 復雜巖石圖像實驗結果詳見圖7、8。

在復雜巖石裂隙圖像中, 單獨使用蟻群算法結果顯示,Eps=6 時, 誤差為0, 部分斷裂特征無法連接且產生部分多余特征(圖7c、8);Eps=10 時, 誤差很大, 大部分特征修復, 同時產生部分多余特征,這說明Eps適中時, 蟻群算法效果也不理想(圖7d、8)。Eps>10 時或Eps過大時, 誤差很大, 斷裂特征基本連接, 但產生很多多余特征(圖7e、f, 8)。

在復雜巖石裂隙圖像中, 使用混合算法時結果顯示,Eps過小時, 誤差為負, 大部分斷裂特征無法連接(圖7g、8,)。Eps=10 時, 誤差為0, 斷裂特征基本連接, 且只有少部分多余特征產生(圖7h、8)。Eps過大, 誤差較大, 特征大部分連接且混合算法產生的多余特征少, 但Eps過大也會導致部分特征不能修復(圖7i、j, 8)。

由圖7、8 可知, 混合算法誤差小, 使用混合算法可以有效避免間斷和去除多余的修復特征, 效果比單獨使用蟻群算法更好。由于混合算法使用Canny 算子改進, 裂隙修復的結果得到了提升, 基本去除多余特征, 更接近理想的結果。總之, 實驗結果證明, 混合算法同樣適用于復雜的裂隙圖像中,具有普適性。

3 結 論

在巖石圖像修復過程中, 為了識別出裂隙等重要結構并降低數(shù)據量, 本文提出基于蟻群算法、DBSCAN 聚類算法和Canny 算子的混合圖像特征修復算法, 來改善蟻群算法在巖石圖像特征修復中存在的效率低、易產生多余特征的缺點, 較準確恢復巖石圖像中缺失的裂隙特征。文中首先通過數(shù)值實驗發(fā)現(xiàn)在7≤β≤9, 0.1≤ε≤0.3, 0.1≤ρ≤0.7 范圍內全局最優(yōu)路徑較小;β對螞蟻搜索時間有顯著影響,大于等于3 時螞蟻搜索時間增加。采用蟻群算法與混合算法得到的修復圖像和期望圖像進行對比發(fā)現(xiàn),混合算法修復裂隙和去除多余特征的能力優(yōu)于蟻群算法。實際應用結果證明, 混合算法也可以用于復雜巖石裂隙圖像, 提高了圖像修復的準確度, 驗證了該算法的有效性。

圖7 不同鄰域閾值Eps 對復雜巖石圖形修復的影響(紅色裂隙為修復產生)Fig.7 The influence of different neighborhood threshold Eps on complex rock image restoration

圖8 基于復雜巖石圖像的蟻群算法和混合算法修復誤差對比Fig.8 Comparison of complex rock image repair errors of ant colony algorithm and hybrid algorithm methods