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

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

(1. 海南大學(xué) 計(jì)算機(jī)與網(wǎng)絡(luò)空間安全學(xué)院， 海口 570228；2. 海南省農(nóng)墾科學(xué)院， 海口 570206；3.中國熱帶農(nóng)業(yè)科學(xué)院 橡膠研究所 海口 571101)

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

1 研究對象與技術(shù)方法

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

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

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

表1 地力評(píng)價(jià)各分類指標(biāo)調(diào)查數(shù)據(jù)調(diào)整前后類別數(shù)

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

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

表3 立地指標(biāo)分組標(biāo)準(zhǔn)表

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

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

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

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

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

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

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

評(píng)價(jià)模型質(zhì)量的方法采用輪廓系數(shù)(Silhouette Coefficient)[10-14]，是聚類效果好壞的一種評(píng)價(jià)方式。最早由Leonard kaufman，Peter J.Rousseeuw 在1986年提出。它結(jié)合內(nèi)聚度和分離度兩種因素。可以用在相同原始數(shù)據(jù)的基礎(chǔ)上用來評(píng)價(jià)不同算法，或者算法不同運(yùn)行方式對聚類結(jié)果所產(chǎn)生的影響。

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

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

(1)

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

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

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

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

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

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

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

(2)

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

2 結(jié)果與分析

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

圖1 初次運(yùn)算的地力指標(biāo)的特征重要度

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

圖2 最終確定的地力指標(biāo)的特征重要度

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

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

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

Note: The most important characteristics are in bold.

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

表5 模型質(zhì)量

3 結(jié) 論

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