農(nóng)作物冠層光譜信息的施肥管理分區(qū)研究

陳 浩，王 熙*，張 偉，王新忠，狄小冬，王 暢

1. 黑龍江八一農(nóng)墾大學工程學院，黑龍江 大慶 163319 2. 黑龍江八一農(nóng)墾大學理學院，黑龍江 大慶 163319

引 言

黑龍江省是中國的重要商品糧食種植基地，種植了大量的玉米和大豆，從國家官方媒體公布的數(shù)據(jù)顯示，2020年黑龍江省糧食總產(chǎn)量已經(jīng)達到了7 541萬噸，玉米和大豆的種植產(chǎn)量連續(xù)數(shù)年位居全國前列。 隨著玉米、 大豆等農(nóng)作物的產(chǎn)量逐年攀升，出現(xiàn)了不同程度的化肥濫用情況，化肥的濫用不僅造成了一定的經(jīng)濟損失，也對種植土壤造成了不可逆的損害，施肥管理分區(qū)概念的引入，正是應對化肥利用率下降的有效方法。 基于管理分區(qū)可以進行有針對性的變量施肥，通過劃分管理分區(qū)，能夠快速辨別各個分區(qū)的作物長勢、 土壤養(yǎng)分情況等差別，根據(jù)分區(qū)間的差異精準管理，提高肥料利用率，減少因濫用化肥導致的生態(tài)環(huán)境破壞[1]。 施肥管理分區(qū)的劃分多以土壤養(yǎng)分的含量、 作物的產(chǎn)量及農(nóng)作物冠層光譜特征作為信息源，利用k-均值(k-means)算法、 模糊c-均值(fuzzy c-means，F(xiàn)CM)算法、 層次聚類(hierarchical clustering)算法以及衍生出的相關優(yōu)化算法進行分區(qū)計算。 在各種劃分算法中，F(xiàn)CM算法的應用最為廣泛。 朱昌達[2]基于農(nóng)作物土壤的特性，利用FCM算法計算并劃分管理分區(qū)，劃分的管理分區(qū)與農(nóng)作物的土壤特性具有較高的相關性，可利用劃分的管理分區(qū)進行差異管理。 陳世超[3]基于種植農(nóng)作物的地形、 土壤、 電導率等因素，使用FCM算法劃分了管理分區(qū)，驗證了可以根據(jù)各區(qū)域的不同特點劃分管理分區(qū)。

隨著光譜技術的發(fā)展，比值植被指數(shù)、 增強型植被指數(shù)、 歸一化植被指數(shù)(normalized difference vegetation index，NDVI)等多種光譜植被指數(shù)被提出，依據(jù)植被指數(shù)來反映農(nóng)作物的生長特性的研究也越來越多，其中NDVI的研究最多[4]。 利用光譜傳感器測定植物冠層對紅光和近紅外光的反射率可以得到NDVI數(shù)據(jù)，從而反映植株的氮素營養(yǎng)情況并表征農(nóng)作物長勢[5]。 植物冠層光譜傳感器有無光源傳感器和有光源傳感器，其中無光源傳感器以日光作為光源生成光譜圖像估計作物生長參數(shù)，傳感器以光譜參數(shù)及像素準確描述植物生長狀態(tài)。 衛(wèi)星和飛機基于圖像光譜傳感器以這種方式工作，但衛(wèi)星易受到云層及大氣因素干擾，無人機雖無此問題，但續(xù)航時間較短。 有源傳感器自帶光源照射作物表面，對反射回來的信號進行測量并進行相應計算，受外界因素干擾小，工作穩(wěn)定，時效性強。 Greenseeker作為農(nóng)作物冠層光譜信息采集傳感器，在各個國家的應用都較為廣泛，比如美國約翰迪爾公司開發(fā)的的綠色之星變量噴施系統(tǒng)就使用Greenseeker傳感器實時檢測反映農(nóng)作物長勢信息的光譜數(shù)據(jù)，并利用光譜數(shù)據(jù)劃分施肥管理分區(qū)進行精準施肥，Greenseeker廣泛應用于指導變量施肥，具有良好的研究基礎和應用前景。

本工作使用有光源的地面遙感植物光譜檢測儀Greenseeker獲取大豆及玉米的NDVI數(shù)據(jù)，再利用基于FCM算法提出的基于模型的模糊c-均值算法(model-based fuzzy c-means，MFCM)對NDVI數(shù)據(jù)進行劃分計算，實現(xiàn)施肥管理分區(qū)的劃分，利用管理分區(qū)劃分效果評價指標輪廓系數(shù)和調(diào)整蘭德系數(shù)對MFCM算法結果進行了評價。

1 實驗部分

1.1 地點

實驗地點位于中國黑龍江省趙光農(nóng)場的第4管理區(qū)中的17作業(yè)站11號和12-2號地塊，如圖1所示。 趙光農(nóng)場位于黑龍江省大興安嶺東側，小興安嶺北側，整個地區(qū)群山起伏，海拔240～330 m。 該地區(qū)的氣候類型是寒溫帶大陸性季風氣候，種植農(nóng)作物的主要土壤類型是黑土，其天然條件有利于農(nóng)作物冠層光譜信息采集實驗的開展。 玉米的NDVI數(shù)

圖1 趙光農(nóng)場第4管理區(qū)17作業(yè)站11號地與12-2號地Fig.1 No.11 and No.12-2 of No.17 operation stationin No.4 management area of Zhaoguang farm

據(jù)采集地塊為11號地，地塊面積為618畝，采集的時間為2019年6月18日； 大豆的NDVI數(shù)據(jù)采集地塊為12-2號地，地塊的面積為700畝，采集的時間為2019年6月22日。

1.2 數(shù)據(jù)采集

通過搭建的農(nóng)作物冠層光譜信息采集平臺進行大豆及玉米的NDVI數(shù)據(jù)采集，采集平臺由6個Greenseeker植物冠層光譜檢測儀、 配備10英寸顯示屏一體機、 天寶公司生產(chǎn)的AG332型衛(wèi)星定位系統(tǒng)、 NDVI數(shù)據(jù)記錄器以及國產(chǎn)運動控制器組成，如圖2所示。 Greenseeker植物冠層光譜檢測儀由美國天寶公司生產(chǎn)，自帶兩個窄帶LED主動光源，可發(fā)射650 nm(誤差為10 nm)的紅光和770 nm(誤差為10nm)近紅外光，利用NTech公司的第二代光學傳感器作為檢測傳感器獲取兩波段下植物冠層光譜反射率信息[6]。 Greenseeker植物冠層光譜檢測儀檢測過程中不需要外界輔助光源，不受云層遮擋和土壤反射光的影響，可以全天候進行植物冠層光譜信息的采集作業(yè)，還具有檢測速度快、 數(shù)據(jù)較為準確的特點，可以實時反映植物的生長狀況，從而有效地指導依據(jù)作物長勢的變量施肥[7-8]。 如圖3所示，實驗中將農(nóng)作物冠層光譜信息采集平臺安裝到凱斯2254型拖拉機上，在中耕期進行大豆及玉米的NDVI數(shù)據(jù)采集，利用NDVI數(shù)據(jù)進行施肥管理分區(qū)的劃分，根據(jù)分區(qū)結果進行有針對性的施肥。

圖2 農(nóng)作物冠層光譜信息采集平臺Fig.2 Spectral information collection platform for crop canopy

圖3 農(nóng)作物冠層光譜信息采集現(xiàn)場Fig.3 Field of spectral information collection of crop canopy

1.3 劃分施肥管理分區(qū)方法

1.3.1 FCM算法

1969年Ruspini在模糊集合理論的基礎上提出了一種模糊劃分的概念。 1974年Dunn又依據(jù)Ruspini所提出的模糊劃分概念，把硬c-均值算法推廣到了模糊劃分方向，形成了最初步的FCM算法，后續(xù)又有眾多學者對FCM算法進行了不斷的優(yōu)化和改進。 FCM算法屬于柔性的數(shù)據(jù)劃分算法，在所有模糊劃分算法中的應用位居前列。 FCM算法在優(yōu)化目標函數(shù)的基礎上可以獲得每個數(shù)據(jù)在所有類別中的隸屬度，從而可以決定數(shù)據(jù)的歸屬。 對于數(shù)據(jù)集合X=[x1, …,xn]，劃分中心個數(shù)為c，模糊權重指數(shù)為m，F(xiàn)CM算法在硬c-均值的基礎上，引入了模糊隸屬度矩陣U，矩陣中包含隸屬度uij(0≤uij≤1；i=1，2，…，c；j=1，2，…，n)，uij的值代表某一個數(shù)據(jù)到某一個劃分的程度，劃分中心為V=[v1，v2, …,vc]，該算法的目標函數(shù)為

(1)

式(1)中，‖xj-vi‖2為第j個數(shù)據(jù)到第i個劃分中心vi的歐式距離，通常1.5≤m≤2.5。 隸屬度矩陣U滿足

(2)

FCM算法的目標函數(shù)達到最小值的條件是

(3)

和

(4)

FCM算法就是基于上述目標函數(shù)(1)，不斷尋找該函數(shù)最小值的方法。 FCM算法的具體步驟為： 首先設置劃分中心個數(shù)c、 停止閾值ε、 模糊權重系數(shù)m，隨機選取初始劃分中心； 然后根據(jù)歐氏距離，更新劃分隸屬度矩陣U； 再更新劃分中心V； 最后判斷劃分中心變化是否小于停止閾值ε，如小于停止閾值ε，則停止劃分程序，輸出U和V，否則再次更新劃分隸屬度矩陣U。 該算法迭代終止后，模糊隸屬度矩陣U對應為樣本數(shù)據(jù)的模糊劃分。

1.3.2 肘部法則

在數(shù)據(jù)劃分分析中一般使用肘部法則來估計劃分數(shù)c，具體方法是利用設置好的劃分數(shù)c來計算代價函數(shù)的函數(shù)值，對比由每個劃分數(shù)c的值計算出的代價函數(shù)值曲線的畸變程度來確定劃分數(shù)c的值。 其中，代價函數(shù)選取的是誤差平方和(sum of squares of errors，SSE)，SSE是指每個樣本點到樣本點所在劃分中心的距離之和，計算公式如式(5)

(5)

式(5)中，c為劃分數(shù)，Ci為第i分區(qū)數(shù)據(jù)集，xj為第i分區(qū)某個數(shù)據(jù),vi為Ci的分區(qū)中心。

1.3.2 MFCM算法

MFCM算法是在FCM算法的基礎上引入了分區(qū)計算模型，在一定數(shù)據(jù)量范圍內(nèi)不使用FCM算法進行迭代劃分數(shù)據(jù)，而是使用分區(qū)計算模型進行數(shù)據(jù)劃分，該方法在滿足基本精度的前提下，可有效降低計算耗時，隨著數(shù)據(jù)量的增加，效果更加明顯。 具體算法為： 首先按照用戶提供的分區(qū)數(shù)隨機設置劃分中心，然后利用FCM算法進行2 000個數(shù)據(jù)量的劃分計算并建立分區(qū)計算模型，最后以2 000個數(shù)據(jù)量為分界點，利用分區(qū)計算模型劃分接下來的數(shù)據(jù)，算法流程如圖4所示。 使用Python函數(shù)庫中的Scikit-Fuzzy模糊邏輯算法編寫FCM算法和MFCM算法，分別對玉米和大豆的NDVI數(shù)據(jù)進行分區(qū)計算。 具體使用到Scikit-fuzzy工具箱中的cmeans工具，將Greenseeker傳感器獲取的NDVI數(shù)據(jù)作為程序的輸入數(shù)據(jù)，分區(qū)數(shù)量設置為4個，隸屬度指數(shù)設置為2，當隸屬度的變化小于0.005時提前結束迭代，最大迭代次數(shù)設置為1 000。

圖4 MFCM算法流程圖Fig.4 Flow chart of model-based fuzzy c-means algorithm

1.4 施肥管理分區(qū)劃分結果評價方法

1.4.1 輪廓系數(shù)

1986年，Rousseeuw提出輪廓系數(shù)(silhouette coefficient，SC)的概念，利用輪廓系數(shù)可以評價數(shù)據(jù)劃分結果的好壞，也可以評價各種算法對同一組數(shù)據(jù)劃分的優(yōu)劣效果，輪廓系數(shù)主要是利用數(shù)據(jù)間的結合聚合程度和分離程度兩種因素來進行評價計算[9-10]。 輪廓系數(shù)的計算公式為

(6)

具體有

(7)

輪廓系數(shù)SC(i)的值在[-1, 1]區(qū)間內(nèi)變化，當輪廓系數(shù)的值接近1的時候，表明數(shù)據(jù)間的結合聚合程度和分離程度都比較好，此時的劃分效果最好。

1.4.2 調(diào)整蘭德系數(shù)

可以把數(shù)據(jù)劃分過程看成一系列決策的過程，列舉數(shù)據(jù)集中的每一對數(shù)據(jù)，判斷每一對數(shù)據(jù)是否應該分到同一個簇中，蘭德指數(shù)(Rand index，RI)衡量了其中正確決策的比例[11]。 假設U為理想劃分結果集合，V為待評價的劃分結果集合。 假設存在4個統(tǒng)計量：a為在集合U和集合V中都是同一劃分的數(shù)據(jù)對數(shù)；b為在集合U中是同一劃分，但是在集合V中卻不是同一劃分的數(shù)據(jù)對數(shù)；c為在集合U中是同一劃分，但在集合V中卻不是同一劃分的數(shù)據(jù)對數(shù)；d為在集合U和集合V中都不是同一劃分的數(shù)據(jù)對數(shù)。 則RI的表達式如式(8)

(8)

式(8)中，0≤RI≤1，當待評價的劃分結果與理想劃分結果完美匹配時，RI的值為1，當待評價的劃分結果與理想劃分結果完全不匹配時，RI的值為0。 但是RI不能保證隨機劃分結果的值接近0，調(diào)整蘭德指數(shù)(adjusted Rand index，ARI)可以解決這個問題[12]。 ARI的值越大說明待評價的劃分結果與理想劃分結果越相似，ARI的定義如式(9)

(9)

式(9)中，E[RI]為RI的期望值，max(RI)為RI的最大值。

2 結果與討論

2.1 肘部法則計算施肥管理分區(qū)數(shù)

利用Python編寫基于FCM算法的肘部法則計算程序，通過程序計算大豆及玉米的NDVI數(shù)據(jù)在不同分區(qū)數(shù)下的SSE值，并繪制成SSE變化曲線圖，如圖5所示。 對于大豆的NDVI數(shù)據(jù)，當分區(qū)數(shù)大于2時，SSE的值不再有明顯的下降趨勢，后續(xù)曲線也比較平滑，說明大豆的施肥管理分區(qū)數(shù)在2以上均是合理的。 對于玉米的NDVI數(shù)據(jù)，當分區(qū)數(shù)大于3時，SSE值沒有明顯的下降趨勢，曲線也變得較為平滑，說明玉米的施肥管理分區(qū)數(shù)為3以上是合理的。 施肥管理分區(qū)的數(shù)量將直接影響施肥作業(yè)的難度，分區(qū)數(shù)的增加也相應伴隨著作業(yè)成本的增加； 本工作綜合考慮農(nóng)作物施肥場地的條件和農(nóng)場以往的施肥經(jīng)驗確定施肥管理分區(qū)數(shù)為4。

2.2 施肥管理分區(qū)劃分結果評價

SC是常用的施肥管理分區(qū)劃分效果內(nèi)部評價指標，反映了數(shù)據(jù)點與所在劃分中心的密集程度，以及與其他劃分中心疏遠的程度，一個好的劃分計算，分區(qū)內(nèi)部的數(shù)據(jù)要足夠的聚集，不同分區(qū)的數(shù)據(jù)應該足夠疏遠，文獻[13-14]利用SC評價了不同聚類的質(zhì)量，取得了良好的效果。 通過計算不

圖5 NDVI數(shù)據(jù)在不同分區(qū)數(shù)下的SSE變化曲線(a)： 大豆； (b)： 玉米Fig.5 SSE change curves of NDVI data under differentpartition numbers(a)： Soybean； (b)： Maize

圖6 不同NDVI數(shù)據(jù)量的輪廓系數(shù)(a)： 大豆； (b)： 玉米Fig.6 Silhouette coefficients of different NDVI data volumes(a)： Soybean； (b)： Maize

同NDVI數(shù)據(jù)量下的大豆和玉米施肥管理分區(qū)劃分評價指標SC，得到SC在不同NDVI數(shù)據(jù)量下的變化曲線，如圖6所示。 從圖中可以看出，對于大豆NDVI數(shù)據(jù)，兩種算法相差最多為0.005； 對于玉米NDVI數(shù)據(jù)，兩種算法相差最多0.022。 總體來看，無論對于大豆還是玉米，兩種算法劃分施肥管理分區(qū)的效果均相差不大。

MFCM算法在計算施肥管理分區(qū)過程中利用少部分數(shù)據(jù)建立分區(qū)模型，肯定會降低計算精度； 為了檢驗在不同數(shù)據(jù)量下對計算精度的影響，計算了不同NDVI數(shù)據(jù)量下的大豆和玉米施肥管理分區(qū)劃分評價指標ARI。 ARI是對RI進行了改進聚類效果外部評價指標，試圖表達怎樣聚類才是正確的，通過計算兩個不同聚類的相似性來評價聚類的效果好壞。 文獻[15]在評估單細胞RNA-sep數(shù)據(jù)聚類方法中，利用ARI評價方法評估了5種聚類方法的穩(wěn)定性，取得了良好的效果。 對比數(shù)據(jù)為FCM算法的施肥管理分區(qū)劃分，計算結果如圖7所示。 通過分析大豆的ARI變化曲線，可以發(fā)現(xiàn)當數(shù)據(jù)量為超過6000后，ARI的值維持在0.8左右，當數(shù)據(jù)量為超過8 000后，ARI的值維持在0.9左右，說明數(shù)據(jù)量超過6 000后，兩種算法的劃分精度幾乎接近。 通過分析玉米的ARI變化曲線，可以發(fā)現(xiàn)當數(shù)據(jù)量達到8 000后，ARI的值維持在0.8左右，但是在數(shù)據(jù)量達到10 000時，ARI的值突降到0.5以下，說明在8 000～10 000的數(shù)據(jù)的值波動較大，造成MFCM算法此前建立的分區(qū)模型不再適應后續(xù)2 000個NDVI值，但是此時輪廓系數(shù)相差僅為0.017，說明兩種算法的劃分結果相差較大，但本身的劃分效果相差不大，具體原因需后續(xù)進行深入研究。

通過計算不同NDVI數(shù)據(jù)量下的大豆和玉米施肥管理分區(qū)劃分耗時，得到兩種算法在不同NDVI數(shù)據(jù)量下的耗時變化曲線圖，如圖8所示。 從圖中可以看出，無論對于大豆還是玉米，MFCM耗時均比FCM耗時少，耗時差值最大時達到0.51 s。

圖7 不同NDVI數(shù)據(jù)量的調(diào)整蘭德指數(shù)(a)： 大豆； (b)： 玉米Fig.7 Adjusted Rand index of different NDVI data volumes(a)： Soybean； (b)： Maize

圖8 兩種算法在不同NDVI數(shù)據(jù)量下的計算耗時(a)： 大豆； (b)： 玉米Fig.8 Time costs of two algorithms under different NDVI data volumes(a)： Soybean； (b)： Maize

2.3 不同算法的大豆施肥管理分區(qū)圖對比

ArcGIS軟件是一款非常優(yōu)秀的地理信息處理及繪圖軟件，可以將農(nóng)作物的NDVI數(shù)據(jù)基于地理信息(經(jīng)緯度數(shù)據(jù))進行插值計算，繪制出基于農(nóng)作物NDVI的施肥管理分區(qū)圖。 當大豆NDVI數(shù)據(jù)量為4 000時，兩種算法的施肥管理分區(qū)如圖9所示，通過分區(qū)圖對比兩種算法的劃分效果。 FCM算法和MFCM算法在施肥管理分區(qū)圖上的差異較為明顯，與施肥管理分區(qū)劃分評價指標ARI值的變化相符合，主要差別體現(xiàn)在圖中的綠色和藍色分區(qū)部分。 出現(xiàn)明顯差異的原因是此數(shù)據(jù)量下的NDVI數(shù)據(jù)樣本較少，說明MFCM算法在作業(yè)初期將會出現(xiàn)施肥管理分區(qū)劃分相對不準確情況，在實際應用中建議早期使用推薦施肥量進行施肥，不需根據(jù)長勢調(diào)整施肥量。

圖10為大豆NDVI數(shù)據(jù)量8 000時兩種算法的施肥管理分區(qū)圖，從圖中可以看出，兩種算法的施肥管理分區(qū)圖相差不大，與施肥管理分區(qū)劃分評價指標ARI的變化相符合，當NDVI數(shù)據(jù)量達到8 000時，ARI的值為0.86，說明當大豆NDVI數(shù)據(jù)量達到8 000時，兩種方法的分區(qū)效果相近。

圖9 大豆NDVI數(shù)據(jù)量為4 000時兩種算法計算的施肥管理分區(qū)圖(a)： FCM算法； (b)： MFCM算法Fig.9 Fertilization management zoning calculated by two algorithms for 4 000 soybean NDVI(a)： FCM algorithm； (b)： MFCM algorithm

圖10 大豆NDVI數(shù)據(jù)量為8 000時兩種算法計算的施肥管理分區(qū)圖(a)： FCM算法； (b)： MFCM算法Fig.10 Fertilization management zoning calculated by two algorithms for 8 000 soybean NDVI(a)： FCM algorithm； (b)： MFCM algorithm

2.4 不同算法的玉米施肥管理分區(qū)對比

當玉米NDVI數(shù)據(jù)量為4 000時，兩種算法計算的施肥管理分區(qū)如圖11所示。 對比兩種算法的施肥管理分區(qū)圖，可以發(fā)現(xiàn)，相較于大豆NDVI數(shù)據(jù)，在4 000數(shù)據(jù)量時兩種算法針對玉米的NDVI數(shù)據(jù)的施肥管理分區(qū)圖相差不大，ARI為0.69。 當玉米的NDVI數(shù)據(jù)量為8 000時，兩種方法的施肥管理分區(qū)圖相差不大，如圖12所示。 此時的ARI為0.915，說明此時兩種方法的施肥管理分區(qū)劃分結果高度相似，MFCM能夠在節(jié)省一定時間的基礎上，有一個較好的劃分結果。

3 結 論

使用Python編寫了FCM算法和MFCM算法，對大豆和玉米兩種作物進行了施肥管理分區(qū)劃分，并利用SC和ARI進行劃分效果評價。 從計算結果來看，當NDVI數(shù)據(jù)量低于4 000時，MFCM算法與FCM算法相比，沒有明顯優(yōu)勢，但此時涉及的作業(yè)面積較小，對整體施肥管理影響不大，建議在此數(shù)據(jù)量下使用推薦施肥量施肥。 當NDVI數(shù)據(jù)量大于6 000時，MFCM算法在劃分耗時上出現(xiàn)明顯優(yōu)勢，

圖12 玉米NDVI數(shù)據(jù)量為8 000時兩種算法計算的施肥管理分區(qū)圖(a)： FCM算法； (b)： MFCM算法Fig.12 Fertilization management zoning calculated by two algorithms for 8 000 soybean NDVI(a)： FCM algorithm； (b)： MFCM algorithm

能夠在節(jié)省計算耗時的基礎上保持一定的劃分精度； 從SC和ARI的結果也能夠體現(xiàn)出來，兩種方法的SC相差最多為0.02，ARI能夠維持在0.8左右。 通過計算不同NDVI數(shù)據(jù)量下的大豆和玉米施肥管理分區(qū)劃分耗時，利用本文提出的MFCM算法可以有效降低劃分耗時，隨著數(shù)據(jù)量的增加，這種優(yōu)勢體現(xiàn)的更加明顯。 但是MFCM劃分算法也有一定的局限性，玉米NDVI數(shù)據(jù)量在8 000～10 000范圍內(nèi)的數(shù)據(jù)劃分計算發(fā)現(xiàn)，此時ARI的值僅為0.475，說明兩種算法劃分的施肥管理分區(qū)結果出現(xiàn)了明顯差異，原因是此段數(shù)據(jù)變化對劃分中心的影響較大，MFCM算法應對數(shù)據(jù)量突變的情況沒有很好的適應性，需要繼續(xù)優(yōu)化程序，后續(xù)將針對此問題進行深入研究。