洪 宸，王猛猛，朱芙瑤

（1. 中國(guó)地質(zhì)大學(xué)（武漢）李四光學(xué)院，湖北 武漢 430074；2. 中國(guó)地質(zhì)大學(xué)（武漢）地理與信息工程學(xué)院，湖北 武漢 430074；3. 中國(guó)地質(zhì)大學(xué)（武漢） 地球物理與空間信息學(xué)院，湖北 武漢 430074）

青藏高原位于我國(guó)西部和西南部，是世界上最高的高原，是我國(guó)甚至亞洲水資源產(chǎn)生、賦存和運(yùn)移的戰(zhàn)略要地。大量研究表明，近年來(lái)青藏高原正在經(jīng)歷著氣候暖濕化過(guò)程，導(dǎo)致多年冰川積雪快速退縮、水資源縮減等現(xiàn)象的發(fā)生[1]。因此，研究青藏高原地區(qū)水體的時(shí)空變化特征并分析其影響因素具有重要意義，可為區(qū)域生態(tài)環(huán)境安全、社會(huì)經(jīng)濟(jì)發(fā)展、管理與規(guī)劃提供科學(xué)依據(jù)[2]。

近年來(lái)遙感技術(shù)成為大尺度研究的主要手段之一，為地表水體的動(dòng)態(tài)監(jiān)測(cè)提供了強(qiáng)有力的支持。Landsat衛(wèi)星影像數(shù)據(jù)因其高空間分辨率，在長(zhǎng)時(shí)序水體監(jiān)測(cè)中得到了較好應(yīng)用[3]。國(guó)內(nèi)外學(xué)者基于遙感數(shù)據(jù)開(kāi)展了水體動(dòng)態(tài)監(jiān)測(cè)研究，如徐涵秋[4]利用歸一化差異水體指數(shù)（MNDWI）提取了水體信息；李生生[5]等提出了基于Landsat8 OLI 數(shù)據(jù)的青海湖水體邊界自動(dòng)提取方法；DENG Y[6]等基于1987—2015 年的Land?sat 衛(wèi)星影像研究了武漢城市群的湖泊范圍時(shí)空變化；RAO P Z[7]等基于MODIS產(chǎn)品對(duì)長(zhǎng)江流域地表水動(dòng)態(tài)變化進(jìn)行了分析；FENG L[8]等利用MODIS 數(shù)據(jù)開(kāi)展了2000—2010 年鄱陽(yáng)湖湖水淹沒(méi)變化監(jiān)測(cè)與評(píng)價(jià)；WANG X[9]等基于1984—2018 年的Landsat 衛(wèi)星圖像，采用SMDPSO方法提取了白洋淀水域面積，并進(jìn)行了時(shí)序變化分析。綜上所述，水體提取方法主要包括水體提取指數(shù)模型（NIR、NDVI、NDWI、SWI、MSWI）、監(jiān)督分類(lèi)（最大似然法、神經(jīng)網(wǎng)絡(luò)、SVM）、非監(jiān)督分類(lèi)、面向?qū)ο蟮取Ｉ鲜龇椒ㄔ诘?、中、高分辨率遙感影像上提取水體均有較好的效果[10-12]。雖然關(guān)于青藏高原地區(qū)水體的研究很多，但仍存在一些不足。由于青藏高原地區(qū)地形復(fù)雜，遙感影像上的陰影容易被誤分為水體，降低了現(xiàn)有方法水體提取的精度；且現(xiàn)有關(guān)于水資源變化的研究大多僅針對(duì)個(gè)別區(qū)域，對(duì)青藏高原其他廣大地區(qū)的變化知之甚少，二者之間的關(guān)系尚不明確。針對(duì)青藏高原地區(qū)水體和氣候的大尺度、多時(shí)相、全方位研究較少，分析區(qū)域水體的時(shí)空變化及其耦合響應(yīng)關(guān)系的研究更是少之又少[13-14]。

為研究青藏高原長(zhǎng)序列地表水體覆蓋面積動(dòng)態(tài)變化特征及其與氣候因子的關(guān)系，本文搜集了青藏高原地區(qū)2000—2020 年Landsat5、7、8 衛(wèi)星TM、ETM+、OLI 遙感影像以及相關(guān)氣候數(shù)據(jù)，提出了一種顧及地形的水體提取方法。首先提取青藏高原近20年的水體面積，并分析其時(shí)空變化特征；再利用統(tǒng)計(jì)學(xué)分析方法，探尋其動(dòng)態(tài)變化的氣候驅(qū)動(dòng)因素，為青藏高原保護(hù)和可持續(xù)發(fā)展提供合理的理論依據(jù)。

1 研究區(qū)概況與研究方法

1.1 研究區(qū)概況

青藏高原屬亞洲內(nèi)陸高原，是中國(guó)最大、世界海拔最高的高原，被稱(chēng)為“世界屋脊”、“第三極”[15]。青藏高原南起喜馬拉雅山脈南緣，北至昆侖山、阿爾金山和祁連山北緣，西部為帕米爾高原和喀喇昆侖山脈，東和東北部與秦嶺山脈西段和黃土高原相接，介于26°00′～39°47′N(xiāo)、73°19′～104°47′E之間；東西長(zhǎng)約為2 800 km，南北寬約為300～1 500 km，總面積約為250萬(wàn)km2，發(fā)育有大量冰川[16]，被譽(yù)為“亞洲水塔”（圖1）。青藏高原是世界海拔最高的一個(gè)巨型構(gòu)造地貌單元，具有獨(dú)特的自然環(huán)境和空間分異規(guī)律，受大氣環(huán)流和高原地勢(shì)格局的制約，形成了獨(dú)特的水熱狀況地域組合，呈現(xiàn)出由東南溫暖濕潤(rùn)向西北寒冷干旱的變化。高原腹地年均溫度在0℃以下，大片地區(qū)最暖月平均溫度也不足10℃。

圖1 青藏高原地區(qū)地形圖

1.2 數(shù)據(jù)來(lái)源

1）Landsat 數(shù)據(jù)。本文采用Landsat5 TM、Land?sat7 ETM+、Landsat8 OLI 衛(wèi)星遙感數(shù)據(jù)提取青藏高原水體。Landsat 數(shù)據(jù)來(lái)源于Google Earth Engine，數(shù)據(jù)獲取時(shí)間為2000—2020年。根據(jù)遙感數(shù)據(jù)質(zhì)量和云量情況，本文共選取100 675期影像數(shù)據(jù)。

2）氣象數(shù)據(jù)。通過(guò)綜合對(duì)比多個(gè)平臺(tái)的數(shù)十個(gè)氣候因子數(shù)據(jù)，本文在國(guó)家青藏高原科學(xué)數(shù)據(jù)中心（http://data.tpdc.ac.cn/zh-hans/）下載青藏高原及其周邊200 km緩沖區(qū)的氣候因子，包括溫度、降水量、蒸散量，所下載數(shù)據(jù)均為年值。其中，溫度、降水量數(shù)據(jù)覆蓋2000—2017年，蒸散量數(shù)據(jù)覆蓋2000—2015年。

1.3 研究方法

1.3.1 顧及地形效應(yīng)的水體提取方法

青藏高原地形起伏變化很大，為減少?gòu)?fù)雜地形條件對(duì)地物提取的干擾，本文提出了顧及地形的水體提取方法。詳細(xì)方法步驟為：

1）為減少?gòu)?fù)雜地形條件對(duì)水體提取的干擾，采用SCS+C模型對(duì)青藏高原地區(qū)進(jìn)行山區(qū)遙感影像地形校正。SCS+C 模型是Scott 于2005 年提出的模型，用于解決SCS 校正中的過(guò)校正問(wèn)題。SCS+C 模型是以SCS 模型為基礎(chǔ)，參照C 模型的建立方法構(gòu)建的[17]，表達(dá)式為：

式中，ρn為地形校正后的反射率；ρ為地形校正前的反射率；α為像元所在平面的坡度角；θ為太陽(yáng)天頂角；i為太陽(yáng)入射角；C為經(jīng)驗(yàn)參數(shù)。

2）采用改進(jìn)的MNDWI 提取各年度影像中的水體，再結(jié)合JRC水體數(shù)據(jù)集去除誤分像元，計(jì)算得到每個(gè)像元被判斷為水體像元的頻率，并輸出各像元被判定為水體的頻率。

式中，ρgreen為T(mén)M/ETM遙感數(shù)據(jù)中的綠光波段數(shù)據(jù)；ρSWIR為中紅外波段數(shù)據(jù)。

3）精度驗(yàn)證。選取兩景影像計(jì)算MNDWI，并以0為閾值分割水體和非水體；然后在影像中選取若干個(gè)隨機(jī)點(diǎn)，以目視的方式判斷隨機(jī)點(diǎn)處分類(lèi)是否正確，統(tǒng)計(jì)分類(lèi)正確和錯(cuò)誤的個(gè)數(shù)，計(jì)算總體精度。

4）選取合理閾值，對(duì)裁剪好的圖像進(jìn)行灰度分割。將提取的水體分為非水體、季節(jié)性水體和永久性水體3 類(lèi)，由各類(lèi)水體面積隨時(shí)間的變化可分析其整體變化趨勢(shì)。

1.3.2 水體與氣候因子相關(guān)性分析方法

1）氣候因子統(tǒng)計(jì)分析。通過(guò)創(chuàng)建時(shí)間序列的方式，采用中心移動(dòng)平均值函數(shù)（跨度設(shè)為2）對(duì)溫度、降水量、蒸散量數(shù)據(jù)進(jìn)行預(yù)處理，以盡可能地消除趨勢(shì)和周期性影響，變成平穩(wěn)性序列。

2）水體與氣候因子相關(guān)性分析。分析不同類(lèi)型水體與3 個(gè)氣候因子（溫度、降水量和蒸散量）的相關(guān)性；并對(duì)兩個(gè)組合分別進(jìn)行多元線(xiàn)性回歸，當(dāng)相關(guān)性顯著度較高時(shí)，對(duì)兩個(gè)變量做一元線(xiàn)性回歸，建立相關(guān)函數(shù)模型，探究其響應(yīng)關(guān)系。

本文研究的技術(shù)路線(xiàn)如圖2所示。

圖2 技術(shù)路線(xiàn)圖

2 實(shí)驗(yàn)與結(jié)果分析

2.1 水體提取與精度評(píng)價(jià)

本文分別利用MNDWI和顧及地形效應(yīng)的MNDWI水體提取方法，對(duì)青藏高原地區(qū)2000—2020年的水體提取結(jié)果進(jìn)行對(duì)比實(shí)驗(yàn)分析。本文選取兩景Landsat影像對(duì)水體提取精度進(jìn)行了驗(yàn)證，首先對(duì)兩景影像進(jìn)行預(yù)處理；然后計(jì)算MNDWI，并以0為閾值區(qū)分水體與非水體；最后在影像中選取若干個(gè)隨機(jī)點(diǎn)，通過(guò)人工目視解譯判斷隨機(jī)點(diǎn)處的分類(lèi)是否正確，并統(tǒng)計(jì)分類(lèi)正確和錯(cuò)誤的個(gè)數(shù)，計(jì)算總體精度。

水體提取精度如表1 所示，可以看出，SCS+C地形校正前，利用MNDWI 提取水體的總體精度分別為96.037 5%和94.308 4%；SCS+C 地形校正后，MNDWI 提取水體的總體精度分別為96.656 2%和95.454 5%，分別提升了0.618 7%和1.146 1%，說(shuō)明顧及地形的MNDWI 方法的分類(lèi)精度較高，通過(guò)地形校正，被誤分為水體的陰影像元被去除，精度得到改善。

表1 水體提取精度評(píng)價(jià)/%

2.2 水體時(shí)空變化分析

2000—2020 年青藏高原地區(qū)水體空間分布如圖3所示，可以看出，其水體分布廣泛，但受氣候和地形等因素的影響，存在明顯的地域性差異；總體上青藏高原地表水體呈現(xiàn)北多南少、西多東少的特點(diǎn)，主要發(fā)育在長(zhǎng)江、黃河、瀾滄江、怒江和雅魯藏布江5 條河流。

圖3 2000—2020年青藏高原地區(qū)水體分類(lèi)圖

首先利用本文提出的水體提取方法計(jì)算得到水體頻率；然后通過(guò)查閱文獻(xiàn)資料和歷史經(jīng)驗(yàn)，以頻率小于25%為非水體、介于25%～75%為季節(jié)性水體、大于75%為永久性水體為標(biāo)準(zhǔn)，對(duì)裁剪好的圖像進(jìn)行灰度分割；最后將提取的水體分為非水體、季節(jié)性水體和永久性水體3類(lèi)[18]。

本文對(duì)青藏高原地區(qū)3 種類(lèi)型水體隨時(shí)間的面積變化情況進(jìn)行了分析，結(jié)果如圖4 所示（圖中僅對(duì)季節(jié)性水體和永久性水體進(jìn)行統(tǒng)計(jì)分析），可以看出，2000—2020 年永久性水體面積由約4.46 萬(wàn)km2變化為約5.50 萬(wàn)km2，呈波動(dòng)上升趨勢(shì)；季節(jié)性水體面積由約0.93 萬(wàn)km2變化為約0.76 萬(wàn)km2，呈波動(dòng)下降趨勢(shì)；總體上2000年青藏高原水體面積（包含季節(jié)性水體和永久水體）約為5.39 萬(wàn)km2，約占青藏高原總面積的2.14%，2020 年水體面積（包含季節(jié)性水體和永久水體）約為6.26 萬(wàn)km2，約占青藏高原總面積的2.48%，水體面積總體呈波動(dòng)上升趨勢(shì)。

圖4 2000—2020年青藏高原地區(qū)水體面積年際變化柱狀圖

2.3 水體與氣候因子的相關(guān)性分析

1）青藏高原地區(qū)的氣候因子統(tǒng)計(jì)分析。根據(jù)年度統(tǒng)計(jì)平均值，將收集的青藏高原地區(qū)的氣候因子數(shù)據(jù)生成溫度、降水量、蒸散量的年際變化圖，如圖5～7 所示，可以看出，青藏高原受到了顯著的氣候變化影響，年均氣溫呈明顯的波動(dòng)上升趨勢(shì)；降水量、蒸散量年際變化呈現(xiàn)起伏變化，擬合曲線(xiàn)較為平穩(wěn)，特別是蒸散量年度變化情況較為穩(wěn)定。

圖5 青藏高原溫度年際變化圖

圖6 青藏高原降水量年際變化圖

圖7 青藏高原蒸散量年際變化圖

2）青藏高原地區(qū)水體變化與氣候的響應(yīng)關(guān)系。本文分析了青藏高原地區(qū)各類(lèi)水體面積與溫度、降水量和蒸散量的相關(guān)性，結(jié)果如表2～4 所示。從影響因子的顯著性分析可知，溫度對(duì)永久性水體的顯著性小于0.05，相關(guān)性較強(qiáng)，對(duì)季節(jié)性水體的相關(guān)性較弱。對(duì)溫度和季節(jié)性水體、永久性水體作一元線(xiàn)性回歸，可以得到溫度與永久性水體滿(mǎn)足關(guān)系式y(tǒng)=0.291x+5.097，即溫度越高、永久性水體面積越大；溫度與季節(jié)性水體滿(mǎn)足關(guān)系式y(tǒng)= -0.001x+0.866，即溫度越高、季節(jié)性水體面積越小。

表2 非水體與氣候的相關(guān)性

表3 季節(jié)性水體與氣候的相關(guān)性

表4 永久性水體與氣候的相關(guān)性

綜上所述，在青藏高原氣候暖化方向發(fā)展的背景下，水體面積變化與氣候因子具有顯著的區(qū)域相關(guān)性，蒸散量和降水量對(duì)水體面積變化的影響不大，相關(guān)性較低，溫度對(duì)水體面積變化的影響較大。青藏高原氣溫每升高1℃，永久水體面積增加0.291 萬(wàn)km2，季節(jié)性水體面積減少0.001萬(wàn)km2。

3 結(jié) 語(yǔ)

本文研究了2000—2020年青藏高原地區(qū)水體的時(shí)空變化特征及其與氣候因子的相關(guān)性，得出的主要結(jié)論為：

1）2000—2020年青藏高原地區(qū)水體面積總體呈波動(dòng)上升趨勢(shì)。其中，永久性水體面積由約4.46萬(wàn)km2變化為約5.50萬(wàn)km2，季節(jié)性水體面積由約0.93萬(wàn)km2變化為約0.76萬(wàn)km2。

2）通過(guò)對(duì)2000—2020 年青藏高原水體面積時(shí)空變化與溫度、降水量、蒸散量進(jìn)行相關(guān)性分析發(fā)現(xiàn)，水體面積變化與溫度的關(guān)聯(lián)性最大，與降雨量、蒸散量的相關(guān)性較??；溫度與永久性水體呈顯著正相關(guān)關(guān)系，與季節(jié)性水體呈負(fù)相關(guān)關(guān)系，即溫度越高、永久性水體面積越大、季節(jié)性水體面積越小。

在今后的研究工作中，將更加深入探索青藏高原地區(qū)水體信息的智能化提取方法以及氣候變化預(yù)測(cè)新技術(shù)，以期為青藏高原地區(qū)生態(tài)環(huán)境保護(hù)治理提供更加有力的技術(shù)支撐。