基于NDVI的喀斯特水資源遙感定量監測及分析——以貴州省為例

賀中華，陳曉翔，梁 虹，黃法蘇，趙 芳

（1.中山大學 地理科學與規劃學院，廣州510275；2.貴州師范大學 地理與環境科學學院，貴陽550001；3.貴州省水文水資源局，貴陽550002；4.貴州省貴陽市白云區職業技術學校，貴陽550014）

喀斯特無論在世界上還是我國，都是一類脆弱的生態環境，已引起國內外學術界的深切關注。喀斯特流域是具有特殊的雙重含水介質，特殊的地表、地下雙重分水嶺，獨特的地貌—水系結構的地域綜合體［1］。喀斯特水資源是大氣降水在喀斯特流域下墊面再分配的表現，根據喀斯特流域的特征，其水資源可分為地表水資源和地下水資源。影響喀斯特水資源的因素很多，除土地利用類型、巖組類型、地貌類型外，其植被類型也不容忽視。例如高大的喬木和低矮的灌木對減小降水對地表的沖擊力不同，不同森林植被通過影響降水在地表的側向流速、降水在流域下墊面滯流的時間及下滲率，進而影響喀斯特流域的水資源量。目前，衡量流域植被好壞的一個重要指標為植被指數。植被指數（VI）是對地表植被活動的簡單、有效和經驗的度量，在一定程度上反映流域下墊面的賦水信息。經過近20a的發展，植被指數已有幾十種，其中歸一化植指數（NDVI）被廣泛地應用在全球與區域土地覆蓋、植被分類和環境變化等研究領域［2－9］。對于喀斯特水資源的研究，課題組曾做過相關工作［10－11］，而基于 NDVI的喀斯特水資源的定量研究，無論在國內還是國外，未曾見有研究報道。本文在貴州省內選取20個具有連續5a觀測水文數據和遙感資料的典型喀斯特流域，利用遙感技術，從TM影像中提取喀斯特流域NDVI，利用現代數學方法，探討喀斯特流域水資源與NDVI的關系，建立水資源監測、預測模型，并利用5個研究樣區進行模型檢驗。

1 數據選取及處理

1.1 水文數據選取

根據貴州省水文總站整編的《貴州省歷年各月平均流量統計資料》以及貴州省水文水資源局整編的《貴州省水資源公報》，選其中都處于相同的氣候帶的20個水文斷面，時間從2005—2010年，流域面積以中小流域為主，目的是為了保證流域下墊面的地質條件能盡可能相同或相近。各流域的9月平均徑流深見表3。

1.2 遙感數據選取

數據選用TM影像的2005—2010年，成像時間分別為每年的9月，保證降雨對流域賦水影響較小，保證每個時段研究樣區云量小于30%。

1.3 遙感影像處理及計算

1.3.1 遙感影像預處理

（1）大氣校正。目前，大氣校正的方法有很多，其中，大氣輻射傳輸模型是大氣校正中精度較高的方法。它是利用電磁波在大氣中的輻射傳輸原理建立的模型來對遙感圖像進行大氣校正的。本研究采用的FLAASH模型是改進的MORTRAN模型，它不僅可以對高光譜數據進行大氣校正，而且還可以對多光譜數據如 Landsat，SPOT，AVHRR，MERIS，IRS和ASTER等數據進行大氣校正。

（2）幾何校正。幾何校正包括圖像對地形圖和圖像對圖像的配準。地形圖是國家基礎地理信息1∶25萬數據，其坐標為地理坐標、采用克拉索夫斯基橢球。影像配準利用多項式中4項式，控制點選擇40個左右，配準精度在5個像素內；精度圖像對圖像配準保證在0.3個像素，為了保證光譜信息，重采樣時選取最近鄰法。

1.3.2 表觀反射率計算

（1）光譜輻射亮度的計算

如果沒有定標參數Gain和Bias的資料，某一波段的L可以根據式（2）計算。

式中：QCAL——某一像元的 DN 值，即 QCAL＝DN；QCALmax——像 元 可 以 取 的 最 大 值 255；QCALmin——像元可以取的最小值。對于Landsat—7來說，式（2）［12－13］可改為式（3）（QCALmin＝1）。

（2）表觀反射率的計算［13－15］。

式中：ρ——大氣層頂（TOA）表觀反射率（無量綱）；π——常量（球面度sr）；L——大氣層頂進入衛星傳感器；D——日地之間距離。根據表1，可以推算全年任何一天的日地距離；ESUN——大氣層頂的平均太陽光譜輻照度，根據表2可查得［12］。

θ為太陽的天頂角，地面站提供的頭文件給出的是太陽高度角，因此θ＝90°－β。另外，可以用式（5）直接求取cosθ［16］。

式中：φ——地理緯度；δ——太陽赤緯；h——太陽的時角。

1.3.3 NDVI指數的計算 根據植被的反射光譜特征，利用紅波段、近紅外波段的反射率和其他因子及其組合所獲得的植被指數來提取植被信息，且這些波段常包含90%以上有關植被的信息。歸一化植被指數是廣泛使用的一種植被指數，由Rouse等人提出［17］。

式中：NIR——近紅外通道反射率；R——紅色通道的反射率。其中，－1≤NDVI≤1，負值表示地面覆蓋為云、水、雪等，對可見光高反射；0表示有巖石或裸土等，NIR和R近似相等；正值表示有植被覆蓋，且隨覆蓋度增大而增大。

表1 隨時間變化日地距離（天文單位，D）

表2 Landsat－7和Landsat－5的大氣層頂平均太陽光譜照度 W／（m2·μm）

本研究選用Landsat—7數據，首先，選用公式（2）計算光譜輻射亮度；其次，選用公式（4）、公式（5），并根據表1、表2計算表觀反射率；再次，利用公式（6）計算光譜輻射亮度的LNDVI、表觀反射率的ρNDVI，得表3。

表3 喀斯特流域研究樣區水文數據及NDVI

2 遙感模型建立

2.1 模型原理

假定喀斯特流域某水文斷面觀測值Y和該流域植被指數X之間關系可用如下模型表示［18］：

其中，b0，b1，b2，b3是未知因素參數；ε—N（0，σ2）隨機變量。為評價回歸方程的精度，需對其進行顯著性檢驗，通常用F檢驗。

因F服從自由度為（m，n－m－1）的F分布，對于指定的α，由F分布表可查Fα（m，n－m－1），如F＞Fα（m，n－m－1），則認為回歸模型適合該組資料稱它為顯著的，否則稱為不顯著即不能使用。

2.2 模型建立

首先，根據表3，借助SPSS和Matlab統計軟件，利用公式（8）計算喀斯特流域水資源與NDVI的相關關系，得到表4；其次，利用公式（7）建立喀斯特水資源監測、預測模型，模型系數如表5所示。圖1，圖2表示研究樣區水資源徑流深與NDVI的擬合效果。

表4 徑流深與NDVI相關系數矩陣

（1）從表4可知，喀斯特水資源與其植被指數相關性都很高，尤其是地物光譜輻射亮度的歸一化植被指數，高達0.857；其次，水資源與地物表觀反射率的歸一化植被指數的相關性也很高（0.652）；另外，反映喀斯特植被覆蓋率的NDVI之間的相關性也很好（0.866）。

（2）從表5可知，分別由地物光譜輻射亮度的歸一化植被指數、地物表觀反射率的歸一化植被指數來擬合喀斯特水資源，其擬合的效果很好，如圖1—2所示，擬合度都很高，尤其是地物表觀反射率的歸一化植被指數來擬合喀斯特水資源，其值達0.971；利用公式（9）對擬合的效果進行F檢驗，其F的最大值達176.832，最小值為69.815，均大于給定的臨界值5.29，說明由此建立的模型高度顯著。

表5 模型系數表

（3）根據公式（7），利用表5，其喀斯特流域水資源監測、預測模型可表達為：

圖1 喀斯特流域研究樣區水資源徑流深與LNDVI擬合效果

綜上所述，在喀斯特地區，由于地表崎嶇，地下洞隙縱橫交錯，水文動態變化劇烈，地表水滲漏嚴重，地下持水保水能力差；土層薄、肥力低、植被生長困難，水土流失嚴重，形成了獨特的、脆弱的喀斯特自然環境，嚴重地制約喀斯特流域的持水、供水能力。喀斯特流域具有特殊的雙重含水介質，形成獨特的地表－地下水系結構，因此，喀斯特流域與正常流域特別是濕潤地區常態流域相比，其流域水資源的形成機制、空間分布規律具有一定的特殊性。流域植被類型及覆蓋率將直接影響降雨在喀斯特流域入滲及徑流，即影響降雨在流域空間的再分配，因此，喀斯特流域植被指數是喀斯特流域賦水狀況的重要性指標。

圖2 喀斯特流域研究樣區水資源徑流深與ρNDVI擬合效果

3 遙感模型檢驗

為了評定監測、預測模型的精度，任選5個喀斯特流域作為研究樣區，按上述方法對研究樣區遙感影像進行處理，分別提取LNDVI、ρNDVI，如表6所示，分別代入模型（10），（11）進行計算，并與實測數值對比見表6。通過計算比較得出，模型（10），（11）相對誤差值都比較小，說明用這兩個模型對喀斯特流域水資源進行監測、預測，效果是比較理想的，且模型（11），效果更好，精度更高。從理論上分析，原始遙感影像的DN是未經過任何校正，包括輻射定標校正，只是進入傳感器中的輻射能的一種數字轉換形式，不能本質地反映地物的輻射特性。L和ρ都經過了輻射定標校正，但是，當ρ再經過大氣校正后，它就是地物的反射率，能本質地反映地物的輻射特性，因此，由它構成的NDVI植被指數最接近地物的NDVI。

表6 模型檢驗結果

4 結論與分析

（1）喀斯特具有特殊的流域下墊面介質結構，其流域產、匯機制復雜，流域賦水影響因素多樣，流域植被覆蓋率起到重要的作用。

（2）利用地物光譜輻射亮度的歸一化植被植指數（LNDVI）和地物表觀反射率的歸一化植被指數（ρNDVI）對喀斯特流域水資源進行監測、預測效果很好，尤其是利用ρNDVI進行監測、預測，精度更高。

（3）適合于喀斯特流域水資源監測、預測數學模型是：

通過方差分析和樣區檢驗，得出很好的預測效果。

［1］ 楊明德，譚明，梁虹.喀斯特流域水文地貌系統［M］.北京：地質出版社，1998.

［2］ 何隆華，儲開華，肖向明.Vegetation圖像植被指數與實測水稻葉面積指數的關系［J］.遙感學報，2004，8（6）：672－676.

［3］ 王福民，黃敬峰，唐延林，等.采用不同光譜波段寬度的歸一化植被指數估算水稻葉面積指數［J］.應用生態學報，2007，18（11）：2444－2450.

［4］ 楊曦，武建軍，閆峰，等.基于地表溫度—植被指數特征空間的區域土壤干濕狀況［J］.生態學報，2009，29（3）：1205－1216.

［5］ 李玉霞，楊武年，童玲，等.基于光譜指數法的植被含水量遙感定量監測及分析［J］.光學學報，2009，29（5）：1403－1407.

［6］ Fabio Maselli1，Antonio Di Gregorio，Valerio Capecchi1 et al.Enrichment of land－cover polygons with eco－climatic information derived from MODIS NDVI imagery［J］.Journal of Biogeography，2009，36：639－650.

［7］ Karnieli A，Nurit A，Rachel T，et al.Use of NDVI and land surface temperature for drought assessment［J］.Journal of Climate，2010，23（3）：618－633.

［8］ Pu Ruiliang，Gong Peng，Tian Yong，et al.Using classification and NDVI differencing methods for monitoring sparse vegetation coverage：a case study of saltcedar in Nevada，USA［J］.International Journal of Remote Sensing，2008，29（14）：3987－4011.

［9］Thenkabail P S.Inter－sensor relationships between IKO－NOS and Landsat－7ETM＋NDVI data in three ecoregions of Africa［J］.International Journal of Remote Sensing，2004，25（2）：389－408.

［10］ 賀中華，梁虹，黃法蘇，等.基于RS流域枯水資源的判讀識別［J］.貴州師范大學學報：自然科學版，2004，22（2）：36－39.

［11］ 賀中華，楊勝天，梁虹，等.基于GIS和RS的喀斯特流域枯水資源影響因素識別：以貴州省為例［J］.中國巖溶，2004，23（1）：48－55.

［12］ Chander G，Markham B L.Revised Landsat－5TM radiometric calibration procedures and post－calibration dynamic ranges［J］.Geoscience and Remote Sensing，2003，41（1）：2674－2677.

［13］ Markham B L，Barker J L.Thematic Mapper bandpass solar exoatmospheric irradianees［J］.International Journal of Remote Sensing，1987，8（3），513－523.

［14］ Roderick M，Smith R，Lodwick G.Calibrating longterm AVHRR－derived NDVI imagery［J］.Remote Sensing of Environment，1996，58（1）：1－12.

［15］ Huete A.Overview of the radiometric and biophysical performance of the MODIS vegetation indices［J］.Remote Sensing of Environment，2002，83（1）：195－213.

［16］ LüS H.Physical Foundation of Remote Sensing［M］.Beijing：Commercial Press，1981：206－234.

［17］ Rouse J W，Haas R H，Schell J A，et al.Monitoring vegetation systems in the great plains with ERTS［J］.Third ERTS Symposium，NASA，1973.

［18］ 張超，楊秉賡.計量地理學基礎［M］.北京：高等教育出版社，1985.