基于MATLAB編程的碎屑巖粒度計算對比:以江西省會昌盆地上白堊統周田組為例

陳念楠, 李滿根* , 宋志杰, 劉東興, 范鵬飛, 吳思楷, 魏廣富, 劉穎

(1.東華理工大學地球科學學院, 南昌 330013; 2.東華理工大學, 核資源與環境國家重點實驗室, 南昌 330013)

粒度分析是通過研究碎屑巖的粒度大小和粒度分布,判斷沉積巖的沉積環境及水動力條件特征。當碎屑巖顆粒所處沉積環境變化時,其粒度特征也會隨之改變[1-2]。而對粒度進行精確有效的測量是該分析實驗的關鍵,激光粒度分析及顯微薄片粒度分析是地質領域較為常用的粒度分析方法[3]。前者對風沙沉積物、黏土等固結程度較差的樣品具有良好的應用。而后者常測量固結程度較高的砂巖、粗砂巖等樣品的粒徑和含量,該方法優點在于原理和設備簡單,缺點也表現為需要人工大量統計顆粒直徑,后期還需計算粒度參數及繪制粒度曲線,整體操作流程煩瑣且效率低下。

為了提升傳統粒度分析法的工作效率,許多學者將計算機強大的計算統計功能引入到粒度分析中,羅章等[4]通過動態圖相法對不規則海灘砂粒度的測量。張曉曉等[5]基于MATLAB程序使用圖解法和矩陣法快速計算了海州灣潮灘沉積物參數值,均通過了置信水平95%的顯著性檢測。陳麥雨等[6]使用ImageJ軟件對塔里木輪南地區砂巖鏡下薄片照片進行處理,取得了良好的實驗結果。方夢陽等[7]使用MATLAB編程獲取了廣東龍須帶水庫上三疊統碎屑巖顆粒輪廓及各類形態學參數信息。Liu等[8]通過QGrain開源軟件的應用,在粒度端元分析中也取得了不錯的效果。

現基于MATLAB編程軟件,以江西省會昌盆地上白堊統周田組碎屑巖為例,統計并計算碎屑巖粒度參數,并將計算結果與傳統薄片粒度圖像分析法進行對比,借助計算機編程語言對碎屑巖薄片進行圖像處理及碎屑顆粒測量能夠較好地提高傳統薄片粒度分析過程中的效率及精度。

1 區域地質背景

會昌盆地位于江西贛南地區,武夷山隆起帶南段西側(圖1),處于華夏板塊陸緣造山帶,沿著NNE向的石城-尋烏斷裂呈NE向近弧型展布的狹長箕狀盆地[9]。在早白堊世晚期,會昌盆地受太平洋板塊北北西向左行走滑的影響下轉為斷陷盆地。盆地主要蓋層為茅店組(K2m)、周田組(K2z)以及河口組(K2h),其中茅店組為一套紫紅色中粗砂礫巖,發育斜層理及細粒層序,主要沉積于盆地東部。周田組巖性為紫紅色鈣質粉砂巖夾黃色泥巖,發育水平層理,主要分布于盆地中西部。河口組出露于盆地西緣,巖性主要為紫紅色礫巖,與盆地邊緣呈斷層接觸[10]。

圖1 會昌盆地地質略圖Fig.1 Geological map of the Huichang Basin

2 粒度分析方法

為了提高粒度分析實驗的準確性和高效性,現以粒度分析為基礎,使用MATLAB軟件進行編程,計算會昌盆地周田組碎屑巖各項參數值,并與傳統薄片粒度圖像分析法所得結果進行對比。MATLAB直接計算法主要包括以下6個步驟(圖2)。

圖2 MATLAB處理流程圖Fig.2 MATLAB processing flow chart

(1)讀入圖片:將拍攝碎屑巖單偏光鏡下圖像以矩陣的形式導入到MATLAB軟件中[圖3(a)]。

圖3 MATLAB圖像處理結果Fig.3 MATLAB image processing results

(2)灰度化:將原始RGB圖像轉換為灰度圖像,這一步驟是為了能更好地識別每一個顆粒的形態學要素,以便于后續的處理[圖3(b)]。

(3)二值化:本文中使用的是最大類間方差法(OTSU法)來將圖像進行二值化,這是一種能夠根據灰度圖像的薪資自動選取最佳閾值,進而通過最佳閾值對灰度圖鏡下二值化[11]。通過OTSU法對圖像進行處理過后,既能將圖像清晰分割,同時也完整保留了局部細節信息[圖3(c)]。

(4)圖像增強:盡管圖像經過二值化處理后,圖像中仍然存在較多黑點(碎屑巖膠結物中的雜基或碎小顆粒),在MATLAB運行環境中將其近似地看作是一種隨機的噪聲,因此還需要通過濾波器將其去除。本文中所使用的濾波器為自適應中值濾波,是在中值濾波的基礎上,根據預設的條件,動態地改變濾波器的窗口尺寸,以達到去除噪聲和保護細節的效果[圖3(d)]。

(5)形態學處理:濾波后的圖像雖然去除了較多噪聲,但仍有許多細小碎屑附著在巖石顆粒表面及間隙中,也會影響后續的計算統計結果。而通過形態學腐蝕及膨脹的操作可以將碎屑顆粒原始形態學要素保留,處理后的圖像顆粒內部光滑,邊界清晰,不僅能夠更好地被計算機識別,同時也還原了原始薄片中顆粒的形態[12][圖3(e)]。

(6)彩色標注并獲取參數:通過MATLAB函數庫中regionprops函數有效獲取碎屑顆粒數目、面積、最大直徑、周長等參數[圖3(f)],可用作于顆粒的沉積環境判別式判斷沉積環境及磨圓度計算。

主程序代碼如下所示。

f=imread(′D: ZX07.jpg′);

img11=rgb2gray(f)

X=graythresh(img11);%X=閾值,通過確定閾值對圖像進行二值化操作

imgbwI=im2bw(img11,X);%二值化圖像

ff =imgbwI;

image_gray = imgbwI;

%以下為自適應中值濾波,ff為自適應濾波后的圖像

ff(:) = 0;

alreadyProcessed = false(size(image_gray));%生成未處理的邏輯矩陣(進行初始化)

Smax = 7;% 迭代.

for k = 3∶2∶Smax

zmin = ordfilt2(image_gray,1,ones(k, k),’symmetric′);% 最小值濾波器

zmax = ordfilt2(image_gray,k*k,ones(k, k),’symmetric′);% 最大值濾波器

zmed = medfilt2(image_gray,[k k],’symmetric′);% 中值濾波器

processUsingLevelB = (zmed >zmin) &(zmax >zmed) &～alreadyProcessed;

zB = (image_gray >zmin) &(zmax >image_gray);

outputZxy = processUsingLevelB &zB;

outputZmed = processUsingLevelB &～zB;

ff(outputZxy) = image_gray(outputZxy);

ff(outputZmed) = zmed(outputZmed);

alreadyProcessed = alreadyProcessed | processUsingLevelB;

if all(alreadyProcessed(:))

break;

end

end

ff(～alreadyProcessed) = zmed(～alreadyProcessed);

I3=imopen(ff,strel(′disk′,15));%使用半徑為 15 的盤形工具對二值化圖像進行開運算

I4=imfill(I3,′holes′);%充填孔隙

[labeled,numObjects]=bwlabel(I4,4);%獲取彩色標注圖像

rgb_label=label2rgb(labeled,@spring,′c′,′shuffle′);

graindata=regionprops(labeled,′all′);

imshow(rgb_label)

3 MATLAB圖像粒度結果統計

3.1 顆粒形態學參數計算

經過MATLAB軟件統計后的碎屑顆粒均被標注及著色[圖3(f)],顆粒粒徑單位為像素。將粒徑統計結果依比例尺換算成真實粒徑d,根據弗里德曼[13]粒度回歸校正方程D=0.381 5+ 0.902 7d將粒徑矯正,使用公式φ=-log2D將粒徑轉換為粒級φ。最后在Origin在軟件中計算累計含量分別為5%、16%、25%、50%、75%、84%、95%對應的φ并繪制頻率累計曲線和概率累計曲線[14]。通過MATLAB圖像粒度共統計得到顆粒數151顆,并按照0.25φ為統計區間,劃分φ粒級(表1)。

表1 薄片粒度統計結果Table 1 Results of particle size statistics of flakes

3.2 粒度結果計算對比

通過Folk-Ward公式[15]計算碎屑顆粒的平均粒級、標準偏差、偏度和峰度等粒度參數(表2)。以直接計算法所得結果作為基準,圖解法所得結果中平均粒度誤差22.57%,標準偏差誤差0.375%,偏度誤差78.97%,峰度誤差0.548%。兩種計算方法結果顯示:標準偏差和峰度的計算結果均較為相近,結果上能夠互相替換。平均粒級和偏度的計算結果偏差較大,無明顯相關性,因此不能互相替換。通過粒度間的散點關系[圖4(a)～圖4(c)]圖解法和直接計算法均顯示周田組分選程度中等,尖銳程度基本為中等。而圖解法顯示樣品偏度呈近似對稱,粗細顆粒組分大抵相等,而直接計算法顯示為負偏態,即沉積物粒度偏細粒一側,顯示周田組以細粒組分為主,弗里德曼離散圖[圖4(d)]顯示周田組主要形成于濱淺湖及水下重力流的環境下,與該區域地質背景判斷沉積環境一致。

表2 砂巖粒度參數計算對比Table 2 Comparison of calculated sandstone particle size parameters

圖4 粒度參數間散點關系Fig.4 Scatter relationship between particle size parameters

從粒度頻率曲線可以看出,樣品的圖解法頻率曲線和直接計算法頻率曲線近似對稱[圖5(a)],呈現中等峰度,沉積物粒徑范圍均處于φ～6φ。圖解法直方圖曲線呈現近似對稱的雙峰馬鞍形狀,反映了沉積過程中存在兩種沉積組分,且二者含量相近[10]。而MATLAB直接計算法直方圖曲線呈現單峰、微負偏態,說明直接計算法分選中等至較好,組分單一?？傮w來看,周田組的沉積環境相對簡單,沉積物搬運方式固定,水流動力相對穩定。兩種計算方法所得累計曲線均顯示出滾動-跳躍-懸浮的三段式[圖5(b)],其中滾動次總體含量少、斜率低、分選性差。跳躍次總體約占85%,跳躍和滾動次總體的截點位于1.25φ,且跳躍總體斜率為40°～50°,分選中等。懸浮次總體斜率相對較高,分選性較好,反映了水流條件強于滾動和跳躍次總體。通過對比圖解法和直接計算法結果,發現兩者非常接近。不過由于圖解法在現有的測量條件下無法統計所有顆粒的粒度值,因此直接計算法借助數字圖像處理技術,可以統計鏡下所有大于15 μm以上的顆粒粒徑,具有更高的統計粒度數目和精度。

圖5 MATLAB計算法與圖解法粒度曲線圖對比Fig.5 Comparison of particle sizes curves between MATLAB calculation method and graphic method

3.3 磨圓度統計計算

磨圓度是指碎屑顆粒的原始棱角被磨圓程度的度量。Powers[16]制作了圓度圖板,用于反映顆粒組成源、搬運距離和沉積環境等信息。本文中采用IPP(image pro plus)軟件中的圓度計算公式(c=L2/4πA,其中c為磨圓度,L為碎屑的周長,A為礦物顆粒面積)來衡量顆粒實際面積與長軸外接圓面積之比。利用MATLAB圖像分析功能獲取和統計碎屑巖圖像中的各碎屑顆粒參數,并批量計算每個顆粒的圓度。在選擇的周田組樣品中,大部分顆粒的圓度范圍為[1,1.5](圖6),對應Powers圓度圖版中的碎屑顆粒呈棱角狀至次棱角狀。根據偏光顯微鏡分析結果顯示,周田組碎屑巖樣品中碎屑的含量較高,整體顆粒呈線狀接觸、顆粒支撐結構,磨圓度較差,以棱角狀-次棱角狀為主。這一結果與MATLAB直接計算法得出的結果一致,從側面反映了MATLAB軟件在磨圓度計算中具有較高的精度[9]。

圖6 周田組碎屑巖顆粒IPP圓度分布直方圖Fig.6 Histogram of IPP roundness distribution of clastic rocks of Zhoutian Formation

3.4 沉積環境判斷

判別分析作為一種多元分析方法,通過統計碎屑巖的粒度大小可以反推當時的地理環境,以及能鑒別不同成因類型的沉積物。目前常用薩胡公式[17]來判斷沉積環境,然而該判斷函數是根據有限的現代沉積物樣本做出的,在判別環境上會存在一定的偏差。李昌志等[18]建立核心區域圖解法和分析參數分布特征和聯系,對泥石流、冰磧和河湖相沉積物的粒度特征參數進行補充。將周田組的粒度參數代入計算,結果顯示,圖解法與直接計算法所得結果均顯示周田組形成于河流沉積環境中(表3),與地質背景判斷沉積環境一致,也證明直接計算法的精度高,適用性強,兩種計算方法可以替代使用。

表3 沉積環境判別公式計算結果對比Table 3 Comparison of the calculation results of the deposition environment discriminant formula

4 討論

利用MATLAB軟件數字圖像處理技術能力,能夠精確統計碎屑巖中顆粒的粒徑,完整還原了顆粒及膠結物的原始面貌。在時效性方面,傳統鏡下薄片統計法耗時超過1 h,而MATLAB 直接計算法運行時間僅30 s。極大地減少了人為工作量,提升了粒度分析的效率,尤其是避免了人為主觀因素的干擾?；谟嬎銠C強大的計算能力,可以快速獲取碎屑巖中顆粒粒徑、周長、面積、質心等參數,這是傳統鏡下薄片統計法所不能比擬的,這些參數也為后續進行沉積巖磨圓度、組分分析、顆粒遷移距離以及判斷沉積環境等研究提供重要依據。盡管計算機編程技術在粒度分析中具有較高的可行性,但仍然存在一些劣勢及不足,主要表現為:受顯微鏡分辨率及計算機識別能力的影響,計算機無法區分出及識別顆粒粒徑小于 15 μm 的基質與膠結物的邊緣形態,因此該方法適用于顆粒粒徑較大的碎屑巖樣品。同時本文選取的樣品為細粒石英長石砂巖,碎屑成分簡單,而對于成分較為復雜的碎屑巖適用性如何,這仍然是一個值得關注和討論的問題[19]。

5 結論

(1)MATLAB直接計算法與圖解法所求得平均粒徑分別為4.008φ和3.103φ;標準偏差分別為0.8φ和0.803φ;偏度分別為-0.195φ和-0.041φ;峰度分別為0.911φ和0.916φ。兩種計算方式中標準偏差和峰度誤差較小,可以互相替代,而標準偏差和偏度誤差程度較大,不可互相替代。

(2)MATLAB直接計算方法與鏡下巖礦鑒定均顯示周田組碎屑顆粒呈棱角狀-次棱角狀,整體磨圓度較差。通過沉積環境判別公式反映周田組形成于河流沉積及湖相沉積環境中,與地質背景及顯微觀測結果一致。

(3)通過MATLAB編程軟件的強大算力,能夠快速、精準地識別碎屑巖中顆粒輪廓及顆粒粒徑、周長、面積、質心等參數,極大地減少人為工作量,提升了粒度分析的效率,同時也避免了人為主觀因素的干擾。