基于MATLAB編程的碎屑巖粒度計(jì)算對(duì)比:以江西省會(huì)昌盆地上白堊統(tǒng)周田組為例

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

(1.東華理工大學(xué)地球科學(xué)學(xué)院, 南昌 330013; 2.東華理工大學(xué), 核資源與環(huán)境國(guó)家重點(diǎn)實(shí)驗(yàn)室, 南昌 330013)

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

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

現(xiàn)基于MATLAB編程軟件,以江西省會(huì)昌盆地上白堊統(tǒng)周田組碎屑巖為例,統(tǒng)計(jì)并計(jì)算碎屑巖粒度參數(shù),并將計(jì)算結(jié)果與傳統(tǒng)薄片粒度圖像分析法進(jìn)行對(duì)比,借助計(jì)算機(jī)編程語(yǔ)言對(duì)碎屑巖薄片進(jìn)行圖像處理及碎屑顆粒測(cè)量能夠較好地提高傳統(tǒng)薄片粒度分析過(guò)程中的效率及精度。

1 區(qū)域地質(zhì)背景

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

圖1 會(huì)昌盆地地質(zhì)略圖Fig.1 Geological map of the Huichang Basin

2 粒度分析方法

為了提高粒度分析實(shí)驗(yàn)的準(zhǔn)確性和高效性,現(xiàn)以粒度分析為基礎(chǔ),使用MATLAB軟件進(jìn)行編程,計(jì)算會(huì)昌盆地周田組碎屑巖各項(xiàng)參數(shù)值,并與傳統(tǒng)薄片粒度圖像分析法所得結(jié)果進(jìn)行對(duì)比。MATLAB直接計(jì)算法主要包括以下6個(gè)步驟(圖2)。

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

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

圖3 MATLAB圖像處理結(jié)果Fig.3 MATLAB image processing results

(2)灰度化:將原始RGB圖像轉(zhuǎn)換為灰度圖像,這一步驟是為了能更好地識(shí)別每一個(gè)顆粒的形態(tài)學(xué)要素,以便于后續(xù)的處理[圖3(b)]。

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

(4)圖像增強(qiáng):盡管圖像經(jīng)過(guò)二值化處理后,圖像中仍然存在較多黑點(diǎn)(碎屑巖膠結(jié)物中的雜基或碎小顆粒),在MATLAB運(yùn)行環(huán)境中將其近似地看作是一種隨機(jī)的噪聲,因此還需要通過(guò)濾波器將其去除。本文中所使用的濾波器為自適應(yīng)中值濾波,是在中值濾波的基礎(chǔ)上,根據(jù)預(yù)設(shè)的條件,動(dòng)態(tài)地改變?yōu)V波器的窗口尺寸,以達(dá)到去除噪聲和保護(hù)細(xì)節(jié)的效果[圖3(d)]。

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

(6)彩色標(biāo)注并獲取參數(shù):通過(guò)MATLAB函數(shù)庫(kù)中regionprops函數(shù)有效獲取碎屑顆粒數(shù)目、面積、最大直徑、周長(zhǎng)等參數(shù)[圖3(f)],可用作于顆粒的沉積環(huán)境判別式判斷沉積環(huán)境及磨圓度計(jì)算。

主程序代碼如下所示。

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

img11=rgb2gray(f)

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

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

ff =imgbwI;

image_gray = imgbwI;

%以下為自適應(yīng)中值濾波,ff為自適應(yīng)濾波后的圖像

ff(:) = 0;

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

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 的盤(pán)形工具對(duì)二值化圖像進(jìn)行開(kāi)運(yùn)算

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

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

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

graindata=regionprops(labeled,′all′);

imshow(rgb_label)

3 MATLAB圖像粒度結(jié)果統(tǒng)計(jì)

3.1 顆粒形態(tài)學(xué)參數(shù)計(jì)算

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

表1 薄片粒度統(tǒng)計(jì)結(jié)果Table 1 Results of particle size statistics of flakes

3.2 粒度結(jié)果計(jì)算對(duì)比

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

表2 砂巖粒度參數(shù)計(jì)算對(duì)比Table 2 Comparison of calculated sandstone particle size parameters

圖4 粒度參數(shù)間散點(diǎn)關(guān)系Fig.4 Scatter relationship between particle size parameters

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

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

3.3 磨圓度統(tǒng)計(jì)計(jì)算

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

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

3.4 沉積環(huán)境判斷

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

表3 沉積環(huán)境判別公式計(jì)算結(jié)果對(duì)比Table 3 Comparison of the calculation results of the deposition environment discriminant formula

4 討論

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

5 結(jié)論

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

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

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