基于HJ 衛星數據的土壤含水量反演及其旱情預測

熊世為， 李衛國， 賈天山， 胡姍姍， 景元書

(1.滁州市氣象局，安徽 滁州239000;2.江蘇省農業科學院農業經濟與信息研究所，江蘇 南京210014;3.南京信息工程大學應用氣象學院，江蘇 南京210044)

農業上將長期無雨或少雨，使土壤水分不足、作物水分平衡遭到破壞而減產的災害定義為干旱。干旱是中國常發的災害之一，嚴重影響社會經濟發展。以2010 年西南旱災為例，此次災害共造成5.03 ×106hm2面積的農作物受災，其中絕收面積達1.12 ×106hm2，直接經濟損失達2.366 ×1010元［1］。蘇北地區是中國重要的糧食生產基地，由于其特殊的地理、氣候條件，冬春季常常遭受干旱影響，因此，實施及時的干旱監測對保護該地區農業生產以及發展精細化農業具有積極意義。傳統的干旱監測是基于臺站式的“點上”監測，不僅耗時、費力，而且難以獲得災情發展的“面上”信息，在時空上具有局限性。作為一種能夠快速獲取大面積地物光譜信息的技術，遙感在干旱監測方面已經得到了廣泛應用［2-3］。目前常用的基于衛星遙感的干旱監測模式主要有熱慣量法、植被指數法和地表溫度法等，每種方法都存在各自的優缺點，如熱慣量法物理機制明確，但僅適用于裸土或者植被生長早期［4］;植被指數法在高植被覆蓋區具有較好的應用效果，卻容易受其他環境因素干擾［5］;地表溫度對土壤水分變化比較敏感，但反演過程復雜，且精度有待提高［6］。耦合多種模式的綜合模型已成為該領域的發展趨勢，其中NDVI-Ts特征空間法耦合了歸一化植被指數(Normalized difference vegetation index，NDVI)和地表溫度(Temperature of surface，Ts)，物理意義明確，該方法衍生出來的植被供水指數(Vegetation supply water index，VSWI)、作物水分虧缺指數(Crop water stress index，CWSI)、植被溫度狀態指數(Vegetation temperature condition index ，VTCI)、溫度植被干旱指數(Temperature vegetation dryness index，TVDI)［7］都已經在干旱監測研究中得到了較好的應用。

HJ-1A/B 系列衛星是中國自主研發的用于環境與災害監測的小衛星，兼具MODIS 數據重返周期短及TM 數據空間分辨率較高的優點，已廣泛應用于作物長勢監測［8］、產量估算［9］等方面，但在大范圍、業務化監測(尤其是干旱監測)的應用仍顯匱乏［10-11］。本研究通過構建基于HJ 衛星數據的NDVI-Ts 特征空間，提取溫度植被干旱指數并結合實測數據將遙感指數轉化為土壤相對濕度，并探討利用HJ 數據的TVDI 法在蘇北地區干旱監測中的適用性，以期為國產衛星數據在資源與環境監測及評估方面的深入應用奠定基礎。

1 材料與方法

1.1 試驗數據的獲取

本試驗研究區為江蘇省宿遷市，介于33°80' ～34°25' N、117°56' ～19°10' E 之間，年均降水量為892.3 mm，受季風影響年際變化不大，但分布不均，易形成春旱、夏澇、秋冬干燥天氣。本研究采用手持式GPS 定位儀在宿遷市區、沭陽縣、泗陽縣以及泗洪縣各選取10 個具有代表性的樣點并進行取土采樣。樣點的選定原則是大面積成片的田塊，地表植被單一(調查期間種植小麥)，在距離道路、水體等至少100 m 的位置取土，樣點間距5 km。調查期間發現部分田塊已出現龜裂、小麥葉片萎蔫等旱象。

采用稱質量法獲取土壤含水量，即土壤水分占土壤干質量的百分率。用0.01g 精度的天平稱取土樣的質量，記作土樣的濕質量(M)，在105 ℃的烘箱內將土樣烘6 ～8 h 至恒質量，然后測定烘干土樣，記作土樣的干質量(Ms)，則土壤含水量(W)可表示為:

田間持水量室內測定采用環刀法。土壤相對濕度即土壤含水量與田間最大持水量的比值:

式中R為相對濕度，W為含水量，W0為田間最大持水量。

取土深度為0 ～10 cm、10 ～20 cm 和40 ～50 cm。為方便討論，將各層次的含水量記為10 cm、20 cm、50 cm 處的土壤含水量。本研究共進行4 次地面取土試驗，并選取4 景與之同步的HJ 衛星數據(表1)。

遙感影像預處理過程包括:首先采用校正過的Landsat/TM 影像對HJ 衛星的CCD 和IRS 影像進行幾何精糾正，保證平均誤差在1 個像元以內;再對CCD 和IRS 影像進行輻射定標，將灰度值圖像轉化為具有物理意義的輻亮度圖像，定標公式及系數均來自影像頭文件;最后采用具有較高波譜還原精度的中光譜分辨率大氣輻射傳輸模式(MODTRAN)對CCD 影像進行大氣校正，以獲取真實地表反射率。

1.2 歸一化植被指數(NDVI)和地表溫度(Ts)提取

NDVI被定義為近紅外波段與紅光反射率之差與兩者反射率之和的比值，表示為［12］:

式中，ρR、ρNIR分別表示為紅光、近紅外通道的反射率值。

由于IRS 數據只有一個熱紅外通道(10.5 ～12.5 μm)，與TM 數據相似，因此可以采用單窗算法進行地表溫度的反演［13］:

式中，TS是地表溫度，TB為地表輻射亮度對應的亮溫，有效中心波長λe=11.511 μm［14］，常數C =1.438 768 69 ×10-2mK，ε 為地表比輻射率。其中，TB可以根據Plank 公式推導［15］:

式中，B為IRS 數據熱紅外通道輻射亮度，K1、K2均 為 常 數，對 于 IRS 數 據:K1= 579.20 W/(m2·sr·μm)，K2=1 245.58 K［13］。地表比輻射率ε 根據覃志豪等［16］提出的NDVI閾值法計算。反演的地表溫度參照文獻［17］進行詳細的驗證和正確性討論。

1.3 NDVI-Ts 特征空間與溫度植被干旱指數(TVDI)

土壤濕度與地表溫度(TS)密切相關，根據地表能量守恒，土壤蒸發、冠層葉片蒸騰作用越小，帶走潛熱能量就越小，地表感熱能量大，TS就高;反之，蒸發、蒸騰作用大，潛熱增加，感熱能量降低，導致TS低。而土壤濕度是影響土壤蒸發及冠層蒸散阻抗的重要因素。研究結果表明，在不同土壤表層含水量和地表覆蓋條件下，TS和NDVI的散點圖呈現出一種三角形空間，在此基礎上定義TVDI為［17］:

式中，TS是特征空間內某給定像素的地表溫度;NDVI為該像素的歸一化植被指數;TSmax和TSmin分別表示該NDVI對應的干、濕邊上的地表溫度的最大值和最小值;a1、b1和a2、b2分別為干、濕邊方程的斜率和截距，可通過線性擬合獲取。因此TVDI是一個歸一化值，介于0 ～1，其值越高表示土壤濕度越低，反之土壤濕度越高。

2 結果與分析

2.1 特征空間構建

根據上述方法進行歸一化植被指數(NDVI)和地表溫度(TS)計算，并構建NDVI-TS特征空間。由于NDVI＜0 對應地表類型為水體(可以認定含水量為100%)，在特征空間構建時不予討論。因此，本試驗在NDVI大于等于0 的范圍內，以0.01 為步長，找出各NDVI值對應的最大地表溫度TSmax和最小地表溫度值TSmin，得到研究區各時相的NDVI-TS特征空間(圖1)。可以看出，TSmax隨著NDVI的增大(NDVI＞0)先有所升高，再逐漸降低，但總體上看是隨著NDVI增大呈減小趨勢;TSmin隨NDVI的增大逐漸升高，但升高的幅度沒有TSmax降低的幅度大。總體而言，NDVI與TSmax和TSmin基本呈線性關系(表2)，隨著植被覆蓋程度增加，地表溫度最大值逐漸減小同時地表溫度最小值逐漸升高并交匯于一點，也證明了TS與NDVI之間確實存在一種三角形的特征空間關系。其中TSmin的變化幅度雖然沒有TSmax的大，但并不平行于NDVI值。不同時相特征空間的干、濕邊的離散程度有所區別，分析原因，應該是不同NDVI值對應的象元量差別導致的。

2.2 特征空間干、濕邊確定

前人在進行類似研究時，大都是直接將特征空間所有非水體象元(NDVI值為0 ～1)對應的地表溫度最大值與最小值進行干、濕邊擬合的。但在NDVI的某些范圍內，其對應的地表溫度最大值、最小值并不符合理論模型的變化趨勢，如本研究中在NDVI處于0 ～0.2時，其對應的地表溫度最大值并沒有隨NDVI的升高而降低。地表溫度最大值應隨NDVI增大而線性減小，這是建立在假設植被覆蓋度與NDVI呈線性相關的基礎上的，但實際情況并非如此。有研究結果表明，當植被覆蓋度不足15%時，NDVI對植被覆蓋度的指示作用很低［18］，也就是說，在植被覆蓋低的區域，NDVI夸大了植被信息，因此本試驗在進行特征空間干、濕邊確定時，對NDVI的范圍進行優化處理，選擇NDVI≥0.2 的植被覆蓋地像元，舍棄NDVI介于0 ～0.2 裸地巖石等低植被覆蓋區域。該方法有別于常規的考慮所有非水體象元，這樣做既合理舍棄了NDVI指示性低的區間，又有效提高了高植被區的反演精度。

表2 各時相特征空間干、濕邊方程Table 2 Dry and wet edge equations of feature space of each time phase

2.3 TVDI 與實測土壤水分含量的關系

根據實測土壤含水量數據，分析TVDI與不同深度土壤相對濕度之間的關系。圖2 反映了2012年3 月26 日各土層實測結果與提取的TVDI線性關系，可以看出，各層土壤相對濕度與TVDI均存在一定的線性負相關關系(其他時相結果類似，圖略)，即隨著TVDI指數增加，土壤相對濕度呈逐漸減小趨勢，這與TVDI的原理一致。此外，圖中數據點相對比較離散，經分析有兩類因素導致誤差，第一類是空間性因素，包含兩點:(1)地面試驗每個樣方采集的五點并不能完全代表樣方的平均水平;(2)幾何校正誤差導致像元與地面樣方不能完全匹配。第二類是時間性誤差，本研究地面調查數據都是兩天內完成的，而遙感影像獲取的是某一時刻的土壤水分信息，二者在時間尺度上不能完全匹配。進一步分析各層實測土壤相對濕度與TVDI指數的相關性，并進行顯著性檢驗(表3)，可以看出，TVDI與每個時相各層實測土壤相對濕度的相關性大多達到顯著或極顯著相關的水平(2012-05-17 和2012-05-28 2 個時相的40 ～50 cm土層除外)，說明采用TVDI作為該區域土壤水分含量的監測指標是合理的。

圖2 TVDI 與不同深度土壤相對濕度的關系Fig.2 The relationship between TVDI and relative soil water content in different depths

從不同土層深度上看，0 ～10 cm 深度土層實測含水量與TVDI的相關系數在-0.45 ～-0.60，10 ～20 cm 土層深度實測含水量與TVDI的相關性則更高，為-0.650 ～-0.800，而40 ～50 cm 土層相對濕度與TVDI的相關性最差。這是由于光譜輻射的反射面在淺層地表，可見光、近紅外波段的輻射能量很難到達50 cm 的土層;由于大氣、植被、人類活動等各種因素的影響，表層土壤(0 ～10 cm 土層)一般都是以疏松、充滿間隙的物理形態出現，光譜透射率高，所以TVDI對10 ～20 cm 土層水分含量的指示性優于0 ～10 cm 土層。

2.4 旱情監測效果

盡管TVDI與實際土壤水分含量不能完全匹配，但二者的相關程度較高，尤其是與10 ～20 cm 土層相對濕度的相關系數達到0.8，因此，以TVDI作為淺層(10 ～20 cm)土壤水分含量的監測指標是可行的。根據前文10 ～20 cm 土層相對濕度與TVDI的擬合方程，將基于HJ 遙感影像的TVDI 圖轉化為相對濕度分布圖，在此基礎上根據中國氣象局農業氣候中心提出的旱情等級劃分標準，得到研究區冬小麥淺層土壤(10 ～20 cm)旱情等級分布信息圖(圖3)。

4 個時相發生最重的干旱等級為中旱，雖然沒有出現重旱，但地面試驗期間發現部分麥田土壤出現干裂、植物葉片萎蔫的現象。統計各等級旱情占全市面積的比率(表4)，4 個時相中旱情發生最嚴重的是2012 年5 月28 日，該時相中旱情發生面積占全市總面積的16.57%，輕旱情、中旱情總的發生面積接近75.00%;其次為2012 年5 月17 日，該時相輕旱情、中旱情總發生面積接近50.00%，但中旱情發生面積比率僅為3.78%;旱情程度最輕的是2012 年4 月28 日，全區域沒有干旱發生，且濕潤地區超過60.00%;2012 年3 月26 日出現了輕度旱情，但發生面積較小，不足20.00%。

同期降水數據顯示，宿遷市2012 年2 月～5 月降水比歷年同期少70.00%，春旱旱情明顯，雖然期間發生過幾次降水過程，但主要集中在2 個時段內:第一個時段為2 月中旬至3 月中旬出現了低溫連陰雨天氣，這次降水過程雖然持續時間長，但是雨量不大，且前期的旱情較重，所以3 月26 日的遙感影像仍然監測出部分地區有輕度干旱發生;第二個時段為4 月下旬至5 月初，分別為4 月20 日～21 日、4月25 日、4 月29 日～5 月1 日，這次的寡照連陰雨天氣持續時間長且降水量大，對應本試驗4 月28 日反演結果全市絕大部分土壤呈濕潤狀態。之后天氣晴好，農田水分蒸散導致干旱程度逐漸增強。

3 討論

利用HJ 數據提取研究區的NDVI、TS影像，構建NDVI-TS特征空間并確定其干、濕邊參數，分析了TVDI對不同土壤深度含水量的估算能力，主要結論如下:

表4 各等級旱情占全區面積的比率Table 4 The pecentages of areas with different soil moistures to the whole area

構建NDVI-Ts特征空間時發現當NDVI處于0 ～0.2 區間時，發現由于NDVI對植被覆蓋程度低的地類指示性不足，導致地表溫度最大值并沒有如理論模型一般隨NDVI增大而減小，因此，在進行干、濕邊方程擬合時，只考慮NDVI＞0.2 的區間，這樣做既合理舍棄了NDVI指示性低的區間，又有效提高了高植被區的反演精度。

TVDI與各層土壤實測含水量數據都有一定程度的相關性，尤其是與10 ～20 cm 土層相關系數為-0.65 ～-0.80，說明TVDI對于土壤含水量尤其是淺層土壤具有較高的指示意義。

根據TVDI與10 ～20 cm 土層相對濕度的線性關系，得到10 ～20 cm 土層相對濕度分布圖，并進行旱情等級劃分，同期降水數據證明反演結果是合理的。

TVDI耦合了植被信息與地表溫度信息，在干旱監測領域有著較大優勢。本試驗初步探討了利用HJ 數據TVDI法在蘇北地區干旱監測中的適用性。然而干旱不僅與土壤含水量有關，還受土壤性質、植被類型、種植制度等因素影響，為了更準確地描述受旱程度，耦合多種要素建立更具機理性的監測指數是該領域的一個趨勢。此外，本研究中選取的HJ-1B 影像的熱紅外通道空間分辨率為300 m，相較MODIS 數據的1 000 m 已有較大進步，但仍然達不到田塊尺度要求，今后的研究可以嘗試中、高分辨率數據融合來提高反演精度。

［1］ 張元元.應用FY-2 地表蒸散產品監測西南特大干旱［J］.氣象，2011，37(8):999-1005.

［2］ 郭 巍，魏 明.2010 年春西南酷旱的監測與機理分析［J］.自然資源學報，2011，26(9):1628-1636.

［3］ 孫 麗，王 飛，吳 全.干旱遙感監測模型在中國冬小麥區的應用［J］.農業工程學報，2010，26 (1):243-249.

［4］ 李衛國.作物旱澇災情遙感監測進展與思考［J］.江蘇農業學報，2013，29(6):1503-1506.

［5］ 仝兆遠，張萬昌.土壤水分遙感監測的研究進展［J］.水土保持學報，2007，21(24):107-113.

［6］ ALDERFASI A，NIELSEN D C.Use of crop water stress index for monitoring water status and scheduling irrigation in wheat［J］.Agricultural Water Management，2001，47(1):69-75.

［7］ 熊世為，李衛國，景元書.基于HJ-1B 衛星遙感影像的小麥冠層溫度的反演［J］.江蘇農業學報，2012，28(6):1466-1470.

［8］ 李衛國，趙麗花.中高分辨率遙感影像在小麥監測中的比較［J］.江蘇農業學報，2011，27(4):736-739.

［9］ 李衛國，王紀華，李存軍，等.冬小麥花期生理形態指標與衛星遙感光譜特征的相關性分析［J］.麥類作物學報，2009，29(1):79-82.

［10］ 馮海霞，秦其明，蔣洪波，等.基于HJ-1A/1B CCD 數據的干旱監測［J］.農業工程學報，2011，27(增1):358-365.

［11］ WAN Z M，WANG P X，LI X W.Using MODIS land surface temperature and normalized difference vegetation index products for monitoring drought in the southern great plains，USA［J］.International Journal of Remote Sensing，2003，24:1-12.

［12］ 李衛國，趙春江，王紀華，等.基于衛星遙感的冬小麥拔節期長勢監測［J］.麥類作物學報，2007，27(3):523-527.

［13］ 趙少華，秦其明，張 峰，等.基于環境減災小衛星(HJ-1B)的地表溫度單窗反演研究［J］.光譜學與光譜分析，2011，31(6):1152-1156.

［14］ JIMENEZ-MUNOZ J C，SOBRINO J A.A generalized singlechannel method for retrieving land surface temperature from remote sensing data［J］.Journal of Geophysical Research，2003，108:4688-4696.

［15］ 楊沈斌，趙小燕，申雙和.基于Landsat TM/ETM+數據的北京城市熱島季節特征研究［J］.大氣科學學報，2010，33(4):427-435.

［16］ 覃志豪，李文娟，徐 斌，等.陸地衛星TM6 波段范圍內地表比輻射率的 估計［J］.國土資源遙感，2004(3):28-32.

［17］ 熊世為，景元書，李衛國.基于HJ-1B 遙感數據的冬小麥旱情監測研究［J］.麥類作物學報，2013，33(1):84-88.

［18］ 姚春生，張增祥，汪 瀟.使用溫度植被干旱指數(TVDI)法反演新疆土壤濕度［J］.遙感技術與應用，2004，19(6):473-478.