順布日， 青 松，2， 包玉海，2， 趙文靜

(1.內(nèi)蒙古師范大學(xué) 地理科學(xué)學(xué)院， 內(nèi)蒙古 呼和浩特010022；2.內(nèi)蒙古自治區(qū)遙感與地理信息系統(tǒng)重點(diǎn)實(shí)驗(yàn)室， 內(nèi)蒙古 呼和浩特 010022；3.生態(tài)環(huán)境部華南環(huán)境科學(xué)研究所近岸海域環(huán)境研究中心， 廣東 廣州 510535)

地表水是地球上重要的淡水資源[1]。據(jù)上世紀(jì)80年代統(tǒng)計(jì)內(nèi)蒙古自治區(qū)境內(nèi)有1 500余個(gè)湖泊及千余條河流[2]。近年來(lái)，入湖徑流量補(bǔ)給不足，加上氣候引起的降水量減少、氣溫升高，人為活動(dòng)引起的礦產(chǎn)面積和灌溉面積的增加，導(dǎo)致湖泊面積與數(shù)量急劇減少[3-5]。幾百年以來(lái)，湖泊水資源一直維持著內(nèi)蒙古草原區(qū)湖泊及其周邊的濕地生態(tài)系統(tǒng)，發(fā)揮著重要的作用[6-7]?？焖俣_地監(jiān)測(cè)湖泊面積變化，對(duì)水資源監(jiān)測(cè)和管理具有重要意義[8]。

近年來(lái)遙感技術(shù)的飛速發(fā)展和普及，越來(lái)越多的國(guó)內(nèi)外研究人員通過(guò)遙感影像監(jiān)測(cè)水體。在目前的光學(xué)遙感中，地表水主要分類方法如下：?jiǎn)尾ǘ伍撝捣?、多波段譜間關(guān)系法、水體指數(shù)法等。單波段閾值法：通過(guò)多光譜遙感影像中對(duì)水體高吸收的波段，將水體與其他地物區(qū)分，最后利用最佳閾值將水體分離出來(lái)[9]。陸家駒等[10]基于Landsat TM影像的近紅外波段在清流河流域提取了地表水。結(jié)果發(fā)現(xiàn)單波段閾值法只能提取面積大于4 000 m2的大型水體，無(wú)法區(qū)分水體和山區(qū)的陰影。多波段譜間關(guān)系法：利用遙感影像的兩個(gè)或兩個(gè)以上的波段參與運(yùn)算，分析各個(gè)波段之間存在的關(guān)系、水體與其他地物的光譜特征，從而建立反映地物特征的邏輯關(guān)系式，提取水體[11]。汪金花等[11]基于Landsat TM影像利用譜間關(guān)系法提取了唐山市陡河水庫(kù)中部水體。結(jié)果表明, 該算法精度高于一般的閾值分類法。水體指數(shù)法：通過(guò)比值運(yùn)算，對(duì)遙感影像中對(duì)水體反射率低的波段作為分母，對(duì)水體反射率高的波段作為分子，從而提高水體和其周?chē)匚锊顒e，建立關(guān)系式對(duì)水體進(jìn)行提取的方法。Mcfeeters等[12]利用Landsat MSS影像的綠色波段和近紅外波段進(jìn)行比值運(yùn)算，建立了歸一化差異水體指數(shù)(normalized difference water index, NDWI)，該指數(shù)可以較好地消除山體產(chǎn)生的陰影。徐涵秋[13]利用ETM+遙感影像的綠光波段和中紅外波段進(jìn)行比值運(yùn)算，提出了修正的歸一化差異水體指數(shù)(modified normalized difference water index, MNDWI)，該指數(shù)能很好地區(qū)分城鎮(zhèn)地區(qū)中的水體，且水體提取精度高。Sun等[14]基于HJ-1A/B衛(wèi)星利用面向?qū)ο蟮姆诸惙椒▽?duì)上海市部分區(qū)域的水體進(jìn)行了提取，該方法分類結(jié)果優(yōu)于單一的光譜分類方法。Lu等[15]利用HJ-1A/B衛(wèi)星影像的近紅外波段和NDVI，NDWI指數(shù)以及地形坡度數(shù)據(jù)相結(jié)合，建立了水體綜合提取方法，并應(yīng)用到3個(gè)復(fù)雜的水域當(dāng)中。Liu等[16]在中國(guó)與亞洲戈壁沙漠接壤的半干旱地區(qū)，基于Landsat衛(wèi)星監(jiān)測(cè)了從1975—2009年期間湖泊面積變化。研究發(fā)現(xiàn)241個(gè)湖泊中有121個(gè)完全干涸。

基于水體指數(shù)法的提取精度和遙感影像的波段數(shù)量的增加，加上全球?qū)λY源研究的需求不斷增加，大尺度范圍提取水體成為熱點(diǎn)。Tao等[4]將每3 a或5 a Landsat遙感影像合為一期影像，研究了蒙古高原1976年至2010年湖泊數(shù)量和面積變化。結(jié)果表明，研究期間整個(gè)蒙古高原湖泊面積急劇減少，其中內(nèi)蒙古地區(qū)最為明顯。Tao等進(jìn)一步定量分析了自然因素和人類活動(dòng)的影響之后發(fā)現(xiàn)，在內(nèi)蒙古的草原地區(qū)，煤炭的開(kāi)采直接導(dǎo)致湖泊資源減少，而在農(nóng)業(yè)生產(chǎn)區(qū)，灌溉成為湖泊資源減少的主要因素。Zhang等[5]利用Landsat遙感數(shù)據(jù)，1970—2013年期間的蒙古高原和青藏高原湖泊進(jìn)行了對(duì)比研究。研究表明，這兩個(gè)高原湖泊面積和數(shù)量向相反的方向變化，即蒙古高原由于溫度升高與降水量的減少導(dǎo)致湖泊面積急劇縮小并且有208個(gè)湖泊消失。相反青藏高原由于氣溫升高導(dǎo)致冰凍圈融化及降水量的增加，新增了99個(gè)湖泊。Zhou Y等[3]基于Google Earth Engine平臺(tái)，利用 Landsat遙感影像繪制了蒙古高原1991—2017年期間湖泊面積變化和數(shù)量變化圖。通過(guò)分析發(fā)現(xiàn)內(nèi)蒙古的湖泊變化比蒙古國(guó)更為劇烈，在2009年之前，自然和人為因素對(duì)湖泊動(dòng)態(tài)都有顯著影響，但在2009年之后，降水量對(duì)湖泊的恢復(fù)起著越來(lái)越重要的作用。

HJ-1 A/B遙感影像具有高時(shí)間分辨率和大范圍監(jiān)測(cè)特點(diǎn)，是對(duì)大尺度動(dòng)態(tài)監(jiān)測(cè)的理想數(shù)據(jù)源[14,17-18]。因此，本文以HJ-1A/B遙感影像作為數(shù)據(jù)源，對(duì)內(nèi)蒙古自治區(qū)地表水面積進(jìn)行動(dòng)態(tài)監(jiān)測(cè)，以期為內(nèi)蒙古地區(qū)生態(tài)環(huán)境保護(hù)以及修復(fù)研究和水資源管理提供依據(jù)。

1 研究區(qū)概況

內(nèi)蒙古自治區(qū)屬于蒙古高原的一部分，位于中國(guó)北部邊疆。有9個(gè)地級(jí)市和3個(gè)盟，共12個(gè)盟(市)。境內(nèi)有嫩江、西遼河、黃河和額爾古納河四大水系，大小河流千余條。有近千個(gè)湖泊，如岱海、達(dá)里諾爾湖、呼日查干淖爾、呼倫湖、貝爾湖等[2,19]。大多數(shù)湖泊屬內(nèi)陸湖泊，近年來(lái)湖泊數(shù)量與面積急劇減少[3-5]。圖1為研究區(qū)概況。在研究區(qū)范圍內(nèi)均勻地選取了4景HJ-1A/B影像，把1景HJ-1A/B影像作為一個(gè)實(shí)驗(yàn)區(qū)，一共A，B，C，D 4個(gè)實(shí)驗(yàn)區(qū)，用于選取地物純像元，得到光譜箱型圖。4個(gè)實(shí)驗(yàn)區(qū)包含4個(gè)GF-2 PMS影像，用于HJ-1A/B影像選取純像元時(shí)作參考以及驗(yàn)證地表水提取算法精度。其詳細(xì)信息如表1所示。所選的區(qū)域具有不同的地表水類型，包括水生環(huán)境復(fù)雜程度各異的湖泊、不規(guī)則的線狀河流、人工水渠通道等。同時(shí)充分考慮了影響地表水提取精度的影響因素。

圖1 研究區(qū)及4個(gè)實(shí)驗(yàn)區(qū)地形圖

表1 4個(gè)典型子區(qū)的詳細(xì)描述

2 數(shù)據(jù)與方法

2.1 遙感數(shù)據(jù)

研究使用的數(shù)據(jù)為從2009—2018年，共217景無(wú)云覆蓋或云量少、數(shù)據(jù)質(zhì)量好的HJ-1A/B 影像(從中國(guó)資源衛(wèi)星應(yīng)用中心獲取，http:∥www.cresda.com/CN/)。由于遙感影像的成像時(shí)間會(huì)直接影響地表水面積，因此，本文挑選的影像成像時(shí)間在6到9月汛期(豐水期)。HJ-1A/B遙感數(shù)據(jù)預(yù)處理方法步驟如下所述。

首先，利用公式(1)將遙感影像DN值轉(zhuǎn)為輻射亮度值。

(1)

式中：LTOA為輻射亮度值(W/m2·μm·sr)；a和L0分別為HJ-1A/B數(shù)據(jù)絕對(duì)定標(biāo)系數(shù)增益和絕對(duì)定標(biāo)系數(shù)偏移量。

其次，利用公式(2)獲取大氣層頂反射率：

(2)

式中：ρTOA為大氣層頂反射率(sr)；d為日地距離；F0為波段平均的大氣層外太陽(yáng)輻照度(W/m2·μm)；θ0為太陽(yáng)高度角。

此外，本文使用了4景GF-2 PMS影像(從中國(guó)資源衛(wèi)星應(yīng)用中心獲取，http:∥www.cresda.com/CN/)用于選取純像元及HJ-1A/B分類結(jié)果精度驗(yàn)證。成像時(shí)間準(zhǔn)同步于4個(gè)實(shí)驗(yàn)區(qū)。

2.2 氣象數(shù)據(jù)

氣象數(shù)據(jù)由國(guó)家氣象科學(xué)數(shù)據(jù)中心(http:∥data.cma.cn/)提供。在內(nèi)蒙古范圍內(nèi)共選115個(gè)氣象站點(diǎn)。通過(guò)年降水量和年月均溫度來(lái)研究?jī)?nèi)蒙古2009—2018年間氣象變化。

2.3 人為因素?cái)?shù)據(jù)

在本次研究中，將煤炭產(chǎn)量、有效灌溉面積作為人類活動(dòng)的兩個(gè)指標(biāo)，分析其在內(nèi)蒙古地表水變化過(guò)程中的影響。其中，內(nèi)蒙古自治區(qū)的煤炭產(chǎn)量數(shù)據(jù)(2009—2016年)、農(nóng)田灌溉數(shù)據(jù)(2009—2017年)來(lái)源于內(nèi)蒙古統(tǒng)計(jì)年鑒。

2.4 方 法

2.4.1 光譜建立及閾值確定 本文利用HJ-1A/B影像的近紅外波段對(duì)水體高吸收的特征，采用NDWI(normalized difference water index)水體指數(shù)對(duì)地表水進(jìn)行提取。計(jì)算公式如下：

(3)

式中：ρGreen和ρNIR分別為HJ-1A/B影像的綠色波段和紅外波段。NDWI可以突出水體信息并且有效地減少植被的影響[11]。

圖2a為NDWI光譜指數(shù)箱型圖。在每個(gè)實(shí)驗(yàn)區(qū)內(nèi)選取水體、云、陰影、建筑物、耕地、裸地、露天礦等7類地物的樣點(diǎn)，每個(gè)地物分別取500個(gè)純像元(樣點(diǎn))。通過(guò)計(jì)算7種類型地物NDWI光譜指數(shù)的統(tǒng)計(jì)量，以確定其取值范圍，從而確定初始閾值。初始閾值為四個(gè)實(shí)驗(yàn)區(qū)閾值的平均值(T=0.14)。如果在初始閾值下提取地表水誤差較大，那么對(duì)初始閾值進(jìn)行調(diào)整，再確定最佳閾值。在最佳閾值下仍然有誤差，是因?yàn)榘迅蓴_因素(如露天礦)錯(cuò)分成水體，那么需要對(duì)干擾因素進(jìn)行掩膜。圖2b為第二波段的箱型圖，通過(guò)HJ-1A/B影像的第二波段可以對(duì)云進(jìn)行剔除，閾值為4個(gè)實(shí)驗(yàn)區(qū)的平均值(T≈0.36)。可以利用簡(jiǎn)單的決策樹(shù)方法將云剔除，提取地表水。對(duì)4景GF-2 PMS影像進(jìn)行目視解譯，提取地表水，作為HJ-1A/B影像分類結(jié)果的驗(yàn)證依據(jù)。

圖2 NDWI及第二波段光譜指數(shù)箱型圖

2.4.2 精度驗(yàn)證 圖3為GF-2 PMS目視解譯結(jié)果與HJ-1A/B自動(dòng)提取結(jié)果的地表水面積散點(diǎn)圖。

從圖3a可以看出，散點(diǎn)在1∶1線周?chē)休^好的分布，湖泊面積大于0.15 km2時(shí)，多數(shù)散點(diǎn)在1∶1線下方，湖泊面積小于0.15 km2時(shí)，較多散點(diǎn)在1∶1線上方。這是由GF-2 PMS高分辨率數(shù)據(jù)能提取小面積湖泊和細(xì)小的河流，然而HJ-1A/B數(shù)據(jù)對(duì)細(xì)小的河流和小面積湖泊的提取能力較差，因此把混合像元錯(cuò)分為水體所造成，精度R2為0.96；從圖3b可知，散點(diǎn)分布于1∶1線上方且較少，這是因?yàn)辄S旗海為季節(jié)性內(nèi)陸湖，常年干涸或積水較少，影像的成像日期恰好處于湖泊干涸狀態(tài)，因此提取的水體較少，R2為0.82。而圖3c和圖3d顯示，散點(diǎn)整體分布于1∶1線上方。這是HJ-1A/B數(shù)據(jù)提取細(xì)小的河流和湖泊淺水區(qū)域面積較少，GF-2 PMS影像提取的水體面積較多導(dǎo)致，R2分別為0.98和0.97。

圖3 實(shí)驗(yàn)區(qū)地表水面積目視解澤結(jié)果和自動(dòng)提取結(jié)果精度對(duì)比散點(diǎn)圖

由此可知，本文所選取NDWI水體指數(shù)方法可用于內(nèi)蒙古地表水的提取。

3 結(jié)果與分析

3.1 內(nèi)蒙古2009-2018年地表水面積時(shí)空分布特征

圖4a顯示了2009—2018年間內(nèi)蒙古地表水面積變化趨勢(shì)。2009—2018年期間內(nèi)蒙古地表水面積為增長(zhǎng)趨勢(shì)，增長(zhǎng)的速率為72 km2/a。2009年的面積為3 682.19 km2，2018年的面積為4 893.20 km2，增加了1 211.01 km2；2009—2013年期間內(nèi)蒙古地表水面積以每年247 km2速率持續(xù)增長(zhǎng)，增加了913.23 km2；2013—2017年期間地表水面積以每年150 km2減少，面積減少了752.12 km2；2018年地表水面積增長(zhǎng)到最高點(diǎn)，為4 872.12 km2。圖4b為內(nèi)蒙古各盟(市)地表水面積變化圖。研究發(fā)現(xiàn)2009—2018年期間各盟(市)地表水變化有不同的趨勢(shì)。通過(guò)MK趨勢(shì)檢驗(yàn)和回歸分析，4個(gè)盟(市)地表水面積顯著增加，即阿拉善盟、烏海市、巴彥淖爾市和呼倫貝爾市，增加的面積分別為129.96，75.83，92.09，350.48 km2。對(duì)于內(nèi)蒙古地表水面積來(lái)說(shuō)，在各盟(市)中呼倫貝爾市地表水占據(jù)比重最大，占內(nèi)蒙古地表水總面積的50%以上。其次為鄂爾多斯市，占11%左右。烏海市占最少，約1%。

圖4 內(nèi)蒙古全區(qū)和分區(qū)地表水面積變化特征

為了便于研究影響內(nèi)蒙古各盟(市)地表水面積變化的驅(qū)動(dòng)力，我們以影響因素的相似性將各盟(市)分為兩大部分，即：自然因素影響地表水變化的盟(市)和人類活動(dòng)影響地表水變化的盟(市)。

3.2 地表水變化影響因素分析

研究區(qū)自然因素選取了月平均溫度和月距平降水兩個(gè)指標(biāo)，將農(nóng)田有效灌溉面積和煤炭產(chǎn)量作為人類活動(dòng)兩個(gè)指標(biāo)。另外還考慮了徑流量和土地利用變化。通過(guò)一元線性回歸模型來(lái)分析內(nèi)蒙古地表水面積變化與當(dāng)?shù)貧庀笞兓腿祟惢顒?dòng)的相關(guān)性(詳見(jiàn)表2)。年降水量和河流的徑流量?jī)蓚€(gè)自然因素是內(nèi)蒙古部分盟(市)地表水變化的主要原因之一。建造水庫(kù)，填埋、截流河流進(jìn)行耕地是主要的人為因素。

表2 內(nèi)蒙古各盟(市)驅(qū)動(dòng)力(P值)分析

由表2及圖5—7可知，受年降水量影響的盟(市)鄂爾多斯市地表水面積在研究期間有輕微增加趨勢(shì)。通過(guò)回歸分析得到年降水量對(duì)地表水變化顯著(p<0.05)。包頭市經(jīng)過(guò)2009，2010，2011年3年的干旱，從2012年起降水量開(kāi)始逐步增加，從而包頭市南部和達(dá)茂旗境內(nèi)湖泊群及艾不蓋河面積增加[22-23]，降水量對(duì)地表水顯著(p<0.05)。赤峰市地表水面積2009—2013年期間增加，2013—2018年期間快速減少。紅山水庫(kù)和達(dá)里諾爾面積變化對(duì)赤峰市地表水面積變化顯著。而降水量影響著紅山水庫(kù)和達(dá)里諾爾面積。受徑流量變化變化影響的盟(市)：興安盟地表水面積在2009—2013年期間持續(xù)增加，2014年開(kāi)始面積減少。綽爾河和洮兒河徑流量變化對(duì)興安盟地表水變化顯著。

圖5 內(nèi)蒙古主要湖泊呼倫湖2009-2018年間面積變化

然而部分盟(市)受年降水量和河流徑流量?jī)蓚€(gè)自然因素的影響。如阿拉善盟地表水面積逐年增加主要由年降水量和居延海面積增加所導(dǎo)致(p<0.05)。居延海面積是由2009—2018年期間黑河徑流量連續(xù)增加而增加[20]。巴彥淖爾市地表水面積在2009—2013年期間緩慢增加，在2014—2018年期間急劇增加，兩個(gè)階段的地表水面積增加跟年降水量和黃河徑流量有關(guān)[21]。呼倫貝爾市地表水面積在2009—2014年期間持續(xù)增長(zhǎng)，在2014—2018年期間緩慢減少。地表水面積主要由呼倫湖(圖5)面積變化導(dǎo)致。呼倫湖面積變化是由入湖徑流的克魯倫河和烏爾遜河以及地下水補(bǔ)給和年降雨量有關(guān)[24]。通過(guò)回歸分析發(fā)現(xiàn)年均溫度對(duì)年均溫度顯著(p<0.05)。這結(jié)果可能隨著氣溫的上升呼倫貝爾地區(qū)地表水蒸發(fā)量增加，從而導(dǎo)致面積減少。

建造水庫(kù)導(dǎo)致地表水發(fā)生劇烈變化的盟(市)：烏海市地表水面積在2009—2018年期間持續(xù)增加，2014年地表水面積急劇增加是由烏海市建造烏海湖引入黃河水導(dǎo)致[25]。錫林郭勒盟地表水面積在2009—2013年期間增加，2013—2018年期間快速減少。主要原因是在東查干淖爾東湖(圖6)上養(yǎng)殖漁業(yè)，將西湖的水源阻斷，導(dǎo)致西湖干涸。在烏拉蓋湖上游建立水庫(kù)導(dǎo)致湖面干涸。

圖6 內(nèi)蒙古主要湖泊呼日查干淖爾2009-2018年間面積變化

部分盟(市)因填埋、截流河流進(jìn)行耕地使地表水面積發(fā)生變化。如呼和浩特市地表水面積變化趨勢(shì)復(fù)雜，即2009—2011年期間地表水面積減少，到2012面積急劇增加，在2012—2015期間面積又開(kāi)始減少，隨后再增加。通過(guò)遙感影像對(duì)比發(fā)現(xiàn)在部分水庫(kù)和河流里填土進(jìn)行耕地減少了地表水面積，如和林格爾縣石咀子水庫(kù)。反而對(duì)大黑河進(jìn)行綜合整治規(guī)劃管理增加了地表水面積[26]。烏蘭察布市地表水在2009—2012年期間增加，2012—2017年期間減少，到2018年開(kāi)始回升。地表水面積減少跟把流入岱海(圖7)和黃旗海的徑流截流耕地有關(guān)[27-30]。

圖7 內(nèi)蒙古主要湖泊岱海2009-2018年間面積變化

通遼市地表水面積在2009—2011年期間連續(xù)減少，2012年出現(xiàn)增加的趨勢(shì)，到2013年地表水面積急劇增加，2013—2018年期間減少。把原有的水庫(kù)或自然湖泊填土進(jìn)行耕作是地表水面積減少的主要原因，如莫力廟和水庫(kù)、他拉干水庫(kù)。已有的研究發(fā)現(xiàn)當(dāng)水資源已成為超采狀態(tài)[30]。人為強(qiáng)烈活動(dòng)是內(nèi)蒙古部分盟(市)地表水面積發(fā)生變化的主要原因。

4 結(jié) 論

(1) 根據(jù)HJ-1A/B遙感數(shù)據(jù)，利用NDWI水體指數(shù)能很好地提取地表水，精度較高。

(2) 內(nèi)蒙古地表水面積在2009-2018年期間整體呈現(xiàn)出增加趨勢(shì)(slope=72 km2/a)；2009—2013年階段地表水面積連續(xù)增加(slope=247 km2/a)； 2013—2017年期間地表水面積急劇減少(slope=150 km2/a)；到了2018年地表水面積回升到最高點(diǎn)。

(3) 2009—2018年10 a期間內(nèi)蒙古各盟(市)地表水面積變化趨勢(shì)有所不同，整體增加的盟(市)有阿拉善盟、巴彥淖爾市和烏海市；以2012年為拐點(diǎn)，地表水面積在2009—2012年期間增加，2012—2017年期間減少，到2018年回升的盟(市)分別為包頭市和烏蘭察布市；以2013年為拐點(diǎn)，地表水面積在2009—2013年期間增加，2013—2018年期間減少的盟(市)有錫林郭勒盟、赤峰市和興安盟。

(4) 本文將影響地表水面積變化的驅(qū)動(dòng)力分為兩大部分，即：自然因素影響地表水變化的盟(市)和人類活動(dòng)影響地表水變化的盟(市)。主要有年降水量和河流的徑流量?jī)蓚€(gè)自然因素影響的內(nèi)蒙古部分盟(市)，如阿拉善盟、巴彥淖爾市、鄂爾多斯市、包頭市、赤峰市、興安盟及呼倫貝爾市；建造水庫(kù)，填埋、截流河流進(jìn)行耕地是主要的人為因素。如烏海市、呼和浩特市、烏蘭察布市、錫林郭勒盟及通遼市。