基于凍融循環的土壤物理狀態的自動判別

韓巧玲，趙 玥，姚立紅

(北京林業大學 工學院，北京 100083)

基于凍融循環的土壤物理狀態的自動判別

韓巧玲，趙 玥，姚立紅*

(北京林業大學 工學院，北京 100083)

以東北典型黑土區土壤為研究對象，采用CT掃描技術與圖像處理相結合的方法，通過灰度共生矩陣和主成分分析法提取圖像特征，計算測試圖像特征向量與訓練圖像特征向量間的歐氏距離，以此為依據，實現對經歷不同凍融循環次數土壤的自動判別。研究結果表明：面向土壤CT圖像數據庫，基于灰度共生矩陣和主成分分析提取的圖像特征，均能實現對土壤的自動判別，但灰度共生矩陣法的判別正確率要高于主成分分析法。

土壤CT圖像；灰度共生矩陣；主成分分析；歐氏距離；判別正確率

土壤物理主要包括土壤固、液、氣三相體系所產生的各種物理現象和過程。土壤物理性質的變化制約著土壤肥力水平，影響土壤養分的保持與運移，是合理耕作和灌排的重要依據[1-2]。對土壤不同物理狀態進行判別分析，不僅有利于加深對土壤內部結構的研究，也可為耕作區植物的健康生長提供科學依據。

在我國東北黑土區，凍融循環是改變土壤物理性質的關鍵因素。在凍融過程中，土壤水分的相變與變遷會改變土壤顆粒組成、團聚體的穩定性和土壤導水率等物理因子，從而影響土壤水熱傳導和養分運移等特性[3-4]。齊吉琳等[5]通過對孔隙的量化分析，發現凍融循環能夠改變土壤孔隙的結構；鄧西民等[6]研究發現，凍融循環會降低土壤容重，增加土壤孔隙度和飽和導水率?？傮w來看，現有成果大多研究凍融循環引起的土壤結構的改變，鮮有人研究基于凍融循環的土壤物理狀態的自動判別。現有的判別方法主要由人工實現，對于判斷者有著極為嚴格的要求，不僅耗費大量的人力物力，而且結果也會因為主觀判斷存在一定的偏差。為了克服這一缺陷，本研究基于土壤CT圖像，運用數字圖像處理技術，完成對土壤結構特征的提取與描述，以期實現對基于凍融循環的土壤不同物理狀態的自動判別，以便于及時了解土壤狀態的變化，為農林業的健康發展和生態環境保護提供基礎數據。

1 材料與方法

1.1 土壤樣本制備

所用土壤選自黑龍江省克山農場，土壤類型以黏化濕潤均腐土為主，采用自制內徑和高分別為10 cm的有機玻璃管于0～40 cm層深的侵蝕溝壁進行原狀土取樣，每層重復取樣3次。將樣本進行飽和水處理后[7]，共得到原狀飽和水(3個)、原狀非飽和水(3個)、填充飽和水(3個)、填充非飽和水(3個)等共計12個樣本，以供后續CT掃描所用。

1.2 土壤CT圖像獲取

選用黑龍江省中醫藥大學的Lightspeed16排螺旋CT掃描儀對原狀樣本進行掃描處理。CT掃描儀選用參數分別為：電壓140 kV，電流60 mA，掃描間隔1.5 s，掃描層厚2 mm，窗寬(顯示CT圖像是所選用的CT值范圍)和窗位(窗寬上、下限CT值的平均數) 均為1 300，對每個土柱樣品進行7次橫斷面掃描。

本研究的單個掃描樣本分別經歷0、1、3、6、9次凍融循環，每次凍融循環后進行1次掃描，每次掃描可得7張土壤CT圖像，故單個樣本凍融循環后共需5次掃描，得35幅土壤CT圖像。本研究合計獲得420幅土壤CT斷層掃描圖像，用于建立土壤CT圖像數據庫，作為測試圖像庫。

1.3 土壤CT圖像特征提取

1.3.1 土壤CT圖像處理

原始土壤CT圖像以醫療圖像的格式存儲，并不便于計算機的后續處理。同時，由于CT機器硬件技術缺陷的存在，使得圖像邊界會產生畸變，而且，數字圖像在獲取和傳輸的過程中也會受到噪聲影響，這就決定了基于原始圖像的特征提取可能會影響試驗結果的精確性[8-9]。因此，須針對土壤CT圖像首先進行一些處理。

如圖1-a所示的原始圖像中包含了與土壤分析無關的CT機器掃描參數的基本信息，因此，本試驗基于Matlab編程軟件，采用最大內切正方形的方法，剪裁得到有效土壤CT圖像，并保存為*.bmp格式文件。通過自適應中值濾波算法消除噪聲對后續特征提取的干擾，以此來抑制無用信息，突出有效信息，得到用于結果分析的試驗用圖像，如圖1-b所示。

a. 原始土壤CT圖像 b. 有效CT圖像a. The original soil CT images b. Effective CT images圖1 土壤CT圖像預處理Fig.1 Preprocessing of soil CT image

1.3.2 灰度共生矩陣法

紋理是土壤CT圖像中每個像素的CT值或灰度值在空間中以某種方式變化而產生的，是圖像客觀存在的一種屬性，反映了圖像中灰度的空間相關特性[10]。灰度共生矩陣法是通過研究灰度的空間相關特性來描述圖像紋理變化的方法，反映圖像灰度關于方向、相鄰間隔、變化幅度的綜合信息，是分析圖像的局部模式和排列規則的基礎。

灰度共生矩陣定義為：灰度值為i的像素沿θ方向，到達與其相隔距離為d的像素點(灰度值為j)的概率。設g(x，y)為圖像中位于(x，y)像素點的灰度值，則其數學表達式為：

P(i，j，d，θ)=#{(x1，y1)，(x2，y2)|f(x1，y1)=i，f(x2，y2)=j，|(x1，y1)-(x2，y2)|=d，∠(x1，y1)，(x2，y2)=θ}。

(1)

式(1)中，#表示根據式(1)的條件，圖像中位于(x1，y1)灰度值是i，位于(x2，y2)灰度值是j的像素對個數。通常，d={1，2，3，4}，θ={0°，45°，90°，135°}。

為了便于二次統計特征的提取，需對灰度共生矩陣進行歸一化處理：

(2)

Haralick等[11]根據紋理信息的特點，提出了用于分析灰度共生矩陣的多個特征，用于挖掘圖像深層信息。常用的特征參數簡述如下。

(1)角二階矩

(3)

式(3)中，k為灰度級[式(4)～(6)同]。角二階矩也稱為能量，是灰度共生矩陣各個元素值的平方和，反映了圖像灰度分布均勻程度和紋理粗細度。

(2)對比度

(4)

對比度，也稱主對角線慣性矩，代表檢測像素灰度值與其領域像素灰度值的亮度反差，反映了圖像的清晰度和紋理溝紋深淺的程度。

(3)逆差矩

(5)

逆差矩用于度量圖像紋理局部變化的多少，反映圖像紋理的同質性，其值大說明圖像紋理的不同區域間缺少變化，局部非常均勻。

(4)熵

(6)

熵是圖像所具有的信息量的隨機性度量，反映了圖像中紋理的非均勻程度或復雜程度。

1.3.3 主成分分析法

主成分分析是通過線性變換將多個指標簡化為少數指標的統計分析方法，用少數互不相關的指標盡可能多地表征圖像信息，以保證原始信息損失量最小[12]。其基本原理是：將原始具有相關性的多個指標X1，X2，…，Xp(p個指標)，重新組合成一組互不相關的綜合指標Fm來代替原來指標。

設X=(X1，X2，…，Xp)T是p維隨機向量，其線性變化為：

(7)

采用新變量F1來表征原始p個向量，F1所包含的信息應盡可能多地反映所有變量信息，則其為第一主成分。若第一主成分不足以代表足夠多的原始信息，則引入第二主成分，以此類推。其中，新變量F1，F2，…，Fm要求滿足以下條件：

(1)變量所提取的信息量由方差度量，即要求方差Var(Fm)盡量大；

(2)Fi與Fj互不相關，即COV(Fi，Fj)=0。

主成分分析的主要目的是簡化數據量，因此，所提取的主成分個數m應小于p。m的值依據各主成分的累計方差貢獻率P(m)來確定。其數學表達式為：

(8)

式(8)中：λ為各個主成分對應的特征值；k為全部的主成分個數；i為選定主成分個數。

一般地，當累積貢獻率大于85%時，就認為足夠反映原來變量的信息[13]，也可根據實際應用的需求自主選擇累積貢獻率的范圍。

1.4 識別正確率的計算

1.4.1 驗證方法確定

由于土壤CT掃描圖像數據庫中單個樣本單次凍融循環的圖像為7幅，數據量相對較少，需選用適合于小樣本的驗證方法。留一驗證法的思想是選擇一幅作為測試圖像，剩余的作為訓練圖像，以此循環，直至所有圖像都參與測試與訓練。該方法對圖像的利用率高，能提高實驗結果的精確度，故本研究選用該方法用于對土壤CT圖像的判別與分類[14]。

1.4.2 歐氏距離計算

歐氏距離指的是在m維空間中2個點之間的真實距離，通過歐氏距離的數值可準確定義2個特征向量間的相似度。

1.4.3 識別正確率計算

若測試圖像對于本組內歐氏距離較其他組小，則為判別正確，反之則為錯誤，最終正確判斷個數與總判斷數的比值即為識別正確率。

2 結果與分析

2.1 特征提取

2.1.1 灰度共生矩陣特征

灰度共生矩陣法的主要變量是方向和特征個數。選定0°、45°、90°、135°，提取能量、熵、對比度、逆差矩共4個特征對土壤CT圖像紋理信息進行描述。由于篇幅有限，隨機選取經歷不同凍融次數的同一土壤CT圖像進行部分特征數據的展示，結果如圖2所示。

能量和逆差距主要反映土壤紋理信息的同質性。由圖2可知，不同方向的能量和逆差距先隨著凍融次數的增加而增加，當凍融次數達到6次后趨于減小，說明土壤內部結構在前6次凍融過程中變化比較強烈，之后變化減弱，內部物質結構漸趨均勻。兩個方向的熵和對比度數值較小，且呈現相反的變化趨勢，這說明當土壤圖像紋理信息趨于隨機時，土壤內部物質之間具有一定相似性。4個特征綜合起來表征了土壤內部物質隨著凍融次數的變化，在第9次凍融循環后內部物質趨于均勻化。

圖2 灰度共生矩陣特征值Fig.2 The eigenvalue of gray symbiotic matrix

2.1.2 主成分分析特征

采用主成分分析方法，對不同凍融循環次數的土壤CT圖像進行主要特征分量的提取。選用4種凍融循環次數下的圖像，分別記為A1、A2、A3、A4，設定各主成分的累計方差貢獻率大于0.85，最多選取16維主成分特征。由于篇幅有限，選取4幅圖像的10維特征的貢獻率進行展示，如表1所示，ai表示提取的第i維特征。

由于圖像所包含的信息不同，圖像提取的每維主成分貢獻率亦不相同，累計方差貢獻率大于0.85的主成分特征維數也不相同，最大有16維，最小12維。為保證圖像識別的正確率，所有圖像均選取16維主成分特征參與后續的計算分析。

2.2 判別正確率

由于篇幅有限，僅舉一例用于說明判別過程。首先，將原狀飽和水樣本經歷1、3、6、9次凍融循環的掃描圖像分別編組，標記為Y1、Y3、Y6、Y9，下標中的數字代表其所經歷的凍融次數；然后，選擇Y3組中的7幅圖像(即經歷3次凍融的原狀土圖像)，依次命名為I1～I7，作為測試樣本，分別測算其每一幅圖片的灰度共生矩陣特征向量與Y1、Y3、Y6、Y9組均值特征向量之間的歐式距離，與哪一組的歐式距離數值最小，則自動判別其屬于該組。如表2所示，屬于Y3組的7幅圖像中，除I4被自動判別為Y9組(即經歷9次凍融)，判別錯誤外，其他判斷均正確，則該測試樣本的判別正確率為(6÷7)×100%=85.7%。

基于本研究所用土壤CT圖像數據庫，對全部樣本進行自動判別，最后測定得出采用灰度共生矩陣法自動判別的正確率是81.4%，主成分分析法自動判別的正確率是61.0%。

表1 主成分特征貢獻率

Table 1 Contribution rate of principal component

圖像Imagea1a2a3a4a5a6a7a8a9a10A101915017870092300519004630043900419003290032400276A302647019800063100483004060038600312002780027200242A602732016290088900757006630045400336003160025800221A902922016830092900651005160036400316002720025500248

表2 樣本歐氏距離

Table 2 The Euclidean distance of samples

組別GroupI1I2I3I4I5I6I7Y1997946699885730921611140599198895257Y360742?50340?56900?6618857590?69309?60481?Y6107891675129323610452311143499628110759Y960932528565772160670?679397030660547

每列中加*的表示所在列圖像與該組的歐氏距離最小。

* indicated the smallest Euclidean distance within the image tested in the column and the corresponding groups.

3 討論

本研究以東北典型黑土區土壤為研究對象，合計獲得420幅經歷不同凍融次數的土壤CT斷層掃描圖像，建立土壤CT圖像數據庫，通過灰度共生矩陣和主成分分析法提取圖像特征，計算測試圖像特征向量與訓練圖像特征向量間的歐氏距離，以此為依據，實現對經歷不同凍融次數土壤的自動判別。對比2種特征提取方法對土壤的自動判別正確率可知，灰度共生矩陣的判別正確率較高，其判別效果優于主成分分析法。這主要是因為：(1)空間分辨率的限制使得土壤CT圖像成像出現部分容積效應，即圖像某點的像素值由鄰域內像素值的均值來體現，即像素點的灰度值會受到鄰域空間內灰度值的影響。由于灰度共生矩陣法正是通過灰度值的空間特性來表征圖像的紋理信息，能夠較好地表征出像素值之間的聯系，因此，用于土壤CT掃描圖像的特征提取時能夠獲得較高的正確率；(2)本研究采用灰度共生矩陣法提取出了四維特征向量，通過將其融合豐富了特征向量的信息，使其能夠全面反映土壤CT圖像的內部結構信息，有助于提高其識別正確率。這說明，多種特征提取方法的融合是后期進行土壤物理狀態判別的有效技術手段。而主成分分析法是用較少的數據量來盡可能地表征最多的信息，且各指標之間不具有相關性，這樣做雖然避免了數據的重復使用，可以極大地減少計算量，具有一定的運算優勢，但也丟失了土壤CT圖像中的部分細節信息。土壤具有多重復雜的結構，表現在圖像上的細微的差異，都會對應土壤內部物質結構的變化，因此，在實際應用中其辨別效率較低。

本研究表明，灰度共生矩陣法能夠表征灰度空間變化關系，其特征具有不同的土壤物理含義，通過建立該特征與土壤凍融次數、土壤含水率等的關系，有利于圖像深層次信息的挖掘，不僅為研究土壤內部結構提供了有效方法，而且還能為后期研究凍融循環和含水率對土壤內部結構的影響奠定基礎。

[1] HORN R， TAUBNER H， WUTTKE M， et al. Soil physical properties related to soil structure[J].Soil&TillageResearch， 1994， 30(2/3/4):187-216.

[2] KONRAD J M. Physical processes during freeze-thaw cycles in clayey silts[J].ColdRegionsScience&Technology， 1989， 16(3):291-303.

[3] GATTO L W. Soil freeze-thaw-induced changes to a simulated rill: potential impacts on soil erosion[J].Geomorphology， 2000， 32(1):147-160.

[4] 安樂生， 劉春， 廖凱華. 一種改進的土壤水分特征曲線模型及其驗證[J]. 浙江農業學報， 2015， 27(5):837-840. AN L S， LIU C， LIAO K H. An improved model for prediction of soil water retention properties and its validation[J].ActaAgriculturaeZhejiangensis， 2015， 27(5): 837-840. (in Chinese with English abstract)

[5] 齊吉琳， 張建明， 朱元林. 凍融作用對土結構性影響的土力學意義[J]. 巖石力學與工程學報， 2003， 22(增刊2):2690-2694. QI J L， ZHANG J M， ZHU Y L. Influence of freezing-thawing on soil structure and its soil mechanics significance[J].ChineseJournalofRockMechanicsandEngineering， 2003， 22(Supp. 2): 2690-2694. (in Chinese with English abstract)

[6] 鄧西民， 王堅， 朱文珊， 等. 凍融作用對犁底層土壤物理性狀的影響[J]. 科學通報， 1998， 43(23):2538-2541. DENG X M， WANG J， ZHU W S， et al. Influence of the freeze-thaw action on physical properties in plow pan soil[J].ChineseScienceBulletin， 1998， 43(23): 2538-2541. (in Chinese)

[7] 夏祥友， 王恩姮， 楊小燕， 等. 模擬凍融循環對黑土黏化層孔隙結構的影響[J]. 北京林業大學學報， 2015， 37(6):70-76. XIA X Y， WANG E H， YANG X Y， et al. Pore characteristics of mollisol argillic horizon under simulated freeze-thaw cycles[J].JournalofBeijingForestryUniversity， 2015， 37(6): 70-76. (in Chinese with English abstract)

[8] 王艷麗， 程展林. CT掃描技術在我國土工試驗中的應用研究進展[J]. 地震工程學報， 2015， 37(增刊1):35-39. WANG Y L， CHENG Z L. Progress in the application of CT scanning technology in Chinese soil tests[J].ChinaEarthquakeEngineeringJournal， 2015， 37(Supp. 1): 35-39. (in Chinese with English abstract)

[9] ZELELEW H M， ALMUNTASHRI A， AGAIAN S， et al. An improved image processing technique for asphalt concrete X-ray CT images[J].RoadMaterials&PavementDesign， 2013， 14(2):341-359.

[10] 劉志剛， 李德仁， 秦前清， 等. 支持向量機在多類分類問題中的推廣[J]. 計算機工程與應用， 2004， 40(7):10-13. LIU Z G， LI D R， QIN Q Q， et al. An analytical overview of methods for multi-category support vector machines[J].ComputerEngineering&Applications， 2004， 40(7):10-13. (in Chinese with English abstract)

[11] HARALICK R M， SHANMUGAN K， DINSTEN I. Textural features for image classification[J].IEEETransactionsonSystem，ManandCybernetics， 1973， 3(6):610-621.

[12] 高惠璇. 應用多元統計分析[M]. 北京：北京大學出版社，2005.

[13] PATRAS A， BRUNTON N P， DOWNEY G， et al. Application of principal component and hierarchical cluster analysis to classify fruits and vegetables commonly consumed in Ireland based oninvitroantioxidant activity[J].JournalofFoodComposition&Analysis， 2011， 24(2):250-256.

[14] CAWLEY G C， TALBOT N L C. Efficient leave-one-out cross-validation of kernel fisher discriminant classifiers[J].PatternRecognition， 2003， 36(11):2585-2592.

(責任編輯 高 峻)

Automatic identification of different soil physical state caused by freeze-thaw

HAN Qiaoling， ZHAO Yue， YAO Lihong*

(SchoolofTechnology，BeijingForestryUniversity，Beijing100083，China)

In the present study， typical black soil in northeastern China was selected as the test object， and the simulated image processing was adopted combined with computerized tomography (CT) scanning. With the extracted image feature by gray level co-occurrence matrix and principal component analysis (PCA)， the Euclidean distance between the feature vector of test image and verify image was calculated， which built the basis for automatic discrimination of different soil physical states caused by freeze-thaw. It was shown that it could realize automatic identification of different soil physical state by image features extracted by either gray level co-occurrence matrix or PCA. And the identification accuracy of gray-level co-occurrence matrix method was higher than that of principal component analysis method for the same soil CT tomography image database.

soil computed tomography images; gray level co-occurrence matrix; principal component analysis; Euclidean distance; identification accuracy

http://www.zjnyxb.cn

10.3969/j.issn.1004-1524.2017.07.18

2017-01-13

國家自然科學基金項目(41501283)；中央高?；究蒲袠I務費專項資金(BLX2015-36)

韓巧玲(1990—)，女，河南安陽人，碩士研究生，主要從事圖像處理與模式識別方向研究。E-mail: hanqiaoling0@163.com

*通信作者，姚立紅，E-mail: yaolihong@bjfu.edu.cn

S1

A

1004-1524(2017)07-1189-06

浙江農業學報ActaAgriculturaeZhejiangensis， 2017，29(7): 1189-1194

韓巧玲，趙玥，姚立紅. 基于凍融循環的土壤物理狀態的自動判別[J]. 浙江農業學報，2017，29(7)： 1189-1194.