基于顆粒識別分析系統的碎屑流堆積物顆粒識別和統計方法研究

陳 達，許 強，鄭 光，彭雙麒，王 卓，何 攀

（成都理工大學地質災害防治與地質環境保護國家重點實驗室，四川 成都 610059）

碎屑流是一種高速遠程潰散性滑坡，具有非常高的運動速度和極大的位移，速度在30 m/s 以上[1]，具有極大的破壞力。為研究碎屑流運動特性及運動機理，有必要進行堆積物粒度定量分析，這也是碎屑流災害研究的重點。傳統顆粒的測量主要依靠人工，如需通過生成高分辨率數字地表模型（Digital Surface Model，DSM）對堆積物粒度進行人工解譯[2-4]，存在效率低、精度不高以及人為因素誤差大等不足。

計算機數字圖像處理技術具有自動化、高效快捷等特點，被廣泛運用于顆粒識別。吳義祥[5]研究了土細觀結構圖像的定量分析系統，可以獲取細小顆粒面積、最長弦長度等基本要素；施斌等[6]使用Videolab圖像處理系統直接分析SEM 等照片，定量處理土顆粒微細形狀、大小和比例關系；涂新斌等[7]通過Qwin軟件處理顆粒二值圖像，獲得了顆粒數、長度、寬度等顆粒數字特征；梁雙華等[8]提出Mapinfo-Photoshop計算方法，利用Photoshop 圖像處理功能和Mapinfo的統計分析功能，對土顆粒圖像進行信息提取，定量表征顆粒、孔隙的個數以及周長、平均孔徑等參數，極大地降低了圖像分析費用，提高效率；Bai Baojun 等[9]通過Imagej 軟件，直接將頁巖原始孔隙圖像進行二值化處理，進而定量統計孔隙、顆粒發育情況，得出孔隙率。可以看出，這些研究都是為了定量化分析圖像的細觀特征，實現圖像自動識別，得出孔隙、顆粒的結構參數和數量特征。但是，這些方法在以下關鍵問題上存在一定的缺陷：（1）設定研究最小顆粒粒徑及較小粒徑顆粒識別方法；（2）一個圖像中可能有上千個顆粒和孔隙，形狀復雜，分布不均，宏觀特征參數已經不能滿足研究需要，為此需要一些特定的參數，分析每個顆粒的形狀、大小，以此來反映顆粒系統特征。

孔隙（顆粒）及裂隙圖像識別與分析系統（Pores and Cracks Analysis System，PCAS），可以定量分析碎屑流堆積物粒度特征，準確獲得堆積物顆粒結構的幾何參數和統計參數，且操作簡單，具有自動化和可重復等優點。PCAS 在顆粒識別中，存在著三個參數的選取問題：閾值（T）、孔喉封閉半徑（r）、最小孔隙面積（S0）。該方法在滑坡堆積物粒度分析中具有很好的效果，但實用范例很少，目前僅彭雙麒等[10-11]通過PCAS 系統研究了碎屑流堆積體，得出碎屑流堆積物顆粒分布的一般規律，但沒有具體闡述PCAS顆粒識別系統的工作原理、參數意義以及參數選取方法。為此，本文通過PCAS 軟件對貴州納雍崩塌碎屑流堆積物顆粒特征進行分析，得出圖像識別系統中閾值、孔喉封閉半徑、最小孔隙面積合適的參數取值，為PCAS 在崩塌堆積物顆粒的定量研究中提供了一個高效可行的方案。

1 方法介紹

PCAS 圖像識別步驟見圖1。

首先將圖像導入PCAS 軟件，調整閾值（T）進行二值化處理，區分出孔隙和顆粒；然后選取合適的孔喉封閉半徑（r）和最小孔隙面積（S0），最后孔隙和顆粒的各種幾何參數將匯集于數據表，包括個數、長度、寬度、面積、形狀系數、定向性等，并統計得到顆粒和孔隙的分形維數、表觀孔隙率、面積概率分布指數等統計參數，實現孔隙和顆粒結構的定量分析。本次主要研究顆粒個數和長度兩個參數。

PCAS 得出的參數都是以像素為單位，可以通過換算得到實際的幾何參數，如式（1）、（2）[12]：

式中：R—圖像的分辨率；

S、C—像素面積和像素周長；

St、Ct—真實面積和真實周長。

PCAS 主要由三個參數確定：閾值（T）、孔喉封閉半徑（r）、最小孔隙面積（S0）。本節主要通過介紹這三個參數的選取意義，分析PCAS的原理與方法。

1.1 閾值（T）

首先獲取崩塌堆積物的正射影像圖，圈定分析的范圍，獲得適當的圖片大小，即開始圖像識別。在圖像識別前，需要進行二值化處理，獲得二值圖像，PCAS 采用閾值分割法進行圖像二值化[13]，所謂閾值分割法，是基于影像圖中各個組成要素的灰度值的不同，通過選擇適當的灰度閾值將影像圖轉換為二值圖像，根據圖像中孔隙與顆粒的灰度值不同加以區分[14]。具體來說，圖片由像素組成，而每一個像素是由紅綠藍（RGB）三原色分量組成，軟件自動選取一個特征顏色值作為灰度值X（RGB），計算圖像中各個部分的灰度值P（RGB）與這個特征灰度值的距離[12]：

與選取的閾值作對比，d小于閾值為黑色，識別為孔隙；大于閾值為白色，識別為顆粒；從而把圖片識別成灰度圖像（圖2b）。

具體步驟為：起初設定較大的閾值，此時較多為顆粒的單元被識別為孔隙（黑色），再慢慢地減小閾值，直到黑色的細小顆粒單元被轉化為白色，能清晰地識別出顆粒與孔隙，此時為較好的閾值。為了減少人為判斷誤差，通常對不同堆積區進行二值化處理，求取閾值平均值，作為選定的最終閾值。

1.2 孔喉封閉半徑（r）

顆粒與顆粒間并不是完全分開，會有連接（圖2b），在二值化過程中，顆粒通過細小連接，會被錯誤識別為一個顆粒。因此，PCAS 在識別過程中，會設定一個特定直徑的孔喉封閉半徑，定義為腐蝕結構元素的半徑r，對圖像做腐蝕運算[15]。當顆粒之間連接的直徑小于2r時，則區分為兩個獨立區域，分開顆粒（圖2c），再將剩余像素歸并到顆粒上，實現顆粒與顆粒的自動區分[16]。

圖2 PCAS 在T/r/S0=170/1/10時的顆粒識別結果Fig.2 PCAS particle recognition results at T/r/S0=170/1/10

1.3 最小孔隙面積（S0）

圖片是由一系列像素點組成，當像素較低時，無法真實表征顆粒的真實形狀，需要將這些像素點及噪聲去除。因此，PCAS 識別程序中設定了可識別分析的最小孔隙面積，即設定一個面積識別下限，小于該值的顆粒面積都識別為雜點而被去除。在圖2（b）可以看到二值圖像中呈現為白色的細小雜點，設置最小孔隙面積，可以消除這些雜點（圖2c），有效識別出顆粒。

PCAS 具有手動編輯的功能，對于顆粒不完整或較多明顯雜點，可以通過界面的“Edit binary image”功能，手動修復，歸并顆粒，區分出顆粒邊界，減小系統誤差。

2 PCAS 最優參數探索

以貴州納雍普灑村崩塌-碎屑流為例，運用PCAS對其粒度進行分析，結合彭雙麒等[2]的實測結果，研究PCAS 在崩塌碎屑流顆粒分析的應用，以及閾值、孔喉封閉半徑和最小孔隙面積的參數選取。

2.1 研究區概況

彭雙麒等[2]通過傳統粒徑統計方法，沿納雍崩塌中部區域的主滑方向，將堆積區分為A1—A13 共13個區，每個區長100 m、寬50 m（圖3）。通過生成地表數字模型（DSM），人工觀測統計，得出每個區域顆粒粒徑及對應數量。為了方便測量，粒徑設定為顆粒對角線最大長度取值。

圖3 崩塌堆積分區示意圖[2]Fig.3 Partition diagram of the collapse accumulation body[2]

分別統計出每個區域0～1 m、1～2 m、2～3 m、3～4 m、4～5 m、5～6 m、6～7 m、7～8 m、8～9 m、9～10 m、10～17 m的粒徑數量（識別最大粒徑16.6 m），取值區間左開右閉。根據統計結果，繪制成累計級配曲線（圖4），縱軸為小于某粒徑顆粒的數量百分比，橫軸為各粒徑。

圖4 分區顆粒粒徑累計曲線[2]Fig.4 Aggregate curve of particle size in subsection[2]

為了更直觀顯示A1—A13 每個統計區域內大粒徑與小粒徑含量之間的關系，分別作出各區域0～2 m、2～4 m、4～6 m、大于6 m顆粒所占比例柱狀圖（圖5）。

圖5 各區域不同粒徑等級所占比例Fig.5 Proportion of different grain sizes in different regions

2.2 PCAS 參數討論

運用PCAS 對每個區域顆粒進行定量分析。不斷調整閾值、孔喉封閉半徑、最小孔隙面積，分別統計不同參數下對應的顆粒分布情況，找到合適的參數取值，分析不同參數條件下對研究結果的影響，驗證PCAS 軟件的準確性。

（1）PCAS 參數選取

導入圖像，不斷調整閾值，對13個區域用PCAS軟件進行二值化處理。通過實際分析本次納雍普灑村崩塌的圖像數據，得出閾值在170時可以清晰地區分出顆粒與孔隙（圖6）。

探究孔喉封閉半徑（r）、最小孔隙面積（S0）的取值對PCAS顆粒識別的影響。不斷調整孔喉封閉半徑、最小孔隙面積合適的取值，并進行相應的組合。如圖7所示，當r/S0=0.5/5（這里為像素，通過與實際半徑、面積換算，r/S0=0.5/5=0.05 m/0.05 m2，即r單位為0.1 m，S0單位為0.01 m2）時，許多雜點及不需要的微小顆粒被識別，且顆粒與顆粒間不能有效分隔開；當r/S0=4/40時，較大顆粒也小于最小孔隙面積，不能識別出來；對比圖2（c），當r/S0=1/10時，識別所得顆粒分布結果更能反映堆積物顆粒各粒徑分布情況。因此，選取孔喉封閉半徑r為1，2，3（單位為0.1 m，下同），最小孔隙面積S0為10，20，30，40（單位為0.01 m2，下同），互相組合（表1），找到合適的參數取值（各參數均以像素為單位）。

（2）PCAS 粒徑統計結果

圖6 PCAS 在不同閾值下的二值化結果Fig.6 Binarization results of PCAS at different thresholds

圖7 PCAS 在r/S0=0.5/5 與r/S0=4/40時的顆粒識別結果Fig.7 PCAS particle recognition results at r/S0=0.5/5 and r/S0=4/40

表1 三個參數的組合情況Table1 Combinations of the three parameters

由各個圖像中顆粒的幾何形態數據，統計出每個區域長軸在0～1 m、1～2 m、2～3 m、3～4 m、4～5 m、5～6 m、6～7 m、7～8 m、8～9 m、9～10 m、＞10 m的粒徑數量（識別最大粒徑11.8 m），取值區間為左開右閉，繪制出顆粒粒徑累計曲線。并分別作出各區0～2 m、2～4 m、4～6 m、＞6 m 粒徑所占比例柱狀圖。注意，PCAS 識別的是像素大小，每個區域圖片像素為465×976，而實際尺寸為50 m ×100 m，由式（1）和式（2）換算，為方便統計，取1 像素等于0.1 m。

分別根據表1中的參數組合繪制顆粒粒徑累曲線（圖8）和各區域不同粒徑等級所占比例柱狀圖（圖9），圖9中粒徑百分比小于1%時沒有標記。

由圖8可知：

（1）A11、A12、A13 圖像中粒徑小于2 m的百分含量明顯高于其他區域，說明運動距離越遠，顆粒碰撞破碎越充分，小顆粒含量比重越大，分選越好。在閾值為170，孔喉封閉半徑為2時，A1 圖像中粒徑小于2 m的百分含量也明顯高于其他區域，曲線較陡，顆粒粒徑級配差于其他區域。

（2）當孔喉封閉半徑為1，2時，隨著最小孔隙面積的變大，顆粒粒徑小于2 m的百分比明顯下降，但總的趨勢不變；表明PCAS 在最小孔隙面積增大時，能去除較大的雜點，區分出有效顆粒，總的百分比就降低。但當孔喉封閉半徑為3時，隨著孔隙面積的變大，大粒徑百分比無明顯變化，小顆粒粒徑百分比變小，小于1 m的粒徑百分比大多在10%以下（除A11、A12、A13）；這表明當孔喉封閉半徑增加到一定程度時，許多顆粒被進一步分割成更微小顆粒，當最小孔隙面積較小時（這里為10），大量微細顆粒就被當作雜點去除，此時即使最小孔隙面積變大，粒徑百分比也無明顯變化。

圖8 不同參數取值下各區域顆粒粒徑累計曲線Fig.8 Accumulation curve of particle size in each region under different parameters

（3）以N2、N6、N10為例，當最小孔隙面積一定時，隨著孔喉封閉半徑的增大，顆粒粒徑小于2 m的百分比明顯下降，且級配變差。N10 條件下，每個區域粒徑小于2 m的百分比分布在50%～90%，而N2、N6 小于2 m的百分比則在80%～95%；這表明隨著孔喉封閉半徑增大，形狀不規則的大顆粒被腐蝕，分割為小顆粒，當小于最小孔隙面積時，被作為雜點去除，所占比重減小。

（4）小于3 m的顆粒百分比在每個區域的占比都在80%以上，而圖4小于3 m的為70%以上，說明PCAS 相比人工計數，顆粒計數更精細，識別更完整，但總的效果相差不大，粒徑級配曲線與實際情況相符，驗證了PCAS 在碎屑流堆積物顆粒識別上的可靠性。

由圖9可知：

（1）由圖N1、N2、N5、N6可得，當孔喉封閉半徑為1，2，最小孔隙面積為10，20時，各區域0～2 m 小顆粒所占百分比都在80%以上，大于4 m顆粒所占百分比都在4%以內；表明當孔喉封閉半徑與最小孔隙面積較小時，各區域較多的小顆粒與雜點被識別，且不能被錯誤識別為單個顆粒。對比N3、N4、N7、N8，當孔喉封閉半徑或最小孔隙面積增大時，0～2 m 小顆粒占比降低，大于4 m的大顆粒比重增加。

（2）由圖N9、N10、N11、N12可得，當孔喉封閉半徑為3時，各區域小于2 m的顆粒占比都具有先減小后增大的趨勢，拐點在A3區；有3個明顯波峰，極點在A2、A9、A12區，有三個明顯波谷，極點在A3、A8、A10區。各區域4～6 m顆粒占比都在2%～4%附近，小于6 m的顆粒占比具有2個波谷，極點在A4 與A8區。各粒徑占比分布情況與圖5對比，基本吻合，接近真實值。

（3）對比圖N9、N10、N11、N12 與圖5，大于6 m顆粒占比在A3、A4、A5 三個區域與實際情況相符。圖5中，0～2 m顆粒占比在A3、A4、A5區呈現波谷形態，但在PCAS 粒徑識別占比中，A1～A10區0～2 m顆粒占比都在52%～69%，表明在PCAS顆粒識別系統中，小粒徑顆粒識別效率高，各區占比較接近。

2.3 PCAS 優劣及可行性討論

根據PCAS 在閾值為170、孔喉封閉半徑為3、最小孔隙面積為30 三個參數下取得的各區粒徑分布結果，結合人工統計數據，作各區粒徑塊數對比柱狀圖（圖10）。結合圖8、圖9可得如下結果：

圖10 PCAS 統計方法在T/r/S0=170/3/30時與人工統計各區顆粒塊數對比圖Fig.10 Comparison chart of the PCAS statistical method at T/r/S0=170/3/30 and artificial statistics

2.3.1 優勢

(1) 高效性：PCAS 統計的細顆粒數量相比人工統計的多，但顆粒數量變化趨勢基本一致，各區數量變化情況基本對應；0～2 m 粒徑塊數在A2、A4、A7 附近出現極值，2～4 m 粒徑數量在A4、A7 附近出現極值。表明PCAS可以有效識別肉眼難以計數的細小顆粒，可以極大提高對小粒徑顆粒（小于4 m）的識別效率，相比人工細小顆粒統計更高效。

(2) 準確性：當顆粒粒徑大于4 m時，人工統計與PCAS 統計結果基本吻合（A8區異常）。人工對大顆粒的識別相比小顆粒準確性較高，這也證實了PCAS對堆積體識別的準確性，人工可以識別的，PCAS 也能自動準確識別。

(3) 規律性：從各區顆粒數量變化趨勢可以看出，PCAS 統計顆粒數量變化呈相當好的規律，曲線圓潤平滑，“波峰波谷”的數量極值與貴州納雍等碎屑流崩塌實際情況對應較好[17-18]，且符合碎屑流堆積體的一般規律[19]。而人工統計由于肉眼的局限性，在顆粒數量、顆粒大小上總會存在判斷失誤，與實際情況會有一定的差異。

(4)可操作性：PCAS顆粒識別系統操作簡單，只需設定閾值（T）、孔喉封閉半徑（r）、最小孔隙面積

（S0）三個參數，即可得到堆積體顆粒分布情況，可節省大量人力物力。

2.3.2 劣勢

PCAS 劣勢主要體現在三個參數自身的特性上，分析如下：

(1)A2、A8區顆粒數量出現奇點；圖10（a）、（b）中，A2區人工統計粒徑小于4 m顆粒數量比PCAS 計數高，主要是因為當顆粒小于PCAS 設置的最小孔隙面積時，就被自動刪除，導致A2區粒徑小于4 m的顆粒數量較低；圖10（c）、（d）中，A8區人工統計粒徑4 m以上顆粒數量較高，主要是因為PCAS 設置的孔喉封閉半徑過大，許多形狀不規則的大顆粒被錯誤分割成多個小顆粒，導致大顆粒數量變少，并且這些誤分割成的小顆粒一部分大于最小孔隙面積，被識別成單個顆粒，導致A8區2～4 m 粒徑顆粒數量較人工統計的多。

(2)PCAS 圖像二值化，是根據圖像中孔隙與顆粒的灰度值不同加以區分的，主要取決于所得堆積體圖像。如果堆積體圖像的分辨率、對比度、飽和度以及拍攝角度等不同時，選取的閾值可能會有差異，圖像二值化結果也可能不同。

2.3.3 可行性

(1) 綜合圖8～10，當PCAS顆粒識別系統中采用不同的參數取值時，對分析崩塌碎屑流顆粒分布情況會有一定的影響。當所得圖像與本文納雍崩塌圖像相似，閾值170、孔喉封閉半徑3、最小孔隙面積30為最優參數選擇，各區不同粒徑數量與人工統計所得結果基本對應，粒徑累計曲線與各區域不同粒徑等級所占比例最接近實際情況；反之，可根據本文所介紹的PCAS 識別方法、步驟以及參數選取的判定依據，依據具體情況以“170、3、30”為基準進行合理的參數調整。

(2) PCAS顆粒識別系統已在貴州納雍崩塌[10]和金沙江白格滑坡[11]中得到應用，并且效果顯著。得到滑坡堆積體圖像后，只需設定三個參數，即可對堆積體進行顆粒識別分析，操作簡單，并且PCAS 具有手動編輯的功能，可以通過界面的“Edit binary image”功能，根據實際情況進行顆粒的調整。在新技術不斷發展的今天，PCAS 不失為崩塌碎屑流粒徑研究的一個高效、可行的途徑，具有應用價值。

3 結論

（1）PCAS 能自動準確地識別碎屑流堆積物顆粒與孔隙，具有高效、準確、可操作性強等特點，可用于堆積物顆粒數字圖像的識別、量化、分析；確定三個關鍵參數：閾值（T）、孔喉封閉半徑（r）、最小孔隙面積（S0）的合適取值后，即可得出粒徑分布情況，節省大量人力物力。

（2）針對顆粒粒徑較大、二值化后孔隙與顆粒區分不明顯、噪點較多的圖像，采用較大的r/S0值更能反應顆粒分布宏觀情況；反之，宜適當減小r/S0進行分析。當所得圖像與本文納雍崩塌圖像相似，閾值170、孔喉封閉半徑3、最小孔隙面積30為最優的參數選擇，可得到有效的碎屑流顆粒分布情況；反之，可以“170、3、30”為基準進行合理的參數調整。

（3）PCAS 具有較高的可行性，所得貴州納雍崩塌堆積物粒徑識別和統計結果與實測值接近，粒徑占比、分布規律基本吻合，即堆積體小粒徑占比較多，0～2 m顆粒粒徑各區占比都在50%以上，且運動距離越遠小粒徑比例越大，各區“波峰波谷”變化趨勢基本對應。