基于多源遙感數(shù)據(jù)的1970—2020 年巴爾喀什湖動(dòng)態(tài)監(jiān)測(cè)

張玉杰， 王寧練,3， 楊雪雯， 邰雪楠

（1.陜西省地表系統(tǒng)與環(huán)境承載力重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710127；2.西北大學(xué)城市與環(huán)境學(xué)院/地表系統(tǒng)與災(zāi)害研究院,陜西 西安 710127；3.中國(guó)科學(xué)院青藏高原地球科學(xué)卓越創(chuàng)作中心，北京 100101）

湖泊對(duì)氣候與環(huán)境變化極為敏感［1］，在維持區(qū)域生態(tài)平衡中起著重要作用［2］，湖泊的水位、面積和水量變化與區(qū)域氣溫、降水、蒸發(fā)及人類活動(dòng)等因素密切相關(guān)［3］，因而其動(dòng)態(tài)變化可間接反映區(qū)域氣候變化以及人類活動(dòng)對(duì)湖泊生態(tài)環(huán)境的影響［4-5］。巴爾喀什湖是中亞干旱區(qū)第三大水體，20 世紀(jì)80、90年代，巴爾喀什湖流域內(nèi)不合理的人類活動(dòng)使其生態(tài)環(huán)境遭到嚴(yán)重破壞，湖水水位大幅下降、湖泊面積萎縮、土地生產(chǎn)力下降［6］。2008年，聯(lián)合國(guó)開發(fā)計(jì)劃署在地理年鑒中表明巴爾喀什湖面臨著完全干涸的可能性［7］，隨后該湖泊的生態(tài)系統(tǒng)持續(xù)吸引著公眾的關(guān)注，因而對(duì)巴爾喀什湖及其環(huán)境的長(zhǎng)期監(jiān)測(cè)至關(guān)重要。

湖泊動(dòng)態(tài)變化主要體現(xiàn)在湖泊面積和水位的變化。傳統(tǒng)的湖泊水位監(jiān)測(cè)依靠水文站點(diǎn)實(shí)測(cè)，這種方法雖然精度較高，但對(duì)人力和物力的消耗較大［8］，且受限于自然條件與空間分布等因素，地域偏遠(yuǎn)的內(nèi)陸湖泊難以獲得長(zhǎng)時(shí)間的連續(xù)觀測(cè)資料。此外，傳統(tǒng)水位監(jiān)測(cè)共享程度較低，增加了數(shù)據(jù)獲取的難度［9］。衛(wèi)星測(cè)高技術(shù)作為1種空間大地測(cè)量技術(shù)，現(xiàn)被廣泛應(yīng)用于內(nèi)陸湖泊的水位監(jiān)測(cè)及分析其對(duì)環(huán)境的響應(yīng)［10］，與傳統(tǒng)的水位監(jiān)測(cè)方式相比，衛(wèi)星測(cè)高技術(shù)能夠大范圍、周期性地監(jiān)測(cè)陸地水體的各種動(dòng)態(tài)變化［11］，多項(xiàng)研究結(jié)果表明衛(wèi)星測(cè)高技術(shù)能夠有效監(jiān)測(cè)湖泊水位的動(dòng)態(tài)變化［12-13］。衛(wèi)星測(cè)高數(shù)據(jù)主要來(lái)源于激光與雷達(dá)2種類型的高度計(jì)。ICESat-1 是世界上第一顆星載激光測(cè)高衛(wèi)星，它利用地球科學(xué)激光測(cè)高系統(tǒng)測(cè)量全球地表高度，其地面足跡點(diǎn)的直徑約為70 m，沿軌跡相鄰足跡點(diǎn)間距約為172 m，能夠?qū)此贿M(jìn)行高精度的測(cè)量［8,14-15］。CryoSat-2 衛(wèi)星由ESA（European Space Agency）于2010 年發(fā)射，搭載的是合成孔徑干涉雷達(dá)高度計(jì)SIRAL，沿軌跡相鄰足跡點(diǎn)間距約為300 m，赤道上的相鄰軌跡間距約為8 km，與傳統(tǒng)的雷達(dá)測(cè)高系統(tǒng)如ERS-1/2、Envisat、GFO以及Jason-1/2等相比，CryoSat-2 軌道更密集，足跡點(diǎn)間距與沿軌間距更小，可以監(jiān)測(cè)更多的小型內(nèi)陸湖泊［16］。ICESat-1的運(yùn)行時(shí)間為2003—2009年，CryoSat-2的運(yùn)行時(shí)間為2010年至今，因而有多項(xiàng)研究結(jié)合這兩顆衛(wèi)星的觀測(cè)數(shù)據(jù)構(gòu)建了長(zhǎng)時(shí)間序列的觀測(cè)資料。宋春橋等［17-18］利用ICESat-1 與CryoSat-2 衛(wèi)星構(gòu)建了青藏高原131 個(gè)湖泊2003—2014 年的長(zhǎng)時(shí)間水位序列，且用ICESat-1 與CryoSat-2 衛(wèi)星對(duì)納木錯(cuò)湖2003—2014 年的水位趨勢(shì)變化進(jìn)行了進(jìn)一步的分析；黎鵬等［8］監(jiān)測(cè)了2003—2017 年洞庭湖的水位變化情況。

巴爾喀什湖水位數(shù)據(jù)實(shí)測(cè)記錄始于1938年，更早之前的延長(zhǎng)資料由前蘇聯(lián)學(xué)者通過野外實(shí)地考察、水量平衡估算和相關(guān)資料考證等多種方式獲得［19］，目前有很多研究采用了巴爾喀什湖水位數(shù)據(jù)。楊川德［20］利用延長(zhǎng)水位資料與實(shí)測(cè)資料討論了1879—1986年巴爾喀什湖水位變化及原因；龍愛華等［21］將巴爾喀什湖1879—2009 年的水位數(shù)據(jù)與該時(shí)段內(nèi)徑流量以及氣溫降水等數(shù)據(jù)相結(jié)合，對(duì)巴爾喀什湖水量平衡進(jìn)行了研究；郭利丹等［22］討論了巴爾喀什湖1810—2006 年的水位周期變化與流域氣候變化特征的關(guān)系；Propastin 等［23-24］利用衛(wèi)星測(cè)高資料和流域氣象站氣象資料，研究了1992—2010年巴爾喀什湖水位變化及其與氣候變化的關(guān)系，并建立了湖泊水位與伊犁河流量監(jiān)測(cè)的簡(jiǎn)易模型。

目前對(duì)于巴爾喀什湖面積的監(jiān)測(cè)主要采用遙感技術(shù)手段。臧菁菁等［25］依據(jù)水體的光譜反射特性，以9個(gè)時(shí)段的Landsat影像為數(shù)據(jù)源進(jìn)行水體信息提取，討論了巴爾喀什湖1975—2014年面積變化特征及其驅(qū)動(dòng)因素；成晨等［26］采用同季相4期Landsat遙感影像，通過目視解譯提取包括巴爾喀什湖在內(nèi)的中亞地區(qū)7 個(gè)湖泊1978—2010 年的面積變化；高彥華等［27］利用Landsat遙感影像，監(jiān)測(cè)了巴爾喀什湖1977—2014年面積及湖岸的動(dòng)態(tài)變化。

以往對(duì)于巴爾喀什湖的研究主要以單一遙感數(shù)據(jù)源為基礎(chǔ)，結(jié)合氣象數(shù)據(jù)分析湖泊近期水文變化，對(duì)于面積、水位及水量的長(zhǎng)期連續(xù)變化缺乏較為全面和系統(tǒng)的分析，且目前大部分水位研究都是采用水文站點(diǎn)的實(shí)測(cè)數(shù)據(jù)或是基于Jason-1、T/P 等高度計(jì)數(shù)據(jù)，在數(shù)據(jù)精度與時(shí)間尺度上都不夠完善，數(shù)據(jù)源較為單一，并未將傳統(tǒng)數(shù)據(jù)、激光以及雷達(dá)測(cè)高數(shù)據(jù)相融合，得到更加連續(xù)完整的長(zhǎng)時(shí)間觀測(cè)數(shù)據(jù)。因此，本文利用ICESat-1、CryoSat-2 衛(wèi)星測(cè)高數(shù)據(jù)及水文站實(shí)測(cè)數(shù)據(jù)，構(gòu)建巴爾喀什湖1970—2020年的水位序列，分析其水位的長(zhǎng)期和季節(jié)變化特征，并結(jié)合MOD09A1 數(shù)據(jù)及相關(guān)資料分析巴爾喀什湖長(zhǎng)期面積和水量變化特征，同時(shí)簡(jiǎn)要探究氣候變化及人類活動(dòng)對(duì)湖泊動(dòng)態(tài)變化的驅(qū)動(dòng)作用，為巴爾喀什湖的湖泊變化長(zhǎng)期監(jiān)測(cè)、水資源合理利用開發(fā)及生態(tài)系統(tǒng)保護(hù)提供一定的科學(xué)依據(jù)。

1 研究區(qū)概況

巴爾喀什湖（73°21′～79°30′E,44°45′～46°44′N）位于哈薩克斯坦東南部［20］，是典型的平原尾閭湖（圖1），湖泊面積寬廣，約為1.83×104km2，最寬處約71 km，總長(zhǎng)約600 km［27］。伊犁河作為巴爾喀什湖的主要水動(dòng)脈經(jīng)伊犁河三角洲注入湖泊西部，對(duì)入湖徑流量的貢獻(xiàn)率約為78%［28］；湖東部入湖徑流少，且蒸發(fā)量大，因而東部湖泊礦化程度明顯高于西部湖泊［22］，該湖北岸水體很少，地勢(shì)較高，多為巖石覆蓋；南岸小型湖沼較多，地勢(shì)低凹，多為沙地［25］。伊犁河—巴爾喀什湖流域是世界上最大的湖泊生態(tài)系統(tǒng)之一［21］，流域內(nèi)匯入湖泊的主要河流有伊犁河、卡拉塔爾河、列普西河以及阿亞古茲河，由于來(lái)自印度洋、太平洋的暖濕氣流無(wú)法深入流域所處的歐亞大陸腹地，因此巴爾喀什湖流域呈現(xiàn)溫帶大陸性氣候，干燥少雨。

圖1 巴爾喀什湖流域示意圖Fig.1 Location map of the Balkhash Lake Basin

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

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

2.1.1ICESat-1 測(cè)高數(shù)據(jù) ICESat-1/GLAS（Geoscience laser altimeter system）發(fā)射于2003 年1 月13日［29］，通過地球科學(xué)激光測(cè)高系統(tǒng)（GLAS）向星下點(diǎn)發(fā)射激光信號(hào)，根據(jù)激光信號(hào)在衛(wèi)星與星下點(diǎn)之間往返所用的時(shí)間，來(lái)計(jì)算二者之間的距離，進(jìn)而獲得星下點(diǎn)的高程［8］，該衛(wèi)星獲取的信息時(shí)間范圍為2003 年2 月20 日至2009 年10 月10 日，共有640 多天的觀測(cè)數(shù)據(jù)［30］，地學(xué)激光測(cè)高系統(tǒng)的數(shù)據(jù)產(chǎn)品共有15 種，本文采用的是GLA14 全球測(cè)高數(shù)據(jù)產(chǎn)品，該數(shù)據(jù)集中包含有激光測(cè)高腳點(diǎn)的經(jīng)緯度、高程以及各項(xiàng)校正參數(shù)等。

2.1.2CryoSat-2 測(cè)高數(shù)據(jù) CryoSat-2 于2010 年4月8 日發(fā)射升空，該衛(wèi)星的測(cè)量模式分別為低分辨率模式（LRM）、合成孔徑雷達(dá)模式（SAR）以及合成孔徑雷達(dá)干涉測(cè)量模式（SIRAL）［31］，CryoSat-2/SIRAL有0級(jí)數(shù)據(jù)、1b級(jí)數(shù)據(jù)、2級(jí)數(shù)據(jù)和高級(jí)數(shù)據(jù)4種數(shù)據(jù)產(chǎn)品，其中2 級(jí)GDR 產(chǎn)品包括了測(cè)量時(shí)間、地理位置、以及經(jīng)過了儀器校正、傳輸延遲改正、幾何改正和地球物理改正之后的高度信息，是單獨(dú)的全球性的全軌道數(shù)據(jù)，即將3 種模式的測(cè)量結(jié)果通過不同的處理過程，按照時(shí)間先后順序進(jìn)行整合，得到格式統(tǒng)一的數(shù)據(jù)記錄。3種模式的數(shù)據(jù)采用了不同的波形重跟蹤算法來(lái)獲得高度值，分別為Refined CFI、UCL 和Refined OCOG［32］，通過對(duì)比3 種算法在巴爾喀什湖湖區(qū)內(nèi)獲得高度異常值的數(shù)目，最終采用GDR 產(chǎn)品的Refined OCOG 重跟蹤算法來(lái)獲取巴爾喀什湖2010—2020年的湖泊水位變化。

2.1.3湖泊面積數(shù)據(jù) 中分辨率成像光譜儀MODIS搭載于TERRA和AQUA衛(wèi)星上，是為全球動(dòng)態(tài)測(cè)量提供大范圍數(shù)據(jù)的重要傳感器，MODIS共有36個(gè)觀測(cè)波段且具有較高的分辨率，因而被廣泛用于水體特征信息的提取，目前已有多項(xiàng)研究表明MODIS數(shù)據(jù)在湖泊面積的監(jiān)測(cè)中具有較高的精度。MODIS數(shù)據(jù)產(chǎn)品極為豐富，其中MOD09（來(lái)源于https://ladsweb.modaps.eosdis.nasa.gov/search/）為MODIS陸地2級(jí)數(shù)據(jù)產(chǎn)品，本文選用的是MOD09 每8 d 合成的數(shù)據(jù)，即MOD09A1，空間分辨率為500 m。

2.1.4氣象數(shù)據(jù) 氣象資料選擇了英國(guó)East Anglia大學(xué)Climatic Research Unit（CRU）發(fā)布的CRUTS v4.05［CRU Time series version 4.05，獲取地址：CRU TS Version 4.05(uea.ac.uk)］。該數(shù)據(jù)集是由CRU 收集整理多種數(shù)據(jù)庫(kù)中的相關(guān)資料，重構(gòu)的一套覆蓋全球范圍的、高分辨率的連續(xù)地表氣候要素?cái)?shù)據(jù)集［33］，數(shù)據(jù)的空間分辨率為0.5°×0.5°，已有研究表明該數(shù)據(jù)可用于討論中亞地區(qū)的氣候變化狀況［34］，由于巴爾喀什湖附近水文站點(diǎn)公布的氣象資料不完整，因此本文使用該數(shù)據(jù)集1970—2019年的氣溫和降水量數(shù)據(jù)作為主要?dú)庀笠蜃觼?lái)分析巴爾喀什湖水位對(duì)氣候變化的響應(yīng)（2020年數(shù)據(jù)尚未公布），同時(shí)，選擇了來(lái)源于NOAA 全球氣象站點(diǎn)數(shù)據(jù)網(wǎng)站（https://www.climate.gov/maps-data）的部分站點(diǎn)數(shù)據(jù)來(lái)驗(yàn)證CRUTS v4.05數(shù)據(jù)集在巴爾喀什湖流域的適用性。

2.1.5其他數(shù)據(jù) 巴爾喀什湖湖泊矢量邊界數(shù)據(jù)來(lái)源于全球湖泊和濕地?cái)?shù)據(jù)庫(kù)（Global Lakes and Wetlands Database，GLWD）1 級(jí)數(shù)據(jù)。巴爾喀什湖1970年起的實(shí)測(cè)水位數(shù)據(jù)來(lái)源于HYDROLARE（http://hydrolare.net/），其中1998 年與1999 年數(shù)據(jù)存在部分月份的缺失。研究中使用的巴爾喀什湖流域灌溉面積來(lái)源于聯(lián)合國(guó)糧食及農(nóng)業(yè)組織FAO（http://www.fao.org/aquastat/zh/countries-and-basins/regionaloverviews/central-asia/）和相關(guān)研究資料。

2.2 研究方法

本文的技術(shù)路線如圖2 所示，首先對(duì)ICESat-1 GLA14 和CryoSat-2 GDR 測(cè)高數(shù)據(jù)進(jìn)行預(yù)處理，提取衛(wèi)星的足跡點(diǎn)高程，篩選出巴爾喀什湖邊界范圍內(nèi)的湖面足跡點(diǎn)，剔除高程異常的湖面足跡點(diǎn)后計(jì)算得到湖泊2003—2009年的衛(wèi)星測(cè)高水位，并根據(jù)實(shí)測(cè)水位數(shù)據(jù)進(jìn)行精度驗(yàn)證，結(jié)合實(shí)測(cè)數(shù)據(jù)以及2種測(cè)高數(shù)據(jù)提取的水位信息得到1970—2020 年巴爾喀什湖的水位序列，水位數(shù)據(jù)的具體組成結(jié)構(gòu)如圖3 所示。此外，采用MOD09A1 數(shù)據(jù)提取2000—2020年湖泊面積數(shù)據(jù)，結(jié)合相關(guān)文獻(xiàn)資料得到長(zhǎng)時(shí)序湖泊面積數(shù)據(jù)，最后根據(jù)公式，估算得到巴爾喀什湖的水量變化。

圖2 技術(shù)路線Fig.2 Flowchart of synthesizing multi-source lake-level data in this study

圖3 水位數(shù)據(jù)組成結(jié)構(gòu)Fig.3 Data composition structure diagram

2.2.1水位數(shù)據(jù)提取

（1）衛(wèi)星測(cè)高數(shù)據(jù)預(yù)處理

ICESat-1 測(cè)高數(shù)據(jù)的預(yù)處理主要包括橢球轉(zhuǎn)換和飽和校正。由于ICESat-1 與CryoSat-2 測(cè)高數(shù)據(jù)基于不同的參考橢球，因此需要將ICESat-1測(cè)高數(shù)據(jù)的參考橢球（T/P 橢球）轉(zhuǎn)換成WGS84 橢球，以消除橢球差異帶來(lái)的影響［8］。此外，ICESat-1 在進(jìn)行高程測(cè)量時(shí)可能因波形飽和現(xiàn)象導(dǎo)致測(cè)量足跡點(diǎn)的高程低于實(shí)際高程值，因此需要對(duì)其進(jìn)行飽和校正。足跡點(diǎn)高程的計(jì)算原理如公式（1）所示：

式中：H為基于EGM96大地水準(zhǔn)面的正高；h為基于T/P 橢球的高程；d_deltaEllip 為T/P 橢球和WGS84橢球之間的差值；d_satElevCorr 為飽和校正參數(shù)，可從數(shù)據(jù)集里獲得；N為當(dāng)?shù)卮蟮厮疁?zhǔn)面差距，可由MATLAB中的geoidheight函數(shù)計(jì)算得到。

CryoSat-2 測(cè)高數(shù)據(jù)的預(yù)處理主要包括各項(xiàng)偏差改正。CryoSat-2 衛(wèi)星雷達(dá)脈沖信號(hào)在傳播過程中發(fā)生的散射或折射會(huì)影響其傳播速度，使觀測(cè)信號(hào)的往返時(shí)間產(chǎn)生延遲，同時(shí)受各種自然因素的影響，距離估計(jì)值會(huì)產(chǎn)生偏差，因此必須加以改正［35］。湖泊水位的計(jì)算公式如下：

式中：H為基于EGM96 大地水準(zhǔn)面的正高；Halt為衛(wèi)星質(zhì)心到參考橢球面的高度；R為衛(wèi)星到湖面的距離；?R為各項(xiàng)誤差改正；N為當(dāng)?shù)卮蟮厮疁?zhǔn)面差距，可由MATLAB中的geoidheight函數(shù)計(jì)算得到。

式中：Dry 為干對(duì)流層改正；Wet 為濕對(duì)流層改正；Ion為電離層改正；Sol為固體潮改正；Pol為極潮改正。

（2）高程異常點(diǎn)移除

數(shù)據(jù)預(yù)處理完成后，對(duì)獲得的所有水位值進(jìn)行處理，以提取湖泊的水位序列，具體的移除步驟如下［5,36-37］：

第一步，對(duì)湖泊邊界向湖心方向做200 m 的緩沖區(qū)，根據(jù)緩沖區(qū)篩選衛(wèi)星足跡點(diǎn)，以確保數(shù)據(jù)點(diǎn)完全落入湖內(nèi)，減少可能與湖岸接觸的高程點(diǎn)對(duì)單天水位數(shù)據(jù)的干擾。

第二步，對(duì)所得的水位數(shù)據(jù)點(diǎn)先進(jìn)行目視解譯，將與大多數(shù)水位值相差幾十米甚至幾百米偏差極大的異常值剔除。

第三步，用3σ準(zhǔn)則剔除單天水位數(shù)據(jù)中的異常值，然后將剩余有效水位值取平均作為日均水位。3σ準(zhǔn)則的具體判別方法如公式（4）所示：對(duì)于采集的數(shù)據(jù)樣本（x1,x2,x3,…,xn），求取算數(shù)平均值xˉ和剩余誤差值vi，即可求得均方根偏差σ：

以ICESat-1在2004年5月23日的過境點(diǎn)為例，經(jīng)3σ準(zhǔn)則剔除異常點(diǎn)之后，高于343.4 m的水位值被剔除，剩余水位點(diǎn)在343.17～343.39 m 之間波動(dòng)（圖4）。

圖4 3σ 準(zhǔn)則剔除異常值前后對(duì)比Fig.4 Comparison before and after removing outliers under the 3σ criterion

第四步，對(duì)所有的日均水位，先目視解譯剔除明顯離群值，再用3σ準(zhǔn)則，進(jìn)一步剔除異常值，最后計(jì)算得到月均水位與年均水位。

（3）數(shù)據(jù)轉(zhuǎn)換

如技術(shù)路線圖（圖2）所示，ICESat-1 數(shù)據(jù)采用了T/P、EGM2008 參考系統(tǒng)，CryoSat-2 數(shù)據(jù)采用了WGS84、EGM96 參考系統(tǒng)，而水文站的實(shí)測(cè)數(shù)據(jù)采用波羅的海（Baltic）高程系統(tǒng)，3 種數(shù)據(jù)采用了不同的高程基準(zhǔn)，因此在構(gòu)建水位序列前需要將其進(jìn)行融合。在本研究中，首先，將ICESat-1 數(shù)據(jù)轉(zhuǎn)換為與CryoSat-2相同的參考系統(tǒng)，通過上文的高程異常點(diǎn)移除方法進(jìn)行數(shù)據(jù)篩選，獲得巴爾喀什湖2003—2020 年湖泊年均水位序列。根據(jù)該水位序列與對(duì)應(yīng)年份實(shí)測(cè)水位間的平均差值，將高度計(jì)獲得的水位減去平均差值轉(zhuǎn)換為與實(shí)測(cè)水位相同的水準(zhǔn)面［36］，進(jìn)而得到巴爾喀什湖1970—2020 年的水位序列。

（4）精度驗(yàn)證

本文參照前人方法，通過相關(guān)系數(shù)（r）、顯著性值（P），平均絕對(duì)誤差（MAE）以及均方根誤差（RMSE）來(lái)計(jì)算雷達(dá)高度計(jì)提取的湖泊水位與對(duì)應(yīng)的實(shí)測(cè)水位間的相關(guān)性［8,17］。相關(guān)系數(shù)（r）是用以反映變量之間相關(guān)關(guān)系密切程度的統(tǒng)計(jì)指標(biāo)。

式中：r為相關(guān)系數(shù)；xi為高度計(jì)獲得的水位值；xˉ為高度計(jì)獲得的水位值的算數(shù)平均值；yi為實(shí)測(cè)水位值；yˉ為實(shí)測(cè)水位值的算數(shù)平均值；n為樣本數(shù)。r>0 時(shí)表明變量之間為正相關(guān)，反之為負(fù)相關(guān)，r的絕對(duì)值越接近于1，表明變量間相關(guān)關(guān)系越強(qiáng)。

顯著性值（P）是用以反映某一事件發(fā)生的可能性大小的統(tǒng)計(jì)指標(biāo)，P<0.05 為顯著，P<0.01 為非常顯著，表明樣本間的差異由抽樣誤差所致的概率小于0.05或0.01。

平均絕對(duì)誤差（MAE）表示提取值和實(shí)測(cè)值之間絕對(duì)誤差的平均值，均方根誤差（RMSE）是提取值與實(shí)測(cè)值偏差的平方與樣本個(gè)數(shù)比值的平方根，2個(gè)指標(biāo)都能夠很好地反映出提取值與實(shí)測(cè)值之間的偏差。

2.2.2面積數(shù)據(jù)提取 巴爾喀什湖1970、1975、1980、1985、1990、1995 年面積數(shù)據(jù)來(lái)源于相關(guān)文獻(xiàn)［19,25］，2000—2020 年逐年遙感影像面積是基于9月MOD09A1 數(shù)據(jù)，采用Mcfeeters［38］提出的歸一化差異水體指數(shù)（Normalized difference water index，NDWI）提取湖泊邊界計(jì)算所得，NDWI 表示如下：

式中：Green 為綠色波段，即MOD09A1 第4 波段；NIR 為近紅外波段即MOD09A1 第2 波段。利用NDWI可提取水體信息，減弱地表土壤、植被等信息特性的干擾進(jìn)行閾值分割［39］，有效地將水體同植被與土壤明顯區(qū)分開。

2.2.3水量變化估計(jì) 根據(jù)內(nèi)陸湖的水量平衡，內(nèi)陸湖的水量變化是由湖區(qū)面積和水位共同決定的。我們使用以下方程來(lái)估算水量變化［40］：

式中：?V為2 個(gè)時(shí)期的蓄水量變化；?H為2 個(gè)時(shí)期的湖泊水位變化；S1和S2為2 個(gè)時(shí)期的湖區(qū)面積。

3 結(jié)果與分析

3.1 巴爾喀什湖水位變化特征

結(jié)合ICESat-1、CryoSat-2 2 種衛(wèi)星測(cè)高數(shù)據(jù)和水文站實(shí)測(cè)數(shù)據(jù)，得到了巴爾喀什湖1970—2020年湖泊水位變化情況（圖5）。在研究時(shí)段內(nèi)，湖泊水位呈先下降后上升趨勢(shì)。1970—1987年，湖泊水位顯著下降（P<0.01），其變化率約為1.23 m·（10a）-1，水位總體下降約2.21 m，且1987年湖泊水位達(dá)到最低值（340.64 m）。1988—2020年湖泊水位總體呈上升趨勢(shì)（P<0.01），在此期間，可以觀察到多次較短的湖泊水位振蕩，增長(zhǎng)區(qū)間分別為：1988—1990 年、1992—1995 年、1999—2007 年、2009—2011 年、2015—2020年，其余時(shí)段巴爾喀什湖水位呈小幅下降趨勢(shì)。其中1999—2007 年湖泊的變化趨勢(shì)與Propastin 等［23］基于TOPEX/Poseidon 和Jason-1 衛(wèi)星在巴爾喀什湖的測(cè)高結(jié)果基本吻合，該時(shí)段為20世紀(jì)90年代后湖泊水位的顯著上升期，持續(xù)上升時(shí)間最長(zhǎng)。2009—2011 年上升速率最高，為1.66 m·(10a)-1，局部峰值達(dá)到342.66 m。

圖5 1970—2020年巴爾喀什湖水位變化Fig.5 Changes of water level in Balkhash Lake from 1970 to 2020

為進(jìn)一步監(jiān)測(cè)巴爾喀什湖水位的年內(nèi)動(dòng)態(tài)變化，本文對(duì)研究時(shí)段內(nèi)各月的水位進(jìn)行了分析（圖6）。結(jié)果表明，巴爾喀什湖的水位變化具有明顯的季節(jié)性特征，1—6 月多年水位均值呈上升趨勢(shì)（P<0.01），且6 月達(dá)到年內(nèi)水位峰值，約為341.99 m，隨后湖泊水位呈下降趨勢(shì)（P<0.01），進(jìn)入9月后，各月平均水位高度值波動(dòng)幅度逐漸變小，開始趨于穩(wěn)定，水位均值介于341.65～341.7 m，說(shuō)明巴爾喀什湖的年內(nèi)水位增長(zhǎng)主要發(fā)生在2 月底至6 月初，且暖季（4—10月）水位波動(dòng)比冷季（11月—次年3月）更劇烈。

圖6 1970—2020年巴爾喀什湖逐月水文變化Fig.6 Monthly hydrological changes of Balkhash Lake from 1970 to 2020

3.2 巴爾喀什湖面積變化特征

1970—2020 年巴爾喀什湖面積的長(zhǎng)期變化如圖7a 所示。1970—2020 年，巴爾喀什湖面積從19996 km2減 少 到16641.93 km2，平 均 變 化 率 為-65.77 km2·a-1（P<0.01），湖面縮減約16.77%，面積最大值出現(xiàn)在1970 年（19996 km2），最小值出現(xiàn)在1990 年（16638.87 km2）。1970—1990 年，巴爾喀什湖的湖泊面積以-159.86 km2·a-1的速度大幅萎縮（P<0.01）。進(jìn)入21世紀(jì)后，2000—2010年湖泊面積有所回升，由16672.89 km2增加到了16896.44 km2，增加約0.013%，面積增加速率約為20.32 km2·a-1（P<0.01），而2010 年之后，湖泊面積呈波動(dòng)減少的趨勢(shì)，萎縮率約為21.31 km2·a-1（P<0.01），至2020 年，巴爾喀什湖面積約為16641.93 km2。

圖7 1970—2020年巴爾喀什湖面積變化Fig.7 Changes in the area of Balkhash Lake from 1970 to 2020

3.3 巴爾喀什湖水量變化特征

本文根據(jù)巴爾喀什湖面積數(shù)據(jù)與水位數(shù)據(jù)，計(jì)算得到巴爾喀什湖多年水量變化（圖8）。1970—2020年巴爾喀什湖水量損失約為12.33 km3，1970—1985 年，湖泊蓄水量以1.84 km3·a-1的速率持續(xù)下降，累積變化量為-27.55 km3，1985—1990年蓄水量大幅增長(zhǎng)，增長(zhǎng)量達(dá)10.39 km3，1995—2007 年湖泊蓄水量以1.86 km3·a-1的速率持續(xù)上升，隨后水量變化呈波動(dòng)下降趨勢(shì)，其中2011—2015年蓄水量減少速率最快（-2.65 km3·a-1），累積減少約10.59 km3。

圖8 1970—2020年巴爾喀什湖水量變化Fig.8 Changes of water volume in Balkhash Lake from 1970 to 2020

3.4 水位數(shù)據(jù)精度驗(yàn)證

ICESat-1 和CryoSat-2 衛(wèi)星的觀測(cè)數(shù)據(jù)與水位站實(shí)測(cè)數(shù)據(jù)在2003—2016 年存在部分重合，其中，ICESat-1 獲取時(shí)間主要集中在2—6 月和9—12 月，與實(shí)測(cè)數(shù)據(jù)的重合值有27 個(gè)，CryoSat-2 為逐月觀測(cè)，從2010 年7 月起，與實(shí)測(cè)數(shù)據(jù)的重合值個(gè)數(shù)為73。因此，為驗(yàn)證衛(wèi)星測(cè)高數(shù)據(jù)的準(zhǔn)確性和可靠性，對(duì)衛(wèi)星觀測(cè)數(shù)據(jù)與實(shí)測(cè)水位數(shù)據(jù)做相關(guān)性分析，結(jié)果顯示（圖9），巴爾喀什湖的ICESat-1 和CryoSat-2 衛(wèi)星測(cè)高水位與水文站實(shí)測(cè)水位之間呈顯著的正相關(guān)關(guān)系，相關(guān)系數(shù)分別為0.874和0.805，顯著性值均小于0.01，且ICESat-1 數(shù)據(jù)與實(shí)測(cè)數(shù)據(jù)相比平均絕對(duì)誤差（MAE）及均方根誤差（RMSE）分別為0.206和0.241，CryoSat-2數(shù)據(jù)與實(shí)測(cè)數(shù)據(jù)相比平均絕對(duì)誤差（MAE）與均方根誤差（RMSE）分別為0.153 和0.199，結(jié)果表明利用2 種衛(wèi)星測(cè)高數(shù)據(jù)來(lái)進(jìn)行巴爾喀什湖長(zhǎng)時(shí)間序列的水位變化監(jiān)測(cè)具有一定的可行性。

圖9 精度驗(yàn)證Fig.9 Precision validation

3.5 影響因素分析

3.5.1氣候變化 巴爾喀什湖流域內(nèi)水文站點(diǎn)公布的氣象資料不完整，僅有1970—2005年的年氣溫與降水量數(shù)據(jù)，為驗(yàn)證CRUTS v4.05 氣象數(shù)據(jù)集在巴爾喀什湖流域內(nèi)的適用性，本文對(duì)流域站點(diǎn)數(shù)據(jù)與CRUTS v4.05 氣象數(shù)據(jù)進(jìn)行了相關(guān)性分析，結(jié)果顯示（圖10），氣象站點(diǎn)獲得的氣溫與降水量數(shù)據(jù)與CRUTS v4.05 數(shù)據(jù)集之間呈顯著的正相關(guān)關(guān)系，相關(guān)系數(shù)分別為0.965 和0.951。表明CRUTS v4.05 數(shù)據(jù)集可用于討論巴爾喀什湖流域的氣候變化狀況。

圖10 CRUTS v4.05數(shù)據(jù)集精度驗(yàn)證Fig.10 Precision validation of CRUTS v4.05 dataset

1970—2019 年巴爾喀什湖流域多年平均氣溫約為6.02 ℃，流域氣溫在研究時(shí)段內(nèi)整體呈顯著上升趨勢(shì)，通過了99%置信度檢驗(yàn)，變化率為0.39 ℃·(10a)-1（圖11a）。M-K 檢驗(yàn)的結(jié)果表明巴爾喀什湖流域的年均溫在1988年發(fā)生突變，開始出現(xiàn)增溫趨勢(shì)，2003 年以后流域的年均溫上升趨勢(shì)顯著，置信度超過95%（圖11b）。巴爾喀什湖流域多年平均降水量為280.584 mm，其變化率為0.06 mm·（10a）-1，未通過90%置信度檢驗(yàn)（圖11c）。M-K 檢驗(yàn)結(jié)果表明巴爾喀什湖流域的年降水量出現(xiàn)多個(gè)突變點(diǎn)，降水量在1978—1995 年波動(dòng)幅度較大，而在1997 年之后開始增多，2015 年之后呈顯著增加趨勢(shì)，置信度達(dá)到95%（圖11d）。

圖11 1970—2020年巴爾喀什湖流域氣候變化特征Fig.11 Characteristics of climate change in Balkhash Lake Basin from 1970 to 2020

對(duì)巴爾喀什湖流域氣溫及降水量變化與湖泊動(dòng)態(tài)變化分時(shí)段進(jìn)行相關(guān)性分析，結(jié)果表明水位、面積及水量變化在1970—1989年與氣溫、降水量均無(wú)相關(guān)關(guān)系，1990—2019 年，僅水位與氣溫呈正相關(guān)，相關(guān)性為0.387（在95%水平上顯著）。因而，巴爾喀什湖流域氣溫及降水量變化與湖泊水位、面積及水量的年際變化在長(zhǎng)期趨勢(shì)上并無(wú)明顯相關(guān)性。

3.5.2人類活動(dòng) 1970 年卡普恰蓋水庫(kù)建成后，伊犁河入湖水量由1950—1969年的148.2×108m3銳減到1970—1985年的116.1×108m3［41-42］，同時(shí)隨著流域灌溉農(nóng)業(yè)的發(fā)展，湖區(qū)及周邊地區(qū)的農(nóng)業(yè)用水大幅增長(zhǎng)［43-44］，最終使得巴爾喀什湖水位、面積以及水量在1970—1987 年均持續(xù)下降。1986 年起哈薩克斯坦有關(guān)部門終止水庫(kù)蓄水計(jì)劃、調(diào)節(jié)水庫(kù)發(fā)電量并改善河道淤積狀況，這一系列措施使得巴爾喀什湖生態(tài)環(huán)境有所改善。1991年蘇聯(lián)解體，哈薩克斯坦經(jīng)濟(jì)受創(chuàng)，國(guó)家終止了對(duì)農(nóng)場(chǎng)的支持，實(shí)施土地私有化，并開始征收水費(fèi)，巴爾喀什湖流域的灌溉面積因此直線下降［21］，許多依靠電力驅(qū)動(dòng)水泵的灌溉區(qū)被完全放棄，農(nóng)業(yè)用水的減少進(jìn)一步促使水位緩慢上升［23］，面積與水量隨之呈現(xiàn)出回升趨勢(shì)，2000—2006 年，流域灌溉面積介于大約2000～3000 km2之間，與20世紀(jì)90年代相比大幅減少（圖12）［45］，在此期間，巴爾喀什湖的水位有了明顯的回升趨勢(shì)。因而，巴爾喀什湖流域內(nèi)強(qiáng)烈的人類活動(dòng)（包括卡普恰蓋水庫(kù)蓄水和哈薩克斯坦境內(nèi)灌區(qū)耕地?cái)U(kuò)張等）是1970—1987年湖泊動(dòng)態(tài)變化的主要影響因素，同時(shí)，人類活動(dòng)也在湖泊動(dòng)態(tài)變化的長(zhǎng)期趨勢(shì)中發(fā)揮著重要的影響作用。

圖12 1970—2006年巴爾喀什湖流域灌溉面積Fig.12 Irrigated area in the Balkhash Lake Basin from 1970 to 2006

以往的研究結(jié)果表明［21,41,44-46］，巴爾喀什湖動(dòng)態(tài)變化受氣溫、降水量、蒸發(fā)量、冰川融水以及人類活動(dòng)等多種因素的綜合影響，定量地探討各因素在不同時(shí)期湖泊動(dòng)態(tài)變化中的驅(qū)動(dòng)作用，系統(tǒng)地掌握和理解巴爾喀什湖動(dòng)態(tài)變化的過程及原因，還需做更進(jìn)一步的研究。

4 結(jié)論

（1）據(jù)ICESat-1與CryoSat-2的監(jiān)測(cè)結(jié)果，巴爾喀什湖水位于1970—1987年快速下降，平均變化速率為1.23 m·（10a）-1，1987 年至今，巴爾喀什湖水位總體呈現(xiàn)上升趨勢(shì)，并伴隨有周期性波動(dòng)。同時(shí)水位具有較為規(guī)律的季節(jié)性變化，春季多年平均水位最高，秋季最低，水位增長(zhǎng)主要發(fā)生在湖泊結(jié)冰期（11 月底到次年4 月初），暖季（4—10 月）水位波動(dòng)比冷季（11月—次年3月）更劇烈；研究時(shí)段內(nèi)，巴爾喀什湖面積呈現(xiàn)出先明顯萎縮后略有恢復(fù)的變化過程，51 a湖面縮減約16.77%；在此期間，巴爾喀什湖水量損失約為12.33 km3。

（2）巴爾喀什湖動(dòng)態(tài)變化受氣候變化與人類活動(dòng)的共同影響，1970—1987 年水位、面積及水量變化主要由卡普恰蓋水庫(kù)蓄水和哈薩克斯坦境內(nèi)灌區(qū)耕地?cái)U(kuò)張等人類活動(dòng)引起，但研究時(shí)段內(nèi)，巴爾喀什湖的動(dòng)態(tài)變化與氣候變化在總體趨勢(shì)上無(wú)明顯相關(guān)性，定量描述二者之間的聯(lián)系，還需進(jìn)一步探討。

（3）測(cè)高衛(wèi)星為湖泊水位監(jiān)測(cè)提供了有效的技術(shù)手段，與常規(guī)水位測(cè)量方法相比，測(cè)高衛(wèi)星不受自然環(huán)境及人力因素的限制，可覆蓋更多湖泊，整合不同衛(wèi)星的測(cè)高數(shù)據(jù)，還可獲得較長(zhǎng)時(shí)序的湖泊水位信息，這對(duì)研究湖泊水位的動(dòng)態(tài)變化及其驅(qū)動(dòng)因素具有重要意義。