基于大數(shù)據(jù)挖掘的橡膠園耕地地力評價

彭金蓮，唐朝勝，蔣菊生，趙春梅

(1. 海南大學 計算機與網(wǎng)絡空間安全學院， ?？?570228；2. 海南省農墾科學院， ?？?570206；3.中國熱帶農業(yè)科學院 橡膠研究所 ?？?571101)

我國耕地地力評價主要方法有模糊數(shù)學法、層次分析法、指數(shù)和法等，我國第2次土壤普查多采用指數(shù)和法評定耕地地力等級。近年來，陳桂芬等[1]應用K-means聚類方法、Johnson粗糙集屬性約簡算法與C4.5決策樹算法優(yōu)化地力評價。蔡麗霞等[2]應用C4.5決策樹、K-means和DBSCAN聚類算法，對土壤養(yǎng)分數(shù)據(jù)進行對比分析，依靠大量土壤肥力狀況以預測未來土壤肥力的變化趨勢。閆一凡[3]認為評價指標篩選目前尚無國際通用的標準，中國已建立了包括7類共64項指標的公用指標體系;當前主流的地力評價方法(如中國農業(yè)部推薦的特爾斐—層次分析法等)仍存在人為隨意性等不足，推薦將不同評價方法(特別是特爾斐—層次分析法和分類與回歸樹模型法的組合)組合起來用于實際工作，可取長補短，進一步提高結果的可靠性。建議建立起統(tǒng)一的地力等級定義，引入分類與回歸樹、人工智能、機器學習等新方法，并將其與現(xiàn)有評價方法(如特爾斐法、多元回歸等)集成，充分挖掘數(shù)據(jù)、實現(xiàn)優(yōu)勢互補，認為是未來耕地地力評價研究的重點和發(fā)展趨勢[4-7]。海南農墾測土配方施肥項目從2006年開始試點，歷時10年，實施范圍為海南農墾原有的92個農場(分10個片區(qū))，面積大約2萬hm2，依據(jù)測土配方調查數(shù)據(jù)進行了橡膠園耕地地力評價。采用了專家打分法和特爾菲法，10個片區(qū)分別建立了耕地地力評價指標體系，10個片區(qū)采用的耕地地力評價指標不盡相同，耕地地力評價指標的選擇與重要度由專家打分確定，人為主觀影響明顯，無法對海南墾區(qū)膠園耕地地力進行統(tǒng)一評價。筆者對海南墾區(qū)膠園大數(shù)據(jù)進行挖掘，旨在篩選出對海南膠園耕地類型有影響的地力指標，以便客觀評價海南墾區(qū)膠園的耕地地力。

1 研究對象與技術方法

1.1 研究對象本研究通過整合海南農墾測土配方施肥10個片區(qū)(項目場)的屬性數(shù)據(jù)庫、空間數(shù)據(jù)庫和測土配方施肥數(shù)據(jù)庫，創(chuàng)建全墾區(qū)統(tǒng)一的耕地地力評價數(shù)據(jù)庫；對海南墾區(qū)測土配方施肥原始數(shù)據(jù)進行數(shù)據(jù)探索，通過數(shù)據(jù)分析，探索土壤內在的分類規(guī)律，篩選出與墾區(qū)耕地地力相關的主要因素，以確定海南墾區(qū)的耕地地力指標；使用4種聚類模型，利用輪廓系統(tǒng)進行適用性評估，最終確定一種聚類算法作為耕地地力分類模型來劃分膠園耕地類型，從而實現(xiàn)對海南橡膠園耕地地力進行統(tǒng)一的綜合評價。

1.2 技術方法使用IBM公司開發(fā)的企業(yè)級大型數(shù)據(jù)挖掘工具SPSS Modeler進行數(shù)據(jù)挖掘，技術路線參照跨行業(yè)數(shù)據(jù)挖掘標準流程 CRISP-DM 模型設計而成，此數(shù)據(jù)挖掘標準流程由問題理解、數(shù)據(jù)理解、數(shù)據(jù)準備、建立模型、模型評估和成果發(fā)布等6個階段構成。

1.2.1 分類型指標的數(shù)據(jù)預處理通過對原始數(shù)據(jù)的審核，發(fā)現(xiàn)能夠作為地力評價潛在的分類指標為:地貌類型、地形部位、坡向、成土母質和質地等。但由于原始數(shù)據(jù)時間跨度大、涉及人員多，許多數(shù)據(jù)沒有嚴格按照操作規(guī)程填寫，導致各指標的類別數(shù)超出規(guī)程規(guī)定，產生了地貌類型71類，地形部位95類，坡向71類，成土母質41類，質地64類，所以數(shù)據(jù)質量較差，因此,必須按照農業(yè)部測土配方施肥野外調查操作規(guī)程要求調整。其中，地形部位原來只有其他部位和山腰，已經沒有數(shù)據(jù)挖掘價值。調整前后各指標及其類別值見表1。調整后，地貌類型調整為3種類別，分別是山地、丘陵、平地。坡向調整為9種類別，分別是北、東北、東、東南、南、西南、西、西北和平地等。成土母質調整為8種類別，分別是安山巖、變質巖、花崗巖、淺海沉積物、砂頁巖、石灰?guī)r、玄武巖和云母片巖風化物。質地調整為6種類別，分別是粘土、粘壤土、中壤土、砂壤土、砂土和其他等。分類型地力指標的缺失值處理依據(jù)不同情況進行了不同的處理。地形部位只有4個片區(qū)農場有調查數(shù)據(jù)，缺失值占74%，已經沒有數(shù)據(jù)挖掘價值。地貌類型的缺失值為11.44%，只有1個片區(qū)農場沒有調查此類數(shù)據(jù)，通過C&RT模型進行缺失值插補。質地指標60%是其他類別，已經沒有分類價值。

表1 地力評價各分類指標調查數(shù)據(jù)調整前后類別數(shù)

1.2.2 連續(xù)型地力指標的預處理連續(xù)型地力指標的預處理相對分類型指標來說要復雜很多，包括離群值、極值、缺失值的數(shù)據(jù)預處理，也包括數(shù)據(jù)分布和數(shù)據(jù)分箱處理。連續(xù)型地力指標共有10個，分別是pH值、有機質、全氮、有效磷、速效鉀、北緯、東經、海拔、常年平均降雨量和有效土層厚度。從描述性統(tǒng)計分析，查看均值、標準差、最大值和最小值。一般情況下，若標準差遠遠大于均值，可粗略判定數(shù)據(jù)存在異常值。采用3倍標準差上最大值替換離群值和丟棄極值的方法。 為了提高數(shù)據(jù)挖掘效率，基于對連續(xù)型地力指標進行離散化處理和適應某些模型只能處理分類型變量的要求，減少變量取值數(shù)以實現(xiàn)樣本量的縮減。參考《海南島土系概論》[8]，對海南農墾膠園土壤肥力的養(yǎng)分指標進行分組，如表2所示,分組給出的區(qū)間范圍是大于左邊的區(qū)間數(shù)和小于右邊的區(qū)間數(shù)。在SPSS Modeler中對連續(xù)數(shù)據(jù)分組的方法有固定寬度、分位數(shù)、等級、平均值/方差和最優(yōu)等5種方法，根據(jù)調查樣本的數(shù)據(jù)分布，確定海拔高度和有效土層厚度采用樣本量固定寬度分級，年平均降雨量、東經和北緯采用數(shù)值固定寬度分組，各立地指標分組標準如表3所示，分組給出的區(qū)間范圍是大于左邊的區(qū)間數(shù)和小于右邊的區(qū)間數(shù)。

表2 土壤養(yǎng)分指標分組

表3 立地指標分組標準表

1.2.3 聚類模型算法聚類分析將一批樣本數(shù)據(jù)，在沒有先驗知識的前提下，根據(jù)數(shù)據(jù)的諸多特征，按照其在特征上親疏自動分組，使得組內個體的結構特征具有較大的相似性，組間個體特征相似性較低。聚類與分類不同，聚類分析的輸入數(shù)據(jù)集是一組未標記的對象，也就是說此時輸入的對象還沒有被進行任何分類，聚類的目的是根據(jù)一定的規(guī)則合理地進行分組或聚類，并用顯式或隱式的方法描述不同的類別。由于分析可以采用不同的算法，所以對于相同的數(shù)據(jù)集合可能有不同的劃分。從聚類的原理角度來看，聚類的算法可以分為:劃分聚類、層次聚類、基于密度的聚類和基于網(wǎng)格聚類等方法。本研究分別使用K-Means，Two-Step，Kohonen和TwoStep-AS等4種聚類算法，對海南農墾測土調查數(shù)據(jù)進行聚類分析，最終，根據(jù)輪廓系數(shù)(Silhouette Coefficient)的大小確定TwoStep-AS聚類模型為海南墾區(qū)耕地進行分類。

1.2.3.1 TwoStep-AS聚類模型算法介紹[9]SPSS Modeler中的TwoStep-AS聚類算法，它的優(yōu)勢至少表現(xiàn)在以下幾個方面：1)可同時基于類別變量和連續(xù)變量進行聚類；2)可自動確定最終的分類個數(shù)；3)可處理大型數(shù)據(jù)集。

TwoStep-AS聚類基本步驟如下：

1)預聚類：采用“貫序”方式將樣本粗略劃分成L個子類。開始階段視所有數(shù)據(jù)為一個大類。讀入一個樣本數(shù)據(jù)后，根據(jù)“親疏程度”決定該樣本應派生出一個新類，還是應合并到已有的某個子類中。這個過程反復進行，最終形成L個類。預聚類過程聚類數(shù)目不斷增加。

2)聚類：在預聚類的基礎上，再根據(jù)“親疏程度”決定哪些子類可以合并，最終形成L類。聚類數(shù)目不斷減少的過程，隨著聚類的進行，類內部的差異性將不斷增大。

在距離測量有“歐氏距離”和“對數(shù)似然”兩種算法作為聚類變量相似度的測量形式；如果聚類變量均為數(shù)值型，采用歐氏距離，若為混合型采用對數(shù)似然距離。在聚類準則有貝葉斯(BIC)和Akaik信息準則(AIC)，作為聚類個數(shù)的判斷依據(jù)。同其他統(tǒng)計方法一樣，TwoStep-AS聚類算法也有嚴苛的適用條件，它要求模型中的變量獨立，類別變量是多項式分布，連續(xù)變量符合正態(tài)分布。

注意，與其他的聚類分析算法一樣，所得到的模型可以在一定程度上依賴于訓練數(shù)據(jù)的順序。重新排序數(shù)據(jù)和重建模型可能會導致不同的最終群集模型。對記錄順序的解決方案的魯棒性可以通過多次擬合一個模型來評估，每次使用不同的隨機順序記錄。

評價模型質量的方法采用輪廓系數(shù)(Silhouette Coefficient)[10-14]，是聚類效果好壞的一種評價方式。最早由Leonard kaufman，Peter J.Rousseeuw 在1986年提出。它結合內聚度和分離度兩種因素?？梢杂迷谙嗤紨?shù)據(jù)的基礎上用來評價不同算法，或者算法不同運行方式對聚類結果所產生的影響。

對數(shù)據(jù)中每一個數(shù)據(jù)樣本與本結果簇中其他樣本間的相似性以及該數(shù)據(jù)樣本與其他結果簇中樣本的相似性進行定量分析，并將定量結果進行某種形式的組合，得到聚類結果優(yōu)劣的評價。

對于第i個樣本，計算該樣本到其所屬簇中所有樣本的平均距離，記為a(i)，對于第i個樣本和不包含該樣本的任意簇，計算該樣本到給定簇的中所有樣本的平均距離，所有簇的平均距離最小者記為b(i)。對于第i個樣本，其輪廓系數(shù)的計算如下公式(1)所示。

(1)

式中：a(i) = average(i向量到所有它屬于的簇中其他點的距離)；b(i) = min(i向量到所有非本身所在簇的點的平均距離)。

a(i)越小，說明樣本i越應該被聚類到該簇。b(i)越大，說明樣本i越不屬于其他簇。

輪廓系數(shù)s(i)=+1時，表示 樣本i與其他簇中的對象相異性較大, 聚類合理。

輪廓系數(shù)s(i)=0時，表示樣本i分類不明顯,在兩個簇的邊界上。

輪廓系數(shù)s(i)= -1時，表示樣本i被分配到一個錯誤的簇中，輪廓系數(shù)的值在-1和+1之間變化。

所有樣本的s(i)的均值稱為聚類結果的輪廓系數(shù)，是該聚類是否合理、有效的一種度量方式。

1.2.3.2 評價指標的重要度算法 考察地力指標的重要性，可使用Pearson樣本相關系數(shù)、似然比(Likelihood Ratio)卡方值和Cramer’V系數(shù)。SPSS Modeler計算Cramer’sV系數(shù)[9]，見公式2。

(2)

式中：χ2是Pearson卡方統(tǒng)計量，R和C分別表示列聯(lián)表的行數(shù)和列數(shù)。Cramer’s V系數(shù)取值范圍在0～1，越接近1，越說明輸入變量與輸出變量的相關性越強，輸入變量越重要。重要度在0～ 0.2 范圍為差，重要度在0.2～ 0.6范圍為中， 重要度在0.6～1范圍為良。

2 結果與分析

2.1 地力指標的篩選結果初次運算，輸入地力特征參數(shù)包括成土母質、地貌類型、坡向、pH值、全氮、速效鉀、有機質、年平均降水量、北緯、東經、有效土層厚度、海拔等在內的14個指標，特征重要度排名前10的指標如圖1所示。年平均降雨量特征重要度為1，全氮特征重要度為0.67，有效土層厚度特征重要度為0.58，緯度特征重要度為0.58，成土母質特征重要度為0.51，有機質特征重要度為0.46，經度特征重要度為0.36,地貌類型特征重要度為0.21,坡向特征重要度為0.18，有效磷特征重要度為0.11。

圖1 初次運算的地力指標的特征重要度

根據(jù)特征重要度的大小，選擇特征重要度>0.2以上的地力指標作為模型聚類的特征值，經過反復迭代篩選，以輪廓系數(shù)最優(yōu)為依據(jù)，最終確定地力指標為年平均降雨量、全氮、有效土層厚度、緯度和成土母質(圖2)，年平均降雨量重要度為1，緯度的重要度為0.9，全氮的重要度為0.75，有效土層厚度的重要度為0.59，成土母質的重要度為0.57。從圖2可以看出，影響海南膠園地力的5個指標中，土壤方面有3個，分別是全氮、成土母質、有效土層厚度，這3個因子一起綜合在一定程度上代表了不同土壤的養(yǎng)分總量，因有機質含量與全氮有較大的相關性，選擇特征重要度>0.2重新迭代時就隱含了；地理方面有1個指標，就是緯度，這是一個影響氣溫和熱量的十分重要的因子；氣象方面有一個指標，就是年均降雨量，影響的是土壤水分?？梢姡?、熱、養(yǎng)分，作為熱帶樹種橡膠樹來說，是其生長發(fā)育和高產高效的前提條件，也就是構成地力的重要因子。

圖2 最終確定的地力指標的特征重要度

2.2 TwoStep-AS模型聚類結果根據(jù)最終確定的地力指標(常年降雨量、緯度、全氮、有效土層厚度和成土母質)建立TwoStep-AS聚類模型。模型參數(shù)信息如下：最小常規(guī)聚類數(shù)2類，最大常規(guī)聚類數(shù)15類。自適應特征選擇，信息標準采用貝葉斯信息標準 (BIC)，距離測量采用對數(shù)似然。最終模型常規(guī)聚類數(shù)為3個，即將海南墾區(qū)十大片區(qū)的耕地劃分為3種類型。3種耕地類型的概要特征如表4所示，每種耕地類型最重要的特征在表中以粗體顯示。第1種耕地類型最重要的特征是年平均降雨量在1 600～1 700 mm，全氮質量含量在0.5～0.65 g·kg-1,有效土層厚度在70～100 cm，成土母質為玄武巖。第2種耕地類型最重要的特征是年平均降雨量大于2 000 mm，全氮質量含量在0.05～0.10 g·kg-1,有效土層厚度在55～70 cm，成土母質為花崗巖。第3種耕地類型最重要的特征是年平均降雨量在1 700～1 800 mm，全氮質量含量在0.65～0.85 g·kg-1,有效土層厚度在70～100 cm，成土母質為花崗巖。

表4 3種耕地類型的聚類中心

注： 最重要的特征用粗體表示

Note: The most important characteristics are in bold.

2.3 模型評估在SPSS Modeler中，提出了2種度量聚類效果的指標，一個是優(yōu)度，另一個是重要度。優(yōu)度是一種測量聚類內聚性和分離性的指標。優(yōu)度在-1～ 0.2范圍為差，優(yōu)度在0.2～0.5范圍為中，優(yōu)度在0.5～1范圍為良。重要度是一種測量聚類內聚性的指標。重要度在0～0.2 為差，重要度在0.2～0.6范圍為中， 重要度在0.6～1范圍為良。輪廓系數(shù)是聚類效果好壞的一種評價方式，輪廓系數(shù)在-1～0.2范圍為差，在0.2～0.5范圍為中，在0.5～1范圍為良。用TwoStep-AS聚類算法為基于海南農墾測土調查數(shù)據(jù)進行分類時，總體模型優(yōu)度(平均輪廓系數(shù))為0.38。各聚類的優(yōu)度和重要度如表5所示。從表5的數(shù)據(jù)來看，各聚類類型的優(yōu)度為中等，而重要度為良好，說明用TwoStep-AS聚類算法對海南農墾測土調查數(shù)據(jù)進行聚類時內聚性一般，而分離性很好。

表5 模型質量

3 結 論

1)本研究使用基于大數(shù)據(jù)的TwoStep-AS聚類模型對海南農墾測土調查數(shù)據(jù)進行聚類分析，將海南墾區(qū)十大片區(qū)劃分為3種耕地類型。第1種耕地類型最重要的特征是年平均降雨量在1 600～1 700 mm，全氮質量含量在0.5～0.65 g·kg-1范圍,有效土層厚度在70～100 cm范圍，成土母質為玄武巖。第2種耕地類型最重要的特征是年平均降雨量大于2 000 mm，全氮質量含量g·kg-1在0.05～0.10 g·kg-1范圍,有效土層厚度在55～70 cm，成土母質為花崗巖。第3種耕地類型最重要的特征是年平均降雨量在1 700～1 800 mm，全氮質量含量在0.65～0.85 g·kg-1范圍,有效土層厚度在70～100 cm，成土母質為花崗巖。2)影響海南膠園耕地的最重要的地力指標分別是年平均降雨量、全氮、北緯、有效土層厚度、成土母質。3)基于數(shù)據(jù)挖掘的地類劃分沒有人為主觀因素的干擾，與作物產量無關，以耕地本身特征內在的凝聚和分離為依據(jù)，排除了人為主觀因素，這種分類將更有助于人們對耕地性質的認識，有利于土地的區(qū)劃和綜合利用。4)本研究結果為建立統(tǒng)一的海南墾區(qū)膠園地力評價指標體系提供了初步的方法和結果，并為下一步開展系統(tǒng)性的研究打下了基礎，其結果的準確性和適用性還有待今后在深入研究和應用中得到檢驗。