黑土區田塊尺度精準管理遙感分區時空格局與成因分析

劉煥軍，殷 悅，鮑依臨，張新樂，馬雨陽，

王夢沛1，孟令華2，宋少忠3※

（1. 東北農業大學公共管理與法學院，哈爾濱 150030；2. 中國科學院東北地理與農業生態研究所，長春 130012；3. 吉林工程技術師范學院信息工程學院，長春 130052）

0 引 言

通過農田分區對土壤和農作物實施變量投入管理是近年來精準農業的研究熱點[1]。分區管理需要將區域劃分為具有同質性的多個管理區[2]，劃分精準管理區有利于防止過量施肥、施藥造成的環境污染[3]，也是精準農業實踐、提高糧食產量的基礎[4-5]。

國內外學者對分區輸入量及分割方法做了大量研究。通過已有研究來看，傳統精準管理分區輸入量以田間采樣所測得的土壤理化性質數據[6-7]、作物產量數據[8-9]為主，分區方法主要為K-means、C-means等聚類算法[10-11]。然而土壤網格取樣需要花費大量人力物力，產量數據一年僅可獲取一次，降低了精準管理效率，近年來精準管理分區輸入量也逐漸過渡為通過遙感數據獲取的土壤及植被信息。

歸一化差異植被指數（Normalized Difference Vegetation Index, NDVI）受氣象、環境等多種因素影響[12]，是植被生長狀況的綜合體現，與生物量、葉面積指數等關鍵作物生長指標密切相關[13]，并能很好地反映作物產量[14]，Gallego等[15]探究了NDVI與產量的關系，利用普通及回歸克里格方法繪制了番茄產量預測圖，因此，NDVI也常被用于劃分農田管理區。Martins等[16]將NDVI、土壤電導率等數據進行主成分分析后運用模糊K-means聚類方法對大豆田劃分管理區，認為NDVI可以反映產量及土壤屬性變異性。劉煥軍等[17]運用面向對象的多尺度分割方法綜合同年作物生長期四期NDVI劃分管理區，并運用莫蘭指數進行評價，結果表明基于多期NDVI的分區精度較土壤有機質插值的精度高。Santos等[18]提出了一種基于無人駕駛飛行器捕獲的作物NDVI運用K-means聚類算法劃分管理區的系統，可以實現對農田管理區實時、快速、準確的劃分。

除關注輸入量與分割方法以提高分區精度外，精準管理分區研究還應關注同一區域精準管理分區結果的規律及其成因，這是本領域關注較少的角度，也是未來的研究趨勢。本研究選取東北典型黑土區田塊，獲取研究區2017－2020年4期6月上旬遙感影像并提取NDVI，運用面向對象分割的方法劃分精準管理區，基于分區結果借助空間轉移矩陣研究精準管理分區格局，并嘗試從微地形的視角研究精準管理分區時空格局成因。

1 材料與方法

1.1 研究區概況

研究區位于黑龍江省東部佳木斯市的友誼農場，地處三江平原腹地，處于松花江第二階地，海拔在59～300 m。地塊中心經緯度為46°74′ N，131°67′ E，面積約為23.42 hm2，屬于黑土區，漫川漫崗。四季分明，屬中溫帶大陸性季風氣候，雨量充沛，降水多集中于夏季。常年平均氣溫3.1 ℃，極端最高氣溫38 ℃，極端最低氣溫-37 ℃，≥10 ℃有效活動積溫2 300～2 600 ℃，平均日照時數2 730 h，年降水量為500 mm左右，無霜期143 d左右。研究區田塊地勢東北及東南部高，西部及中部低。田塊內部作物長勢差異顯著，且多年格局相似，適宜進行分區管理，因而作為本次研究對象，所選研究區具有典型性。研究區位置如圖1所示，2017－2020年以種植玉米供試，且田間管理措施相同。

1.2 數據源

1.2.1 遙感影像獲取及處理

研究區田塊6月上旬作物長勢差異顯著，且多年具有規律性，故通過歐空局網站（https://sentinel.esa.int/web/sentinel）獲取研究區2017－2020年6月上旬共4幅Sentinel-2A遙感影像，采集日期分別為2017年6月9日、2018年6月9日、2019年6月2日、2020年6月13日，下載產品級別為Level-1 C，為經過正射校正和幾何精校正的大氣表觀反射率產品。通過歐空局(ESA) 提供的SNAP軟件對影像分別進行輻射定標、大氣校正及重采樣，處理后為Level-2 A級別產品，空間位置可與通過RTK點位校正的無人機DEM、實地采樣RTK定位點相對應，重采樣后的空間分辨率為10 m。在ArcGIS 10.2中根據研究區范圍對影像進行裁剪并計算NDVI。

1.2.2 地形數據的獲取與處理

2019年4月10日由大疆-精靈Phantom 4 RTK無人機拍攝獲得研究區DEM。選取視野開闊，地勢平坦地區作為無人機起降場，于拍攝區均勻布設相紙作為控制點，并通過定位精度為厘米級的RTK進行定位以進行后期影像幾何校正。參照《低空數字航空攝影測量外業規范》，設置無人機航拍的航向重疊度為90%，旁向重疊度為75%，飛行高度110 m。航拍路線設置“N”字航線，等時間間隔拍攝，時間間隔為2 s，飛行速度為7.9 m/s。去除拍攝姿態角不佳的圖像，將影像導入專業三維建模軟件PhotoScan進行拼接；根據POS點位置并結合臨近影像的同名點對影像進行匹配，生成三維點云數據；之后對單景影像進行裁剪拼接，并根據控制點位進行校正，誤差可控制在一個像元以內，最后生成研究區的無人機正射DEM，并重采樣為與衛星影像相同分辨率。同時在ArcGIS 10.2中運用3D Analyst模塊的slope工具提取坡度。DEM及坡度圖如圖2。

1.2.3 土壤溫濕度及氣象數據獲取

土壤溫濕度數據由HyDra probe II地埋式無線土壤溫濕度傳感器測得。地埋式傳感器共4個，于2019年到實地埋設，數據采集深度20 cm，通過RTK定位并記錄經緯坐標，埋設位置如圖2 a。土壤水分測試精度±1%，土壤溫度測試精度±0.3℃，無線網絡傳輸，每小時回傳數據，全天候工作。氣象數據由國家氣象科學數據中心（http://data.cma.cn/）查詢獲得。

1.2.4 精準管理分區方法

精準管理分區輸入量選擇研究區多年作物生長初期NDVI，分區方法采用面向對象多尺度分割算法，該算法在Definiens（http://www.definiens.com）平臺開發的eCognition 8.7中實現。通過調節“Scale”（NDVI同質性標準的最大對象分割標準偏差）、“shape”（形狀異質性）、“compact”（緊致度）3個參數來調整分區結果，并提取同質性指數Mean（GLi）和異質性指數STD（ΔGL），分別計算各參數下的分割評價指數(Segmentation Evaluation Index, SEI)。計算公式如下：

式中σ為同質性指標，n為分區內所有像素的個數，GL為灰度均值，GLi為像素i的灰均值，ΔGL為與鄰域的平均差分的絕對值，N為與當前分區鄰接對象的個數，L為周長，Li為與第i相鄰分區公共邊長，GLi為第i個相鄰分區的灰度平均值，SEI為分割評價指數，m為分區輸入量層數。

為了對不同分割尺度下的分區進行比較，引入面積變量計算分區SEI的均值平均分割評價指數（Average Segmentation Evaluation Index, ASEI），在一定尺度下存在一組參數集合（scale, shape, compact）為最優分割參數，其對應的ASEI為最大值[19]。

式中A表示整個研究地塊內所有分區的總面積（m2），Ai表示第i個分區的面積（m2），n表示分區總數量，SEIi表示第i個分區的分割評價指數。

通過計算精準管理分區前后NDVI變異系數（Coefficient of Variation，CV）來證明分區有效性，變異系數CV＜0.1時，為弱變異程度；0.1＜CV＜1時，為中等變異程度；CV＞1時，為強變異程度[20]。變異系數為NDVI標準差σ與平均值μ之比，公式如下：

1.2.5 精準管理分區單元定級

分別確定研究地塊4期精準管理分區邊界后，根據不同單元內作物苗期長勢確定各精準管理分區單元的含義。計算各分區內NDVI平均值，運用Jenks自然斷裂點法將各期NDVI由低到高分為1、2、3三個等級。Jenks自然斷裂點法基于數據內部的內在聯系自然分組，目的是最大化組間差距并最優化組內相似值，分組點選在數據變量值出現相對最大變化處，通過不斷迭代計算分級，使得數據級內變異最小，級間變異最大[21]。計算公式如下[22]：

式中A為1個數組，包含n個數值，1≤i≤j≤n；Ei-j是Ai～Aj一列數的均值。

1.2.6 精準管理分區時空格局分析方法

借助空間轉移矩陣方法對精準管理分區時空格局進行分析。通過計算精準管理分區分級轉移矩陣可將4期分區結果進行關聯分析與可視化，對分區格局進行定量表述，進而探究不同年份同一作物生長期精準管理分區結果之間的規律性。空間轉移矩陣來源于系統分析中對系統狀態與狀態轉移的定量描述[23]，可以定量識別某一要素不同等級在某一時間間隔的空間格局變化，除了反映其不同等級面積變化，還能直觀反映各等級面積轉入轉出情況[24]。它能夠從時間角度揭示出精準管理分區結果的年際轉移變化方向及特征，反映某級別的分區單元在一定時間間隔下，從t時刻向t+1時刻狀態轉移的過程。在ArcGIS10.2中對4期NDVI進行疊加分析，得到2017－2020年逐年間精準管理分區不同等級單元的轉化情況。

2 結果與分析

2.1 基于NDVI的精準管理分區結果

對不同參數下的平均分割評價指數進行計算，確定各期NDVI的最優分割參數。分割尺度、分區數、評價指數及分割前后NDVI變異系數如表1所示，分割結果如圖3。2017－2020年分區數分別為27、25、33、28，從耕作單元角度出發，分區的結果符合當代耕作單元的發展進程，利于實施耕作[25]。由表1可知，對研究區各期NDVI進行分割后，區內NDVI的變異性均達到了弱變異程度，分區后變異系數較分區前降低了70%以上。由圖3可以看出，分區內NDVI相對均一，同一時期作物長勢及精準管理分區結果相近，并且NDVI的空間分異特征與地形分異特征近似。

表1 四期NDVI最優分割尺度、分區數、平均分割評價指數及分割前后變異系數Table 1 NDVI optimal segmentation scale, number of partitions, average segmentation evaluation index and coefficient of variation before and after segmentation

2.2 精準管理分區單元分級與定義

圖4為2017－2020年基于精準管理分區的NDVI分級結果。由于不同年氣候條件及影像采集日期略有差異，故4期結果各級NDVI范圍不同，各期1～3級NDVI范圍依次為：0.152～0.202、＞0.202～0.248、＞0.248～0.338；0.401～0.455、＞0.455～0.502、＞0.502～0.571；＞0.091～0.144、＞0.144～0.223、＞0.223～0.298；＞0.179～0.236、＞0.236～0.289、＞0.289～0.340。根據遙感影像目視分析結合實地考察，對各級精準管理分區單元含義進行界定。NDVI分級為1的分區單元，此分區等級區域內作物長勢較差亦或出苗較小，NDVI值較低，地勢也相對較低且坡度較平緩；NDVI為2級分區內作物長勢中等，也多處于田塊地勢中等高度處；NDVI為3級的分區內部作物總體長勢最佳或出苗較早。

2.3 精準管理分區時空格局分析

通過查詢氣象數據可知，研究區2017－2020年分割前均無極特殊氣候條件（表2），且田塊4a均采取相同的種植及管理模式，因此對4a分割結果的關聯分析可排除年際間氣候及管理要素的影響。綜合圖2～圖4分析可明顯看出，在4a作物生長初期，等級為1級的區域均位于研究田塊西北至中部地區，此部分區域地勢較低、起伏平緩；分區等級為2級的區域主要位于田塊東北部及東南小部分區域，區內地勢稍有起伏；3級的分區集中于南部、西南部及東部，地勢較高、起伏較明顯。4a總體分區分級空間分布格局具有相似性，且與地形特征空間分布格局相吻合。

表2 分割前5 d主要氣候要素Table 2 Main climatic factors in 5 days before segmentation

對4期NDVI分級結果進行疊加分析，得出轉移矩陣如表3所示，矩陣表示縱向向橫向的轉化。從轉移矩陣來看，4年中分區等級為1級向次年同級的轉化面積遠多于向其他兩級的轉化面積，說明在相同時期即作物生長初期長勢較差的位置4a基本穩定在相同的區域。2017－2018年、2018－2019年2級向同級的轉化面積也均大于向其他兩級的轉化，即長勢處于中等水平的區域在年際間也大體保持穩定。2017－2018年、2019－2020年3級只有北部小部分轉向1、2等級，絕大部分仍保持在3級，說明長勢較好的區域年際間也基本穩定。2018－2019年3級轉向2級的面積最大，且2019－2020年2級大部分轉向3級，由表2可知2019年分割前5d降水量最大，且平均氣溫最低，田塊存在積水情況，作物長勢差或出苗較小，導致該年作物整體長勢情況不佳，因而精準管理分區分級的轉移情況略有不同。總體來看，研究區精準管理分區情況在年際間相同時期基本相似，即作物長勢好及長勢差的位置基本處于穩定狀態。

2.4 精準管理分區時空格局成因分析

通過作物長勢及精準管理分區結果可以看出，精準管理分區格局與地形特征分布格局具有一定的相似性。因此計算分區前后高程及坡度的變異系數CV，對比分區對地形因子變異性的影響。表4展示了分區前后高程和坡度的變異程度變化，可以看出分區后高程的變異系數降低了50%左右，坡度的變異系數降低了30%左右，由此可見分區內部不僅作物長勢趨于均一，地形特征也趨于均一，同時也說明研究區作物長勢及精準管理分區結果受地形特征影響較大。

表3 2017－2020年精準管理分區分級轉移矩陣Table 3 Classification transfer matrix of precise management zone from 2017 to 2020 m2

表4 分區前后高程及坡度變異性Table 4 Variability of elevation and slope before and after the partition

通過分析物聯網傳感器埋設壟線土壤溫度、水分及NDVI隨地形分布規律來探究地形如何影響作物長勢及精準管理分區結果。提取傳感器6月土壤平均溫度及水分（百分比%）、壟線高程及NDVI變化情況，繪制傳感器月均土壤含水量及溫度變化（圖5）、4期NDVI隨地形變化情況（圖6）。由圖5a可以看出4個傳感器分別位于陽坡、坡頂、陰坡、陰坡坡底4個不同坡位。由圖6分區邊界可發現，該壟線上不同時期下精準管理分區邊界均位于陰坡上的某兩處，說明該壟線上不同地形條件下作物長勢情況不同，且不同時期的分區情況相似，也即代表壟線上作物長勢分布情況多年相似。由圖5b及圖6發現，陽坡土壤溫度及土壤水分均最高，4期NDVI的最高值也均位于陽坡上；坡頂土壤溫度較高，土壤水分最低，坡底土壤水分較高但土壤溫度最低，NDVI值自坡頂至坡底逐漸降低。這是因為6月上旬東北玉米處于出苗期，作物長勢受主要受土壤溫度和水分影響，溫度的影響更大，水分次之，溫度高出苗快，土壤水分充足則會加速出苗速度[26-28]，陽坡及坡頂接受太陽輻射較多，溫度較高且土壤水分充足，土壤溫暖而濕潤故作物出苗較快，而地勢低處接受太陽輻射少于陽坡及坡頂，且地勢低易存在積水情況，導致土壤相對冷而濕，作物出苗條件較陽坡及坡頂差。因此地形的起伏造成了土壤溫度和水分的重分配，從而造成研究區連續4年作物生長初期精準管理分區1、2、3級區域分別位于地勢較低且平緩處、稍有起伏及地勢較高處的原因。

3 討 論

東北地區漫川漫崗，田塊內部作物長勢差異明顯，田間精準管理是高效利用耕地的有效途徑。關于精準管理分區輸入量其他學者已做大量對比研究，NDVI也被證明與產量高度相關[14]，是常用的分區數據源，因此本研究選取NDVI作為輸入量具有可靠性，能夠準確反映作物長勢情況。從研究角度來看，目前已有精準管理分區研究多關注分區輸入量的對比及分區方法的技術創新，對于同一田塊分區時空格局及影響因素的探究較少，而本研究對具有地勢起伏特征的同一田塊連續4 a作物出苗期進行分區，并在分區后進行時空格局及影響因素分析，研究側重點較以往對于分區的研究有所不同。Qiu等[29]雖關注了分區格局，但僅進行了定性描述而未實現定量分析，本研究基于分區結果進行不同分區的分級與定義并借助空間轉移矩陣對田塊4年分區時空格局進行定性與定量分析，具有一定的創新性。分析發現，在東北玉米苗期，同一田塊內部作物長勢差異及精準管理分區結果多年相似，而分區后區內作物長勢、地形都趨于均一，說明作物長勢及精準管理分區結果受地形特征影響，這也與劉煥軍等[30]觀點達成一致，同時雷斯越等[31]研究也表明植被長勢與地形具有一定相關性。通過對傳感器所在壟線進行分析發現，地形起伏造成區域太陽輻射及水分分布差異，影響區域內溫度分布[32]，從而影響精準管理分區時空格局，這一結果也與Qiu等[29]的研究結果相吻合。

本文僅探究了作物生長初期的情況，此時作物長勢主要受土壤溫度和水分影響，下一步將重點關注在之后的作物生長階段中其他影響作物生長的因素對作物長勢的影響，以及不同作物生長階段精準管理分區結果特征，深入探究各種微地形特征在水肥運移中的作用，并考慮在不同作物生長階段精準管理分區結果的動態性和時效性，以實現精準管理分區與精準施肥等各種田間變量管理措施的精準對接。高，分區后NDVI變異系數降低了70.690%～76.420%。

2）作物生長初期，在其他管理措施相同以及未受極端天氣影響的情況下，對于同一田塊，多年精準管理分區時空格局相似。

3）地形是影響精準管理分區格局的重要因素，分區后高程及坡度變異系數分別降低42.857%～57.143%、30.723%～34.940%；同一壟線上4期NDVI最高值均位于陽坡，且自坡頂至陰坡坡底NDVI逐漸降低，因為地形影響土壤水分和溫度分布，進而影響作物長勢及精準管理分區結果。研究成果為精準管理分區后與精準施肥、施藥等田間變量管理措施銜接提供時空框架。

4 結 論

本文以典型黑土區友誼農場田塊為研究對象，構建了適于該區的精準管理分區方法，刻畫了多年精準管理分區格局，揭示了寒地漫川漫崗黑土區精準管理分區格局的成因。得出主要結論如下：

1）基于遙感手段劃分精準管理區方便快捷且準確性