基于NDVI遙感反演的半干旱沙區耕地地表溫度異質性研究

周 建，張鳳榮，徐 艷，邱孟龍，謝 臻

（1. 陜西師范大學西北國土資源研究中心，西安 710119；2. 中國農業大學土地科學與技術學院，北京 100193）

0 引 言

半干旱沙區是指在干濕氣候區劃上為半干旱區且土壤質地以沙質為主的區域[1]，荒漠化是該區域最突出的生態環境問題[2]。相關研究表明，自然因素和人類活動因素共同導致了土地荒漠化，但是，各區域引起土地荒漠化的主導因素存在差異[3-4]。就半干旱沙區而言，土地資源廣袤，自然降水為半干旱，種植業與養殖業在經濟上存在互補[5-6]，以上因素使得廣種薄收的土地利用方式是土地荒漠化的主要原因。因此，研究該區域耕地利用特征對實現土地荒漠化的進一步防治具有重要意義。

當前，相關研究基于區域視角，對半干旱沙區耕地面積的變化及其空間差異、驅動機制等進行了大量研究[7-11]，指出在人口增長背景下該區域耕地面積呈現增長的趨勢。在此基礎上，相關文獻對半干旱沙區土地耕作的適宜性及其空間差異[12]、耕地質量的變化[13]、耕地資源的利用特征[14]及如何實現高效利用進行了研究[15]，指出不適宜開發利用為耕地的土地被開發利用引起了土地荒漠化[16]。同時，相關研究對半干旱沙區耕地利用的穩定性及不穩定耕地的類型進行了分析[1,17]。盡管相關文獻對半干旱沙區耕地變化及其利用狀況進行了較多的分析，但是，很少有研究基于耕地地塊的尺度，研究半干旱沙區的耕地利用特征。尤其是在中國土地荒漠化防治呈“整體逆轉、局部擴張”的背景下[18-19]，研究耕地地塊的利用特征，對進一步防治土地荒漠化具有重要參考。同時，當前中國農用地分等均以耕地圖斑作為評價單元[20-22]，研究半干旱沙區的耕地利用特征，也可為該區域耕地分等評價單元的確定提供參考。

基于此，本研究基于地塊尺度分析半干旱沙區耕地地表溫度特征，在此基礎上探討耕地地表溫度特征形成的原因，據此提出進一步防治土地荒漠化、確定該區域耕地分等評價單元的建議。

1 研究區概況

科爾沁左翼后旗（科左后旗，121°30′ E－123°42′ E，42°40′ N－43°42′ N）位于中國面積最大的科爾沁沙地的東南部，屬于半干旱沙區。面積為1.1×104km2，海拔在88.5～208.4 m之間。氣候為溫帶大陸性季風氣候，多年平均降水量為415 mm，主要集中在6－8月份。年平均風速3～4 m/s之間，且大風天氣主要集中在冬、春季節。土壤類型以風沙土為主，占土壤總面積的68.9%，沙性土壤使土地易風蝕沙化。科左后旗地形平緩，固定沙丘、半固定沙丘、流動沙丘相結合，坨甸相間分布是其主要的地貌特點，坨甸相間的微地形導致了自然降水的重新分布，使得土壤水分地形分異顯著。科左后旗屬于國家級貧困縣，農牧業仍為當地農民的主要收入來源。2013年人均土地面積和耕地面積分別為2.81、0.74 hm2，均遠高于中國0.70、0.09 hm2的平均水平。耕地作物以旱地玉米為主，2016年全旗糧食作物播種面積1.89×105hm2，其中玉米播種面積1.62×105hm2，占85.7%。2016牧業年度牲畜頭數304.29 萬頭（只、口）。50年代末期至80年代末期，荒漠化土地面積占土地總面積的比例由 22%增至48%[2]。圖1為2013年科左后旗遙感影像。

圖1 2013年科左后旗遙感影像及其擦除范圍Fig.1 Remote sensing images of Horqin Left Back Banner in 2013 and deleted parts

2 數據與研究方法

2.1 數據與處理

本文所用數據包括遙感影像、土地利用/覆被、大氣水汽含量、地形高程和土壤化驗數據。

利用2013年9月12日（行編號條帶號為030/120）和2013年9月5日（行編號條帶號為030/119）兩幅Landsat－8遙感影像反演科左后旗的地表溫度，兩幅遙感影像覆蓋科左后旗的范圍如圖1。由于2013年9月5日的030/119幅遙感影像存在一定的云量，為避免云對地表溫度反演產生影響，對有云的區域進行了刪除[23]。

2013年土地利用/覆被數據來源于對2013年9月12日（行編號條帶號030/120）、2013年9月5日（行編號條帶號030/119）、2013年7月3日（行編號條帶號030/119）共3幅遙感影像的解譯[5]，土地利用/覆被類型包括耕地、林地、草地、建設用地、沙地、水域。此外，土地利用數據還有2012年科左后旗土地利用變更調查數據。

大氣水汽含量為 NASA網站（https://ladsweb.nascom.nasa.gov/data/search.html）的Modis05_L2數據，利用 HEGTool（Hdf-Eos to GIS Conversion Tool）將Modis05_L2近紅外大氣水汽含量數據的坐標轉換為UTM投影，并將其1 km空間分辨率利用最近鄰法重采樣為30 m。

地形高程為GDEMV2 30 m（http://www.gscloud.cn/）數據，經過鑲嵌、裁剪等操作，獲得科左后旗的DEM。

2015年 10月采集科左后旗邊布拉村耕地地塊不同位置的表層土壤樣品，經實驗室分析獲得土壤有機質含量數據。2014年10月8日至2014年10月16日、2015年7月20日至2015年7月31日對科左后旗進行了有關土地利用、農業生產、牲畜養殖等方面情況的實地調研。

2.2 研究方法

2.2.1 地表溫度遙感反演

對沙區而言，沙的比熱容小，受熱后溫度上升較快。因此，沙區地表覆蓋狀況對地表溫度影響明顯，地表溫度對植被長勢敏感[24-26]。科左后旗玉米生長期為 5月 5日至 9月 17日[27]，本文選用玉米生長末期（2013年 9月12日、2013年9月5日）的兩幅Landsat－8遙感影像反演地表溫度，以更好地體現玉米生長差異。本文利用Rozenstein 等提出的劈裂窗算法（式1）反演地表溫度[28]。

式中 Ts地表溫度，K；A0、A1、A2是由大氣透射率（atmospheric transmittance）和地表比輻射率（land surface emissivity）決定的系數，T10、T11是TIRS傳感器10波段和11波段的亮度溫度。TIRS傳感器10波段和11波段的大氣透射率利用中緯度夏季大氣剖面的透射率公式進行計算。對于TIRS傳感器10波段和11波段的地表比輻射率，首先根據土地利用/覆被類型將地表劃分為水體、自然表面和城鎮3類，水體和沙地（裸土）的地表比輻射率按照宋挺等給定的參數進行取值[29]，城鎮和自然表面的比輻射率按覃志豪等提出的方法進行計算[30]。TIRS傳感器10波段和11波段的亮度溫度是在計算其輻射亮溫的基礎上進行計算。地表溫度的具體計算過程詳見文獻[23]。

2.2.2 耕地地塊地表溫度及其邊緣地表溫度

從2012年土地利用變更調查數據庫中提取耕地，并與 2013年遙感解譯的耕地進行空間疊加，以獲得 2013年耕地地塊邊界。利用ArcGIS軟件的區域統計功能計算各地塊溫度范圍。

在分析耕地地塊邊緣溫度時，為避免耕地地塊面積過小對結果造成影響。首先，選擇已獲得2012年耕地邊界的面積大于8100 m2的2013年遙感解譯的耕地地塊。其次，利用緩沖功能向地塊內部緩沖30 m，以獲得地塊邊緣部分。最后，利用ArcGIS的區域統計功能分析地塊邊緣部分的地表溫度。

2.2.3 高程指數

利用高程指數展示坨甸相間的微地形狀態（式 2）。首先，將科左后旗劃分為2 km×2 km的格網，然后統計每個格網的高程平均值，最后計算高程指數。某柵格高程指數越小，表示該柵格在某一格網中地形越低。

式中EI為高程指數；E為高程，m；ˉE為格網的高程平均值，m。

3 結果與分析

3.1 耕地地表溫度總體狀況分析

行編號條帶號為 030/120部分的耕地地表溫度最小值為294.03 K，最大值為308.61 K，相差14.58 K。就耕地地表溫度的分布而言，地表溫度分布在[299 K，300 K)范圍內的面積最大，占耕地總面積的 26.278%（表 1）；占耕地總面積的比例在10%以上的溫度范圍包括[297 K，298 K)、[298 K，299 K)、[299 K，300 K)、[300 K，301 K)、[301 K，302 K)共 5個溫度范圍，面積占耕地總面積的92.492%，跨度5 K。030/119部分的耕地地表溫度最小值為291.43 K，最大值為310.62 K，相差19.19 K。就耕地地表溫度的分布而言，地表溫度分布在[302 K，303 K)范圍內的耕地最多，面積占耕地總面積的30.582%（表1）；占耕地總面積的比例在10%以上的溫度范圍包括[300 K，301 K)、[301 K，302 K)、[302 K，303 K)、[303 K，304 K)共 4個溫度范圍，面積占耕地總面積的 82.519%，跨度4 K。以上結果表明，耕地地表溫度分布范圍較大。

表1 耕地地表溫度Table 1 Land surface temperature of cultivated land

3.2 耕地地塊溫差分析

就行編號條帶號為 030/120部分的耕地地塊溫差的分布而言，地塊溫差小于1 K的耕地面積占耕地總面積的比例僅為9.687%，地塊溫差小于2 K的耕地面積占總面積的35.656%（表2）。但是，地塊溫差大于2 K的耕地面積占總面積的64.344%，接近耕地總面積的2/3；地塊溫差大于3 K的耕地面積占總面積的36.895%，超過耕地總面積的1/3；地塊溫差大于4 K的耕地面積占總面積的16.841%。就行編號條帶號為030/119部分的耕地地塊溫差的分布而言，地塊溫差小于1 K的耕地面積占耕地總面積的比例僅為3.299%，地表溫差小于2 K的耕地面積占總面積的19.526%。相反，地表溫差大于2 K的耕地面積占總面積的80.474%；地表溫差大于3 K的耕地面積占總面積的57.971%；地表溫差大于4 K的耕地面積占總面積的41.571%。反映出科左后旗耕地地塊溫差較大，且地塊溫差較小的耕地面積占總面積的比例較小，地塊溫差較大的耕地面積所占比例較大。

表2 耕地地塊溫差Table 2 Temperature difference of cultivated land parcel

3.3 耕地地塊邊緣溫度分析

在行編號條帶號為030/120部分面積大于8100 m2的22107塊耕地中，有13268塊耕地的邊緣地表平均溫度較地塊平均溫度高，數量占耕地總塊數的60.02%，面積占耕地總面積的73.71%（表3）。在行編號條帶號為030/119部分的面積大于8100 m2的1644塊耕地中，有866塊耕地的邊緣地表平均溫度較地塊平均溫度高，數量占耕地總塊數的 52.68%，面積占耕地總面積的 53.80%（表 4）。反映出耕地地塊邊緣的地表平均溫度高于地塊平均溫度。

表3 行編號條帶號030/120部分耕地地塊邊緣平均溫度與地塊平均溫度Table 3 Average temperatures at cultivated land parcel margin and cultivated land parcel in row/path 030/120 part

表4 行編號條帶號030/119部分耕地地塊邊緣平均溫度與地塊平均溫度Table 4 Average temperatures at cultivated land parcel margin and cultivated land parcel in row/path 030/119 part

4 討 論

4.1 耕地地表溫度范圍較大原因分析

科左后旗位于半干旱沙區，耕地利用以旱地玉米為主，自然降水量決定了玉米的生長狀況。根據相關文獻，科左后旗自然降水量由東部的 462.8mm下降到西部的358 mm，降水量空間分異明顯[23]，如此便導致了玉米生長空間分異顯著，耕地地表溫度范圍較大。

4.2 耕地地表溫度非均質性分析

坨甸相間是科左后旗主要的地形，坨地為海拔相對較高的高地，甸地為海拔較低的洼地。沙質土壤持水性能差，坨甸微地形對自然降水進行再分配，導致微地形不同部位的土壤水分狀況差異顯著[31]，甸地土壤水分條件相對最多，坨地土壤水分條件相對最少[32]。地形相對越高，作物生長所需水分滿足程度越低，生長越差，地表裸露，地表溫度相對越高。就行編號條帶號為030/120部分的耕地地表溫度與高程指數而言，在[294 K，302 K）范圍內，隨著高程指數的增加，耕地地表溫度逐漸升高（圖2a），此部分耕地面積占030/120部分耕地總面積的93.729%（表1）。就030/119部分的耕地地表溫度與高程指數而言，在[291 K，307 K）范圍內，隨著高程指數的增加，耕地地表溫度逐漸升高（圖 2b），此部分耕地占030/119部分耕地總面積的99.948%（表2）。表明，在坨甸微地形下，隨著相對海拔的升高，耕地地表溫度在升高。坨甸微地形的普遍存在，使得耕地地塊并非完全平坦而是存在一定的高差，造成耕地地塊水分條件和作物生長存在差異，使得耕地地塊溫差較大。

在行編號條帶號為030/120部分當地表溫度大于302 K和行編號條帶號為030/119部分當地表溫度大于307 K時，隨著高程指數的降低，耕地地表溫度反而出現下降的趨勢。此部分耕地占耕地總面積的比例較小，分別為6.271%和0.052%（表1）。這主要是由于大陸性季風氣候條件下坨甸微地形低洼處匯集水分較多，耕地作物生長在一定程度上受瀝澇影響，作物生長較差[1]，植被覆蓋度低，地表溫度較高。

盡管《農用地分等規程（GB/T 28407—2012）》中提出耕地質量的評價單元可以采用疊置法、地塊法、網格法、多邊形法等方法進行確定。但是，目前在農用地分等的實際工作中，耕地質量評價單元的確定基本上采用地塊法，即以土地利用現狀數據庫中的耕地圖斑作為耕地質量的評價單元[20-22]。本研究發現在半干旱沙區，耕地地塊溫差較大，耕地地塊作物生長差異顯著，耕地質量存在明顯差別。如果采用地塊法（耕地圖斑）勢必會消除耕地地塊內部的質量差異，建議在半干旱沙區進行耕地質量評價時采用網格法確定評價單元。

圖2 耕地地表溫度與高程指數Fig.2 Land surface temperature of cultivated land and elevation index

4.3 耕地地塊邊緣效應分析

通常情況下地塊邊緣種內競爭較地塊內部弱，作物生長獲得的養分等資源多，生長狀況一般較好[33]，地塊邊緣溫度較地塊平均溫度低。但是，對科左后旗的實地調研發現，科左后旗耕地承包分配時均是土壤條件較好的耕地，此后農戶普遍以承包地地塊為基礎，開墾周邊的土地以擴大耕種面積[6]。根據耕層土壤有機質分析結果，1號點土壤有機質含量為 44.33 g/kg，2號點為25.92 g/kg。3、4號點分別為24.29 g/kg、23.41 g/kg，5、6號點分別為1.88 g/kg、2.56 g/kg（圖3），反映出耕地邊緣的土壤有機質含量較耕地內部的土壤有機質含量低，即以承包地地塊為基礎向外擴展的耕地土壤質量較差。其次，質量較高的承包地一般位于地形較低的區域[34]，向外開墾的土地地勢較高[35]，耕地的土壤水分條件較差。第三，由于承包地周邊開墾的耕地質量差，農戶對開墾耕地的利用管理更加粗放[6]。以上因素導致了半干旱沙區耕地地塊邊緣作物生長較內部差，地塊邊緣溫度較地塊內部溫度高。由于耕地地塊周邊作物生長差，地表裸露，土地易于風蝕沙化，阻礙了生態環境的進一步恢復。為此，承包地地塊周邊開墾的新增耕地應當進行退耕，以進一步防治土地荒漠化和恢復生態。

圖3 耕地表層土壤采樣點Fig.3 Soil sampling point at cultivated land surface

4.4 遙感反演地表溫度的準確性及其反映作物生長狀況的有效性分析

本文分析了耕地的地表溫度特征及其形成原因。在遙感反演地表溫度的準確性方面，本文利用遙感反演的方法計算的是科左后旗2013年的地表溫度，由于無法對2013年的地表溫度進行實測，使得無法利用實測地表溫度對遙感反演的地表溫度進行驗證。但是，相關研究表明本文所采用的劈裂窗算法反演地表溫度的精度較高，Rozenstein等采用本方法反演的地表溫度的均方根誤差為0.93 ℃[28]。宋挺等對利用Landsat－8遙感影像反演地表溫度的精度進行分析，表明Jiménez-Mun?oz劈裂窗算法和Rozenstein劈裂窗算法的精度較高，誤差僅為0.7 K左右[29]。遙感反演地表溫度反映作物生長差異的有效性方面，首先沙區土壤質地為沙質，受熱后地表溫度上升快。其次植被生長狀況對沙區地表溫度上升具有顯著影響[24-26,36]。植被覆蓋度高的土地較植被覆蓋度低的土地地表溫度上升小。相關文獻對科爾沁沙地坨甸地形耕地作物生長狀況與耕地地表溫度進行研究，表明隨著坨甸地形的升高土壤水分條件變差，耕地作物生長也在變差，耕地地表溫度隨之升高[23]。因此遙感反演地表溫度可以表征半干旱沙區耕地作物的生長差異。但是，一方面本文所采用的Landsat－8遙感影像的獲取時間為2013年，另一方面由于未在遙感影像獲取的時間進行作物生長狀況的采樣，導致無法對遙感反演地表溫度與作物生長狀況的關系進行分析，未來將對沙區遙感反演的地表溫度與作物生長狀況之間的關系進行深入研究。

5 結 論

1）就耕地地表溫度的總體狀況而言，耕地地表溫度范圍較大。在行編號條帶號為030/120部分，占耕地總面積比例在 10%以上的溫度范圍包括[297 K，298 K)、[298 K，299 K)、[299 K，300 K)、[300 K，301 K)、[301 K，302 K)，跨5個溫度范圍，面積占耕地總面積的92.492%；在行編號條帶號為030/119部分，占耕地總面積的比例在10%以上的溫度范圍包括[300 K，301 K)、[301 K，302 K)、[302 K，303 K)、[303 K，304 K)，跨4個溫度范圍，面積占耕地總面積的82.519%。

就耕地地塊溫差而言，耕地地塊溫差較大，且地塊溫差較小的耕地面積占總面積的比例較小，地塊溫差較大的耕地面積所占比例較大。

2）耕地地塊邊緣的地表平均溫度高于地塊平均溫度。

3）基于以上研究，認為在半干旱沙區進行耕地質量評價時，評價單元的確定宜采用網格法；建議對耕地地塊邊緣的耕地進行退耕，以進一步防治土地荒漠化和恢復生態。