基于正交偏光序列消光特征的巖石薄片顆粒分割與孔隙提取

張 欣，張 棟，楊 倬，高世臣，張 艷

(1. 中國地質大學(北京) 數理學院，北京 100083； 2. 中國石油長慶油田分公司 勘探開發研究院，陜西 西安 710018；3. 中國地質大學(北京) 地球物理與信息技術學院，北京 100083)

利用巖石薄片進行鏡下巖石顆粒特征分析、巖石成分鑒定、孔隙類型識別及其發育程度的表征是復雜巖性油氣藏沉積特征、成巖作用以及儲層綜合評價等研究工作中極其重要的技術手段(Moorhouse，1986； 林培英，2005)。近年來隨著巖性油氣藏勘探開發工作的深入發展，大量的取心資料一方面需要快速地完成鑒定分析工作，另一方面沉積、成巖極為復雜的致密儲層評價對巖石薄片鑒定在定量化、精確化方面又提出了更高的要求(張厚福等，2008； 李富恒等，2009)。傳統的巖石薄片鑒定通常是專業人員交互式的手工操作進行識別和統計，以定性-半定量為主，且耗時較多。自20世紀60年代至今，數字圖像處理領域一直在發展與推進，已在醫學和空間項目的基礎上用于更廣泛的范圍，并且在地質領域的巖相學問題、油氣勘探與預測、油氣儲量計算等方面也發揮了重要作用，早在1993年的油氣開發領域中便已得到應用。利用巖石切面或者巖石薄片的圖像分析，進而獲取巖石參數或進行巖石特征推斷也進一步促進了數字圖像處理技術在地質學領域的應用。數字圖像處理技術的引入使得巖石礦物的鏡下鑒定和分析在一定程度上得以簡化，工作效率也大為提高(吳擁等，2013； Chauhanetal., 2016)。但是巖性油氣藏致密儲層的成巖極為復雜，通常壓實作用很強，膠結作用嚴重，原生巖石顆粒之間接觸緊密、邊緣模糊，而后期溶蝕作用的改造使得部分巖石顆粒輪廓模糊，與孔隙之間的邊緣變得復雜。巖石顆粒溶蝕殘余、孔隙內的自生填隙物等在鏡下薄片圖像中形成較多的噪聲點，而傳統的邊緣提取算子、閾值分割技術大多利用灰度信息及鄰域信息進行邊緣提取，在面對高信息量、高復雜度的數據時具有局限性，巖石顆粒分割和孔隙提取效果較差。本文提出了一種新的巖石薄片圖像采集和數字化處理分析方法，以期改善巖性油氣藏復雜儲層巖石薄片顆粒分割和孔隙提取效果，提高復雜巖石薄片的鑒定效率和可靠性。

1 角度域正交偏光序列圖采集方法

巖石薄片的制作通常是將垂直巖石層理方向的小塊巖樣磨制成厚度約0.03 mm的薄片，固定于載玻片上，并覆蓋蓋玻片。通過光學顯微鏡觀察可鑒定分析巖石薄片中的礦物、孔隙等特征(葉潤青等，2011； Freezetal., 2012)。在單偏光下觀察時，光線通過位于載物臺下的下偏光鏡，射入薄片樣品的平面偏光為單一振動方向，多數透明礦物顆粒無明顯的顏色顯示，巖石顆粒之間的邊緣特征不顯著；如果增加位于物鏡之上的上偏光鏡，在正交偏光觀測模式下，巖石礦物顆粒顯示出明顯的光學特征差異，例如礦物顆粒的干涉色、紋理特征、顆粒邊緣凸起差異等。從晶體光學原理分析，這種差異性取決于礦物的成分、顆粒方位等因素，通過旋轉載物臺，連續改變巖石礦物顆粒光率體方位與正交光之間的夾角，可以跟蹤巖石顆粒的消光變化特征，進而綜合多角度位置上的特征判別顆粒邊緣及其礦物成分，以提高巖石薄片的鑒定可靠性。通過與顯微鏡連接的照相機進行多角度正交偏光的巖石薄片圖像采集，進而采用圖像數字化處理技術進行巖石顆粒分割、粒度分析、礦物成分鑒定和孔隙提取等儲層特征綜合研究工作。

測試巖石樣品為鄂爾多斯盆地某區長6段砂巖薄片，前期人工鑒定主要巖石顆粒成分有石英(26.43%)、長石(44.31%)及巖屑(29.26%)，面孔率24.20%。采取連續偏光圖像采集的方法，以10°為間隔連續旋轉載物臺，進行360°圖像采集，共得到36張正交偏光序列圖像，根據礦物在正交光鏡下的消光機理可知，消光周期為90°(Moorhouse，1986； 林培英，2005)，因此36張圖像為4個完整的消光周期(圖1)。圖像采集過程中，顯微鏡放大倍數為4倍，通過圖像角度方位校正對齊處理，獲取的圖像像素為1 000×1 000的圓形視域。

2 巖石礦物顆粒角度域消光特征分析

選取典型的石英、長石、火山碎屑、孔隙各2塊，對目標顆粒0°～360°偏光角度范圍內的同一視域進行標定，視域中心點位置坐標對應圖1中20°偏光圖中的標注點。計算各偏光角度的視域內灰度均值，形成礦物的角度域灰度特征曲線，統一曲線的初始相位，繪制3種類型礦物灰度的角度域特征圖譜(如圖2a所示)，并進一步對礦物顆粒圖像R、G、B三通道特征在角度域的變化特征進行分析(圖2b～2e)。

圖 1 正交偏光顯微鏡下0°～90°周期內巖石薄片序列圖Fig. 1 Multi-angle (0°～90°) microscopic images of thin section (crossed nicols)Q—石英； F—長石； V—火山碎屑； P—孔隙Q—quartz； F—feldspar； V—volcanic debris； P—pore

顆?；叶却嬖诮嵌扔虻闹芷谛宰兓乙?0°為一個消光周期，周期的角度域信息完整表達了不同礦物顆粒的特征。對目標薄片的3種礦物顆粒多次選擇并統計，普遍情況下，長石顆粒的灰度值低于石英顆粒，即鏡下光度較大，巖石顆?；叶戎迪鄬Υ笮№樞驗椋菏?長石>火山碎屑，石英曲線的峰值主要分布在200～255之間，長石曲線的峰值主要分布在75～150之間，火山碎屑的灰度最低，其峰值主要分布在25～65之間。除不同類型顆粒間灰度差異外，顆粒R、G、B通道值存在明顯不同的相互關系。石英顆粒的R通道峰值相比于G、B兩通道明顯較低，G、B通道峰值接近但B通道振幅略?。婚L石與火山碎屑的三通道值差距較大，由高到低按照B、G、R依次排列，其中長石顆粒的G、R通道差異大，而火山碎屑的B、G通道差異大。巖石顆粒R、G、B三通道特征顯示了不同礦物顆粒色調的差異，而對于同一顆粒而言，隨著消光角度的變化，其光度呈現有規律的明暗變化，但其色調始終是保持不變的。巖石顆粒正交偏光鏡下特征在角度域的變化、不同巖石顆粒變化特征的差異性是進行巖石薄片數字化鑒定分析的基礎。

3 正交偏光圖像巖石顆粒分割和孔隙提取

3.1 常規方法提取巖石顆粒邊緣

巖石薄片序列圖中巖石顆粒表面紋理復雜，顆粒內部以及邊緣附近的噪點較多，嚴重影響了邊緣提取的效果，造成邊緣零散，提取結果較粗，巖石顆粒形態、孔隙結構和巖石顆粒溶蝕特征等難以顯現。雙邊濾波在消除顆粒內部及邊緣噪點的同時，對顆粒邊緣特征及孔隙特征弱化程度較低，多數顆粒邊緣更加清晰，因此選擇雙邊濾波的去噪方法對巖石薄片序列薄片進行了預處理。

圖 2 不同巖石組分角度域正交偏光消光譜特征Fig. 2 Spectroscopies in angular domains of different components under crossed nicols microscopya—各巖石組分灰度； b—石英顆粒； c—長石顆粒； d—火山碎屑； e—孔隙內部a—gray level of component； b—quartz； c—feldspar； d—volcanic debris； e—pore

3.1.1 Sobel、Canny算子的單偏光角度邊緣檢測

采用常用的圖像邊緣提取方法Sobel算子和Canny算子進行薄片巖石顆粒邊緣進行提取。圖3是薄片在20°和70°角度時正交光鏡下的巖石顆粒邊緣及輪廓特征圖。

傳統的巖石薄片邊緣提取技術，多基于巖石薄片在單一角度位置的正交光鏡下圖像，提取結果通常無法準確提取全部巖石顆粒的輪廓。從圖3可以看出，幾乎所有礦物顆粒的輪廓特征都會在某些消光角度的圖像中有所顯現，但是同一顆粒的輪廓特征極少在所有消光角度的圖像中均有明顯顯現。因此，采用單一角度正交偏光圖像對巖石薄片進行輪廓提取存在一定的缺陷：單一角度正交偏光圖像中總會有部分巖石顆粒的邊緣特征不明顯，難以提取完整的邊緣信息，顆粒劃分粗糙。對于某一角度正交偏光鏡下，兩相鄰顆粒光度接近時，常規方法很難精確提取顆粒邊緣和輪廓。

3.1.2 基于角度域序列圖像的邊緣融合

針對單一角度正交偏光圖像只能較好地刻畫出部分巖石顆粒的邊緣特征，且某一巖石顆粒的邊緣總會在某些角度的圖像中有明顯顯現的特點，采用角度域序列圖像邊緣提取結果疊加的方式進行邊緣特征的融合。圖4為巖石薄片在0°～90°一個消光周期內的10張圖像R、G、B三通道參數邊緣提取結果各自的疊加圖以及三通道疊加融合圖。

圖 3 基于灰度的單偏光圖像邊緣提取Fig. 3 The edge diagram extracted in gray level of plainlight imagea—Sobel，20°； b—Sobel，70°； c—Canny，20°； d—Canny，70°a—Sobel at 20°； b—Sobel at 70°； c—Canny at 20°； d—Canny at 70°

相比于傳統的單一角度的正交光圖像邊緣提取的結果，基于角度域和顏色通道的邊緣提取和疊加融合方式實現了顏色通道信息和消光角度的縱橫信息聯合拓展。對比Sobel和Canny兩種邊緣提取方法在R、G、B和灰度參數中提取邊緣的疊加融合結果，可以看出不同通道參數展現的顆粒邊緣特征較為一致，而多角度序列圖像的疊加，展現了全部顆粒的邊緣和輪廓，且較為清晰。但是這種多圖像邊緣的直接疊加也導致多圖像噪點的疊加，給巖石顆粒的準確分割帶來困難。

圖 4 Sobel及Canny算子周期內三通道的邊緣提取與疊加Fig. 4 The overlay edge maps of three channels in a cycle with Sobel and Canny operatorsa～d—Soble算子邊緣； e～h—Canny算子邊緣； 按列依次為R、G、B通道，最后一列為三通道疊加結果a～d—the edge diagrams with Sobel operator； e～h—the edge diagrams with Canny operator； each column is followed by R, G, B, and three channel overlay results

3.2 基于角度域消光特征的相關分析方法

3.2.1 序列消光圖像相關分析方法

簡單地對不同角度的邊緣提取結果進行疊加，是對各角度邊緣提取結果取并集的結果，當各角度邊緣圖參數設置不佳時，對最終提取結果有較大影響，若保留較多單角度邊緣信息，會造成最終疊加圖的過分割。該方法雖然在最大程度將邊緣復刻至最終結果圖，但并未充分考慮各角度域圖像之間的內在聯系，無法表達巖石顆粒角度域消光特征的變化規律，無疑是一種有效信息的損失。因此，構建像素點一個消光周期內的灰度向量A(x,y)作為該像素點的完整消光特征表達：

A(x,y)=vec(f0(x,y),f1(x,y),…f9(x,y))

其中(x,y)表示像素點位置，fi(x,y)表示(x,y)在第i幅序列圖的灰度值。

以Pearson相關系數來度量像素點之間的相似性，選定目標像素點，計算目標像素點與其鄰域內的各像素點灰度向量間的相關系數。圖5是對圖像不同顆粒和邊緣101×101像素鄰域相關系數的計算結果，圖5a～5f分別對應圖1中70°偏光圖像中的標注區域。圖中越接近紅色的區域與中心目標點融合效果越好。顏色相近的區域為強相似區域，可作為融合后的分區，一般表征同一顆粒、邊緣或噪聲點的內部。不同顆粒和不同邊緣組合位置的相關性特征如下：

(1) 同一顆粒內部的點之間相關系數大(大于0.75)，呈現明顯的獨立分區，且巖石顆粒分割特征顯著(圖5a)；

(2) 邊緣區域與顆粒內部在相關系數上存在明顯差異，相鄰顆粒形成明顯的分區，而顆粒之間邊緣附近的相關系數則明顯較低(小于0.3)；低亮度的巖石顆粒邊緣特征偏弱，分區特征不明顯(圖5b、5c)；

(3) 巖石顆粒之間的孔隙內部相關特征明顯，孔隙和顆粒之間存在漸變特點(圖5d、5e)；

(4) 溶蝕的長石顆粒邊緣輪廓不明顯，內部結構復雜，溶蝕顆粒殘余部分和完全溶蝕成為孔隙的部分相關系數較高，而孔隙填物分布區相關系數較低甚至負相關(圖5f)。

3.2.2 基于角度域相關的邊緣提取方法

分別選取石英顆粒內部、長石顆粒內部、孔隙結構內部及邊緣點作為目標點，圖6為邊緣提取示意圖，圖6b中的每一個黑色方框區域即為指定的窗口大小。以采用5×5窗口大小的圖6a為例，其中的數字5為目標點，9個數字分別對應于圖7中每個圖形的9條線，最終產生圖7a～7d中4種情形下9點的消光特征雷達圖。

邊緣點作為相鄰兩顆粒的公共邊，所選鄰域跨越兩顆?；蛞活w粒，包含顆粒內部點，因此部分點與目標點相似程度差，為表征這種特點，采用目標點與鄰域內其他點的相關系數均值、標準差作為邊緣提取的評判指標，并進一步構建相關系數均差比指標，方法如下：

圖 5 不同顆粒和邊緣組合位置的鄰域相關性Fig. 5 Neighborhood correlation of different grains and edge patternsa—顆粒內部； b—相交邊緣； c—單側低亮邊緣； d、e—粒間孔隙； f—長石溶蝕區域a—inside grain； b—intersection edge； c—gray edge on one side； d, e—intergranular pores； f—feldspar dissolution zone

圖 6 基于相關的邊緣提取示意圖Fig. 6 Diagram of edge extraction based on correction

向量之間的相關系數計算公式為：

r(A(x,y),A(m,n))=Cov(A(x,y),A(m,n))/

σ(A(x,y))·σ(A(m,n))

目標點相關系數均值：

目標點相關系數標準差：

目標點相關系數均差比指標：

rms(x,y)=rmean(x,y)/rstd(x,y)

其中(x,y)為目標點，δ為正方形窗口，k為窗口大小，(m,n)為窗口內非目標點。

由圖7可知，石英和長石內部均表現為強相關；孔隙內部相關性有擾動，為非均質結構；邊緣相關性變化最大，表現為相關均值較低，相關標準差較高，均差比明顯小于顆粒內部。

圖 7 5×5鄰域9點的消光特征雷達圖Fig. 7 Radiation diagrams of extinction spectroscopies at 9 points in 5×5 neighborhoodsa—石英內部； b—長石內部； c—孔隙內部； d—邊緣a—quartz； b—feldspar； c—pore； d—edge

3.2.3 基于角度域相關的邊緣提取及巖石顆粒分割和孔隙提取

采用上述方法，對巖石薄片在0°～90°一個消光周期內的10張圖像進行相關性分析，并提取相關系數均值和相關系數標準差(圖8a、8b)。對比于傳統的Canny、Sobel算法，基于相關分析的邊緣提取結果輪廓的閉合性更好，弱邊緣的表達更為清晰，孔隙結構更為明顯。

利用相關系數均值圖對巖石顆粒分割和孔隙進行提取。其中巖石顆粒和孔隙內部的相關系數均大于0.75，孔隙部分在正交偏光鏡下圖像的灰度值低于25，且R、G、B三通道參數相近，顆粒與孔隙的邊緣的相關系數均值小于0.3。

圖 8 巖石薄片顆粒分割成果圖Fig. 8 Grain segmentation results of thin sectiona—相關系數均值； b—相關系數標準差； c—巖石顆粒分割圖； d—孔隙分布圖a—mean coefficient of correlation； b—correlation coefficient standard deviation； c—grain segmentation map； d—pore distribution map

3.2.4 效果評價

從邊緣提取效果層面分析，比照圖3，相比于傳統的單偏光方法，從多偏光序列圖像中提取的邊緣信息更加完整；相比于圖4的Canny、Sobel算子，該方法對邊緣的粗細刻畫更為準確，精確到像素級，且有效降低了圖4中算法本身造成的毛細邊緣太多而產生的過分割現象，為進一步的巖石礦物成分識別、定量分析及參數統計打下了良好的基礎。

從定量分析效果層面分析，統計圖8c和圖8d巖石顆粒和孔隙分割圖像中孔隙部分所占比例，即面孔率，為25.36%；采用人工鏡下網格交點法統計的面孔率為24.20%，兩者偏差為1.16%，總體結果較為一致。網格法作為基于人為經驗判斷邊緣進而估計所占區域大小的傳統方法，其統計精度約為網格大小的一半，對微小孔隙、溶蝕顆粒模糊邊緣等區域的統計易產生較大誤差，在精細定量化方面傳統方法略顯薄弱。另一方面，多偏光序列圖像相關分析法盡管在統計精度、定量化和自動化等方面優點突出，但是對巖石成分和孔隙類型的準確人工鑒定依然是對該方法分析結果進行標定和校正的基礎。通過兩種方法的有機結合和相互印證，可以獲得更為準確、更加精確的結果。

4 結論與建議

(1) 巖石薄片角度域序列正交偏光鏡下圖像的數字分析技術，一方面可以提取并分析巖石礦物在完整的消光周期內的連續變化特征，更為全面地揭示出巖石礦物的光學特征，另一方面巖石薄片圖像的RGB三通道分析則可以提供巖石礦物在顏色通道維度上的變化特征，很好地改善了巖石顆粒邊緣特征的刻畫效果、礦物成分鑒定和識別效果，并提高了工作效率。

(2) 綜合角度域正交偏光消光特征的巖石薄片圖像相關分析法，綜合了巖石薄片圖像的角度域消光變化特征和像素點鄰域相關性，與Soble和Canny等邊緣提取算法相比，明顯提升了巖石顆粒邊緣的提取效果，更為準確地刻畫了巖石顆粒的輪廓和孔隙分布特征。

(3) 巖石薄片角度域正交偏光消光譜相關性分析表明，在同一礦物顆粒和孔隙內部的消光譜相似性較強，鄰域相關性較高，多大于0.75，而在礦物顆粒之間的接觸邊緣、顆粒與孔隙相鄰部位則相關系數較低，多小于0.3，甚至負相關。以此為據可以較為準確地分割出礦物顆粒和孔隙。

(4) 巖石薄片正交偏光鏡下的消光特征在角度域上的連續變化，為研究巖石礦物提供了一個新的維度和視角，巖石礦物在角度域上的消光譜特征不僅僅可以用于刻畫巖石顆粒形態、孔隙特征等，同時還提供更為豐富的巖石顆粒形態、巖石礦物成分、孔隙類型及其內在填隙物等信息(兩個巖石顆粒形態)，為巖石薄片的自動化、數字化、定量化的鑒定分析提供了數據和信息基礎。巖石薄片正交偏光消光譜的研究及應用，在沉積巖石粒度和產狀分析、礦物成分定量識別、孔隙結構參數提取、孔隙填隙物識別及巖石沉積、成巖特征研究等方面必將發揮極為重要的作用。

(5) 在對巖石薄片的邊緣提取的研究中，顆粒邊緣的各向異性是一個顯著的信息，且邊緣和顆粒的敏感窗口大小相距較大，即關于邊緣的方向信息與窗口大小對計算效率的提升在本文的研究角度和技術方法中并未探討，基于巖石薄片信息利用的完整性的探討和研究仍需進一步展開。