黑土區田塊尺度下地形影響作物長勢機理分析

劉煥軍 鮑依臨 孟祥添 徐夢園 錢 蕾 張新樂

(1.東北農業大學公共管理與法學院， 哈爾濱 150030； 2.中國科學院東北地理與農業生態研究所， 長春 130012； 3.中國農業大學土地科學與技術學院， 北京 100094)

0 引言

了解農作物空間變異規律可以提升糧食產量的模擬精度。根據田間作物長勢，有針對性地進行施肥管理，可以同時改善土壤質量、降低土壤侵蝕帶來的影響，對農業可持續發展具有重要意義。隨著遙感技術在農業領域的應用與普及，歸一化植被指數(Normalized differential vegetation index，NDVI)在農業生產中的作用得到了充分的肯定[1]，NDVI與土壤有機碳和總氮含量的富集程度呈正相關[2]，并作為重要指標用于生物量分布及產能評估[3]。NDVI的動態變化不僅受氣象驅動變量(如溫度、陽光和降水)的影響，不同的環境因素(如海拔、表征水分脅迫、農業管理)也與之相輔相成[4]。已有研究表明，地形因素對作物的影響極為重要，土壤性質和地形(坡度和海拔)的綜合作用可較大程度影響植被生物量的變化[5]。因此，探究NDVI的空間變異規律對于提供作物和農場狀況的詳細信息、指導田間管理具有重要意義[6]。

諸多學者已開展了地形對植被長勢作用的相關研究，其研究尺度更傾向于大范圍的復雜地形[7]。高志強等[8]采用相關分析方法探究中國植被指數的時空變化，水、熱和地表植被覆蓋類型對植被指數影響和控制的主導地位因時因地而不同，只能從概念上定性地描述植被指數和3個主導因子間的關系。魏偉等[9]運用疊置分析法研究了高寒草原區不同海拔和地表覆蓋下的NDVI變化情況，表明NDVI與地形具有很強的相關性。縮小尺度范圍可以將地形對作物長勢的影響從定性到定量描述。張少良等[10]利用USLE模型預測黑土區域土壤侵蝕量的空間分布格局，較好地模擬了不同坡位的土壤流失狀況。在此研究的基礎上，梁佳輝等[11]在同一區域選定4條剖面，探究不同耕層的土壤速效磷在水平和垂直方向的分布規律，證明了侵蝕、沉積與施肥等管理因素對不同坡位各土層剖面速效磷變異性的影響。基于上述研究可提出假設：地形對NDVI的作用程度與研究區尺度有關，縮小研究范圍有助于減少空間的不確定因素，更好地表達地形與NDVI之間的關系。

東北黑土區是我國重要的商品糧基地，其產能與我國糧食安全密切相關，但土壤侵蝕對作物生產構成嚴重威脅[12]。現階段關于田塊尺度的地形與作物長勢差異的研究并不多見，田塊內部作物種植結構單一，且田塊尺度與實際管理相對應更有利于指導農業生產。因此，本研究獲取高精度地形信息，分析同一地塊不同年份的作物長勢狀況，分別從坡型、坡度與曲率分析不同地形條件對于作物長勢的影響機理，以探究田塊尺度下地形對作物長勢的影響規律，以期為精準農業管理提供理論支撐。

1 材料與方法

1.1 研究區概況

研究區位于黑龍江省中部海倫市東興農機合作社，地處松嫩平原東北端，小興安嶺西麓，平均海拔239 m，地塊中心處47°25′N、126°55′E，面積約19.35 hm2，屬寒溫帶大陸性氣候，冬季寒冷干燥，夏季溫熱多雨，年均氣溫1.5℃，年均降水量530 mm，年均有效積溫2 450℃，年均日照時數為2 600～2 800 h，無霜期為125 d。種植作物主要以大豆、高粱輪作為主，一年一熟制，種植期間采取統一的田間管理方式。地勢從東北到西南依次呈梯形逐漸降低，地貌屬于小興安嶺山地向松嫩平原的過渡地帶。土壤類型為典型黑土中的黃土質黑土，漫川漫崗地帶，地形起伏較大，作物長勢空間差異性顯著，研究區位置如圖1所示。

1.2 數據源與研究方法

1.2.1地形數據的獲取和預處理

2016年7月15日，使用定位精度厘米級的海星達iRTK2對研究區地塊進行實地測量，得到193個精確的坐標和高程點，利用ArcGIS 10.2將高程數據生成TIN圖層，將TIN圖層轉成空間分辨率為4 m的數字高程模型(Digital elevation model，DEM)柵格數據。2018年，研究區的中部位置種植玉米，故壟線的提取分布在東西兩側，沿等高線間隔由西至東提取15條斷面壟線，壟線方向均由北向南(圖2)。依據地形起伏的相似性，將15條壟線劃分為5類坡型，依據陰坡坡型形狀特征，可將坡型依次定義為：復合坡、坡腰上部微凸型、坡腰上部中凸型、坡腰強凸型、坡腰中凸型，分別以A(壟線1～2)、B(壟線3～6)、C(壟線7～8)、D(壟線9～12)、E(壟線13～15)表達，構建直線以連接坡底與坡頂，以便更好地對比各類坡型的特征。

圖1 研究區位置Fig.1 Study fields location

圖2 斷面壟線空間分布圖Fig.2 Ridge spatial distribution map

1.2.2SPOT影像獲取及處理

研究表明，作物收獲前42～56 d的NDVI信息有助于實現高精度的產量預測[13]。因此，使用大豆生長期2016年7月19日SPOT-6遙感數據(多光譜波段的空間分辨率為6 m)，在ENVI 5.1中對數據進行大氣校正、幾何校正，在ArcGIS 10.2中按研究區范圍進行裁剪并提取NDVI。

1.2.3無人機影像獲取及處理

圖3 高程及NDVI空間起伏度示意圖Fig.3 Spatial relief diagrams of DEM and NDVI

本實驗采用大疆M600Pro型六旋翼飛行器作為遙感平臺，搭載MicaSense RedEdge TM3型多光譜相機，該相機有紅、綠、藍、紅邊、近紅外5個波段，成像分辨率為1 280像素×960像素，配備光強傳感器和一塊固定反射率校正白板。在飛行過程中，光強傳感器可以校正太陽光線變化對影像所產生的影響，白板具有固定的反射率信息，可以對影像進行校正，生成反射率影像。數據采集時間為2018年7月19日，天氣晴朗無云，風速小于1級，適合飛行。將相紙交叉擺放，以交叉中心點作為影像事后處理的空間校正點，利用iRTK2采集地面控制點坐標。為保證影像的完整和準確，通過地面站設置飛行的航向重疊度為80%，旁向重疊度為79%，主航線角99°，邊距-1.8 m，航線9條，飛行高度為110 m，空間分辨率可達0.074 m，航拍完成后采用自動返航方式降落。對無人機采集的影像進行篩選，去除姿態角異常、成像存在問題的影像。將篩選后的原始影像和POS數據輸入Pix4D mapper軟件，利用軟件的農業多光譜影像模塊進行影像拼接。軟件運行后自動生成連接點并利用內封空三計算模型計算像元位置以及參數，經點云加密后，生成數字表面模型，在此基礎上，通過數字微分糾正對影像進行正射校正和鑲嵌，最終得到帶有地理坐標的正射影像和數字表面模型(圖3)。

由于空間分辨率不同，在ENVI 5.1中采用重采樣方法減小尺度效應對地形及NDVI的影響。

1.2.4地形因子的計算

坡度：坡面的垂直高度與水平方向距離之比，反映地表單元陡緩程度。

曲率：通過微分定義，計算一定距離上坡度的變化，以表明曲線偏離直線的程度。

基于DEM數據獲取坡度及曲率數據，利用ArcGIS 10.2軟件3D Analyst模塊slope及 curvature工具生成。

坡型凸凹程度：以復合坡(A)為例，定義多種指標對地形進行表達(圖4)。

圖4 坡型表達因子示意圖Fig.4 Slope type expression factor diagram

直線C距離：構造直線C以連接坡底與坡頂，反映坡長在單位距離內高程的變化，C由高程差a與距離差b通過勾股定理求取。

凸凹距離：應用角檢測法選取凹坡與凸坡內距離直線C最遠的點，該點到直線C的豎直距離以L1、L2表示，該點到直線C的垂直距離定義為坡距，分別以m1、m2表示(圖4)。以凹坡為例，y、y1分別為坡底至坡頂、坡底至凹陷程度最大處的斜率，則

x=y1D1

(1)

x+L1=yD1

(2)

y1D1+L1=yD1

(3)

L1=D1(y-y1)

(4)

L1cosθ=m1

(5)

根據同位角原理，可得

m1=L1sin(90°-θ)

(6)

式中x——坡底至凹坡最低處距離

D1——陰坡凹坡最低處所對應的直線距離

θ——坡頂與坡底夾角

凸凹面積S：反映坡型凸起或凹陷的程度。將凸起范圍按垂直x軸的方向分成n個小梯形，依據黎曼積分原理，對函數進行積分，當n→+∞時，Δx的最大值趨于0，將所有梯形面積相加即矩形面積和，計算式為

S=f(x1)Δx1+f(x2)Δx2+…+f(xn-1)Δxn-1

(7)

因此，凸凹程度由單位直線內的凸凹面積(S/C)與凸凹距離(M/C)表達，S/C或M/C越大，說明坡型變化程度越大。

1.2.5植被長勢分析

計算壟線NDVI的平均值估測作物整體長勢，NDVI的標準差反映作物長勢差異；依據等距原則，將不同空間分辨率的信息進行重采樣處理，計算不同坡面(陰坡、陽坡)NDVI標準差并與地形因子(斜率k、凸凹面積S、凸凹距離M、凸凹面積與長度比S/C、凸凹距離與長度比M/C)進行對比，確定用于表達NDVI變化的最適指標。

1.2.6Logistic回歸模型

Logistic線性回歸模型可以反映因變量影響自變量發生變化的概率。本文選取二分類Logistic進行回歸模擬，探究研究區內地形因子(DEM、坡度、曲率)導致NDVI發生變化的作用程度。將高于壟線NDVI均值的柵格標記為1，其余為0，將因變量的范圍限定在[0，1]范圍內，數據以7∶3的比例分為建模集和驗證集，進行建模及精度驗證。

1.2.7數據處理

影像處理及重采樣、地形因子的相關計算、相關性分析、多元逐步回歸及制圖分別在ArcGIS 10.2、ENVI Classic、Matlab、SPSS stastic22、R Studio、OriginPro中實現。

2 結果與分析

2.1 坡型與NDVI變異規律

2.1.1不同類坡型分析

由表1可知，坡型凸凹面積、凸凹距離、單位直線內凸凹面積、單位直線內凸凹距離均與對應壟線上NDVI的變化趨勢一致，運用上述指標表達地形凸凹程度可行。坡型凸凹程度越大，NDVI標準差越大。值得注意的是，A為復合坡，其復雜的坡型結構使其較其他單一坡型而言，NDVI空間差異性更為顯著。

NDVI的標準差與M/C指標在0.01水平上呈現極顯著相關(表2)，陰坡相關系數R達0.992，陽坡相關系數R達0.997，達極顯著相關。因此，該指標可以更好地反映坡型的凸凹程度，并以此判斷NDVI空間變異是否顯著。

2.1.2同類坡型植被長勢分析

DEM柵格數據與2016年NDVI信息的空間分辨率一致，為了探究坡型差異對NDVI的影響，繪制2016年15條壟線的NDVI隨DEM變化曲線(圖5)。

表1 不同類坡型凸凹程度與NDVI標準差關系Tab.1 Relationship between concave and convex degrees of different slopes and NDVI standard deviation

表2 NDVI標準差與坡型凸凹指標的相關系數Tab.2 Correlation coefficient between NDVI standard deviation and slope concave and convexity

圖5 NDVI隨高程變化曲線Fig.5 Curves of NDVI varied with DEM

對比發現，同類坡型內NDVI的變化趨勢一致，且在坡腰處NDVI的變化相似。同時，為探究同類坡型內不同凸凹程度對NDVI的影響，提取每條壟線的坡型信息。由于陽坡受到的太陽輻射量較強，光熱條件好，因此受到更多的氮素脅迫和水分脅迫[14]，增加了影響作物生長的不確定性，故僅對陰坡進行分析，結果見表3。由表3可知，同類坡型中，陰坡坡型凸凹程度與NDVI的標準差大體成反比。通過對B類和D類坡型分析可知，在同類坡型內存在閾值，決定NDVI標準差與凸凹程度的相關性。

2.2 坡度及曲率對NDVI的影響

2.2.1Logistic模型預測結果

Logistic模型是定量評價模型，預測結果中，ROC(Receiver operating characteristic)曲線用于反映擬合數據與實測數據之間的關系，曲線縱軸為真正類率，即作物長勢的優劣；橫軸為假正類率，即對長勢的錯誤判斷，ROC曲線越接近縱軸表示模型的模擬效果越優異；AUC(Area under the ROC curve)是ROC曲線與坐標軸圍成的面積，是度量分類模型性能的標準，其值越接近1，模型的分類精度越高[15]。以地形因子(DEM、坡度、曲率)為輸入量評估NDVI，結果見圖6。由圖6可知，大豆NDVI的AUC為0.737，高梁NDVI的AUC為0.692，計算結果表明，研究區內的NDVI受地形影響較為顯著，且大豆對于地形變化更敏感。

2.2.2一定坡度范圍內作物長勢分析

由圖7可知，坡度在±0.03范圍內作物長勢相對均一，因此，選定0.03為閾值進行分析。由圖8可知，坡度在±0.03范圍內NDVI的變化與作物的類型相關，對于大豆而言，坡度在±0.03范圍內NDVI均值高于陰坡均值、陽坡均值以及整條坡均值，且該部分的作物長勢相對均一；相比整條壟線，該部分的坡度變化較為穩定[16]。而高粱對于地形的敏感性相對較低，在研究區西部NDVI均值處于陰坡與陽坡均值之間，東部作物的空間分異規律與大豆相同。

2.2.3高坡度對作物長勢的影響

坡腰下部至坡底處坡型陡峭，植被長勢差異較大。根據研究區坡型特征，選取坡度絕對值大于0.04部分作物長勢進行分析，為探究不同年份、不同作物類型植被長勢的共性規律，計算兩期NDVI的相關系數，結果見表4。由表4可知，超半數壟線的兩類植被長勢在0.01水平上極顯著相關，以坡度絕對值高于0.04的坡度均值與對應直線的比例為自變量，NDVI標準差為因變量，構建多元逐步回歸模型，結果見圖9。由圖9可知，坡度絕對值高于0.04的坡度均值與其直線的比例可以解釋大豆NDVI決定系數0.965 2、高粱NDVI決定系數0.888 3的空間變異性，由于作物種類及生長所需環境不同，線性擬合的結果存在差異，大豆的生長更易受到地形的影響，這也進一步印證了Logistic回歸的結果。

2.2.4NDVI極值所處位置分析

曲率可以在一定程度上反映坡面微地形變化，依據坡型位置特征將整個坡面分成4部分：陰坡坡底至陰坡坡腰、陰坡坡腰至坡頂、坡頂至陽坡坡腰、陽坡坡腰至陽坡坡底。每10個像元中選出1或2個NDVI的極值點并標記對應曲率，發現NDVI的極值(某一階段極大值或極小值)易發生在曲率為0及其相鄰像元位置。因此，將曲率無限接近于0及其鄰域柵格設定為符合，由于分辨率不同，NDVI選取的個數存在差別，統計結果見表5。由表5可知，2018年NDVI極值的總個數約為2016年選取個數的2倍，通過對比總符合率可知，2018年NDVI總符合率比2016年多18.99個百分點，像元尺度對評價結果具有較大的影響。

表3 陰坡NDVI標準差與坡型凸凹指標分析Tab.3 Analysis of NDVI standard deviation and slope convex index

圖6 地形因子對NDVI模擬結果Fig.6 Simulation results on terrain factor and NDVI

圖7 坡度與NDVI變化曲線Fig.7 Changing curves between slope and NDVI

圖8 NDVI均值、標準差與坡位的關系Fig.8 Relationship between NDVI mean,standard deviation and slope position

3 討論

探究地形對作物長勢及水肥運移的規律有助于實現精準施肥及追肥。趙越等[17]分析了泰和縣耕地土壤養分空間分布及影響因素，指出地形、成土母質等結構性因素是影響土壤養分含量的關鍵因素；李亨偉等[18]對丘陵區微地形下的土壤有機質(Soil organic matter，SOM)的空間變異性進行了研究，將丘體劃分為丘頂、丘腰上部、丘腰中部、丘腰下部、丘腳5種類型并對不同位置SOM含量進行了分析，結果表明，不同坡型下凹型坡SOM含量高于凸型坡，且SOM含量隨坡度的增大而減小。由此可見，土壤養分含量的變化與地形密切相關，且地形對植被長勢的影響并不隨著空間尺度的減小而降低。

黑土區為我國重要的商品糧基地，確保田塊內部的精細管理有助于產量的提升。因此，本研究在田塊尺度上開展了不同坡型間的對比，明確微地形對作物長勢的作用，為使研究成果更具魯棒性，選取不同年份的NDVI綜合分析，使結果在不同作物類型下適用性更廣。

不同坡位的土壤侵蝕程度存在差異，地表水流輸沙能力較強，下坡移動降雨的土壤侵蝕大于上坡[19]，本研究區坡度在±0.03范圍內處于坡頂至坡肩部位，相比坡肩以下部位，該區域溫度、太陽輻射較為均一，地表徑流流速小，水肥保持效果好，因此作物長勢好。由于高粱對微地形的敏感性較低，而研究區東側侵蝕溝溝面較寬，溝底深，溝內幾乎無作物生長，因此NDVI均值低，且作物長勢差異性顯著。在相同土地利用類型下，隨著坡度的增加并大于一定范圍時，土壤侵蝕顯著增加[20]；同時，土壤表面粗糙度作為影響水文侵蝕過程的重要因素之一，其隨坡度的增加呈指數增加[21]。因此，研究區內坡度絕對值高于0.04的位置土壤表面粗糙且土壤侵蝕嚴重，不利于作物的生長。在進行變量施肥時，需分析地塊內部坡度對植被長勢的影響，對于坡度起伏較大區域進行額外的關注，有利于合理分配人力物力投入，節約管理成本。

同類研究中，不同坡型及坡度分析大多是建立在連續坡面，而微地形對于土壤養分的空間分布與作物長勢造成很大的影響[22-23]。2018年NDVI極值位置的總符合率高于2016年，與無人機影像較高分辨率有關，多項研究表明，衍生屬性的質量很大程度上取決于元數據的來源、規模和分辨率，研究的“最佳”窗口與該地區滑坡的平均大小直接相關[24]，因此，較粗的單元可能無法識別最小滑坡區域[25]，由于本研究是基于田塊尺度的微地形的探究，故空間分辨率越高，越有利于把握微地形對作物長勢的作用。

本文對研究區內的坡型進行定義，并采用多種地形因子對NDVI空間變異性進行描述。同時，劃定了不同坡度閾值以分析及模擬植被長勢的空間變異程度，可為田塊尺度內坡型的研究提供一定的思路。然而研究區內土壤類型單一，不同區域的成土母質差異明顯，因此坡度閾值的選取可能不適用其他土壤類型。因此，在未來的研究中，應綜合分析不同類型土壤下的空間變異特征，使該類研究能夠更好地為精準農業管理提供服務。

4 結論

(1)綜合考慮了坡型的形狀特征及凸凹程度，對黑土區微地形坡型進行了更細致的定義。不同類坡型凸凹程度越大，NDVI的空間變異性越大；同類坡型中NDVI的變化趨勢一致，陰坡NDVI的空間變異程度與坡型凸凹程度大體呈負相關。復合坡中NDVI差異更為明顯，這是由其特殊的坡型決定的，復合坡型中陽坡凸起與凹陷的交界處，NDVI表現為大幅度驟減。

表4 不同壟線兩期NDVI相關系數Tab.4 Pearson correlation values between NDVI in two stages and different lines

圖9 NDVI標準差實測值與預測值關系Fig.9 Relationship between NDVI standard deviation measured value and predictive value

表5 NDVI極值位置符合設定的統計Tab.5 NDVI extreme position matched set statistics

(2)坡度在±0.03范圍內的作物長勢好，且空間變異性小。坡度絕對值大于0.04的坡度均值與其直線距離比例可以解釋大豆NDVI決定系數0.965 2、高粱NDVI決定系數0.888 3的空間變異性。

(3)研究田塊內微地形對作物長勢影響顯著，NDVI的極值易出現在曲率為0及相鄰像元處，且總符合率與元數據尺度呈正相關。