基于NDVI的喀斯特水資源遙感定量監(jiān)測(cè)及分析——以貴州省為例

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

（1.中山大學(xué) 地理科學(xué)與規(guī)劃學(xué)院，廣州510275；2.貴州師范大學(xué) 地理與環(huán)境科學(xué)學(xué)院，貴陽(yáng)550001；3.貴州省水文水資源局，貴陽(yáng)550002；4.貴州省貴陽(yáng)市白云區(qū)職業(yè)技術(shù)學(xué)校，貴陽(yáng)550014）

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

1 數(shù)據(jù)選取及處理

1.1 水文數(shù)據(jù)選取

根據(jù)貴州省水文總站整編的《貴州省歷年各月平均流量統(tǒng)計(jì)資料》以及貴州省水文水資源局整編的《貴州省水資源公報(bào)》，選其中都處于相同的氣候帶的20個(gè)水文斷面，時(shí)間從2005—2010年，流域面積以中小流域?yàn)橹?，目的是為了保證流域下墊面的地質(zhì)條件能盡可能相同或相近。各流域的9月平均徑流深見(jiàn)表3。

1.2 遙感數(shù)據(jù)選取

數(shù)據(jù)選用TM影像的2005—2010年，成像時(shí)間分別為每年的9月，保證降雨對(duì)流域賦水影響較小，保證每個(gè)時(shí)段研究樣區(qū)云量小于30%。

1.3 遙感影像處理及計(jì)算

1.3.1 遙感影像預(yù)處理

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

（2）幾何校正。幾何校正包括圖像對(duì)地形圖和圖像對(duì)圖像的配準(zhǔn)。地形圖是國(guó)家基礎(chǔ)地理信息1∶25萬(wàn)數(shù)據(jù)，其坐標(biāo)為地理坐標(biāo)、采用克拉索夫斯基橢球。影像配準(zhǔn)利用多項(xiàng)式中4項(xiàng)式，控制點(diǎn)選擇40個(gè)左右，配準(zhǔn)精度在5個(gè)像素內(nèi)；精度圖像對(duì)圖像配準(zhǔn)保證在0.3個(gè)像素，為了保證光譜信息，重采樣時(shí)選取最近鄰法。

1.3.2 表觀反射率計(jì)算

（1）光譜輻射亮度的計(jì)算

如果沒(méi)有定標(biāo)參數(shù)Gain和Bias的資料，某一波段的L可以根據(jù)式（2）計(jì)算。

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

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

式中：ρ——大氣層頂（TOA）表觀反射率（無(wú)量綱）；π——常量（球面度sr）；L——大氣層頂進(jìn)入衛(wèi)星傳感器；D——日地之間距離。根據(jù)表1，可以推算全年任何一天的日地距離；ESUN——大氣層頂?shù)钠骄?yáng)光譜輻照度，根據(jù)表2可查得［12］。

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

式中：φ——地理緯度；δ——太陽(yáng)赤緯；h——太陽(yáng)的時(shí)角。

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

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

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

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

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

表3 喀斯特流域研究樣區(qū)水文數(shù)據(jù)及NDVI

2 遙感模型建立

2.1 模型原理

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

其中，b0，b1，b2，b3是未知因素參數(shù)；ε—N（0，σ2）隨機(jī)變量。為評(píng)價(jià)回歸方程的精度，需對(duì)其進(jìn)行顯著性檢驗(yàn)，通常用F檢驗(yàn)。

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

2.2 模型建立

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

表4 徑流深與NDVI相關(guān)系數(shù)矩陣

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

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

表5 模型系數(shù)表

（3）根據(jù)公式（7），利用表5，其喀斯特流域水資源監(jiān)測(cè)、預(yù)測(cè)模型可表達(dá)為：

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

綜上所述，在喀斯特地區(qū)，由于地表崎嶇，地下洞隙縱橫交錯(cuò)，水文動(dòng)態(tài)變化劇烈，地表水滲漏嚴(yán)重，地下持水保水能力差；土層薄、肥力低、植被生長(zhǎng)困難，水土流失嚴(yán)重，形成了獨(dú)特的、脆弱的喀斯特自然環(huán)境，嚴(yán)重地制約喀斯特流域的持水、供水能力?？λ固亓饔蚓哂刑厥獾碾p重含水介質(zhì)，形成獨(dú)特的地表－地下水系結(jié)構(gòu)，因此，喀斯特流域與正常流域特別是濕潤(rùn)地區(qū)常態(tài)流域相比，其流域水資源的形成機(jī)制、空間分布規(guī)律具有一定的特殊性。流域植被類(lèi)型及覆蓋率將直接影響降雨在喀斯特流域入滲及徑流，即影響降雨在流域空間的再分配，因此，喀斯特流域植被指數(shù)是喀斯特流域賦水狀況的重要性指標(biāo)。

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

3 遙感模型檢驗(yàn)

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

表6 模型檢驗(yàn)結(jié)果

4 結(jié)論與分析

（1）喀斯特具有特殊的流域下墊面介質(zhì)結(jié)構(gòu)，其流域產(chǎn)、匯機(jī)制復(fù)雜，流域賦水影響因素多樣，流域植被覆蓋率起到重要的作用。

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

（3）適合于喀斯特流域水資源監(jiān)測(cè)、預(yù)測(cè)數(shù)學(xué)模型是：

通過(guò)方差分析和樣區(qū)檢驗(yàn)，得出很好的預(yù)測(cè)效果。

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

［2］ 何隆華，儲(chǔ)開(kāi)華，肖向明.Vegetation圖像植被指數(shù)與實(shí)測(cè)水稻葉面積指數(shù)的關(guān)系［J］.遙感學(xué)報(bào)，2004，8（6）：672－676.

［3］ 王福民，黃敬峰，唐延林，等.采用不同光譜波段寬度的歸一化植被指數(shù)估算水稻葉面積指數(shù)［J］.應(yīng)用生態(tài)學(xué)報(bào)，2007，18（11）：2444－2450.

［4］ 楊曦，武建軍，閆峰，等.基于地表溫度—植被指數(shù)特征空間的區(qū)域土壤干濕狀況［J］.生態(tài)學(xué)報(bào)，2009，29（3）：1205－1216.

［5］ 李玉霞，楊武年，童玲，等.基于光譜指數(shù)法的植被含水量遙感定量監(jiān)測(cè)及分析［J］.光學(xué)學(xué)報(bào)，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流域枯水資源的判讀識(shí)別［J］.貴州師范大學(xué)學(xué)報(bào)：自然科學(xué)版，2004，22（2）：36－39.

［11］ 賀中華，楊勝天，梁虹，等.基于GIS和RS的喀斯特流域枯水資源影響因素識(shí)別：以貴州省為例［J］.中國(guó)巖溶，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］ 張超，楊秉賡.計(jì)量地理學(xué)基礎(chǔ)［M］.北京：高等教育出版社，1985.