近30年黃河流域水體面積變化與影響因素分析

朱凝嫣，肖作林，夏鴻萱

(重慶師范大學地理信息系統應用研究重慶高校市級重點實驗室，重慶 401331)

0 引言

作為中國的母親河，黃河的水量水質變化一直都是社會關注的重點和熱點。整個黃河流域水文特征復雜且多變，上游源頭部分雖儲水量豐富但依舊有時令河的存在，中游流經的黃土高原是黃河泥沙和污染物的主要來源區，這就導致了下游淤泥堆積河床抬高形成懸河，河勢游蕩多變，防洪形勢嚴峻[1]。黃河地表水的變化及其與洪水事件，與人類活動以及長期氣候變化的聯系十分緊密[2]，而這恰恰最能體現在地表水體面積的變化上。近年來，遙感技術在研究長時間段黃河地表水的變化趨勢上發揮著越來越重要的作用。

水體提取方法的歷史可以追溯至上個世紀70年代，幾十年來各國學者致力于提高水體提取精度和效率，不斷更新優化水體提取方法。早在1977年Bartolucci L.A.等人就提出MSS波段提取出的水體效果最好[3]。1985年，Jupp等人開始使用LandsantTM 7波段通過閾值方法提取水體[4]。RundquistD等人在1987年曾使用單波段(近紅外波段)閾值方法提取水體[5]。Mcfeeters[6]在1995年提出了歸一化差分水體指數(NDWI)，這一水體指數曾在長時間內在世界范圍內得到廣泛使用，而后徐涵秋基于NDWI提出來改良的歸一化差異水體指數(MNDWI)，MNDWI能夠用于更好地揭示水體的一些細小特征，這些細微之處就包括水中懸浮物和沉積物的分布以及水質發生的一些小變化等[7]。另外，MNDWI還可以有效識別陰影信息對水體的影響，一定程度上解決了水體提取中陰影難以消除這一難題[8]。

目前對黃河流域長時序水體變化的研究有一定的開展[9-11][21]，但多集中于某一流域或某一地區，如黃土高原地區[12]，黃河灘區[13]等，且未能揭示黃河流域豐水期和枯水期的水體面積差異特征。鑒于此，本文以黃河流域為研究區域，基于Google earth engine平臺，針對不同河段、不同季節遴選最優水體識別方法，分析黃河流域1990—2020年間黃河流域不同河段、不同典型流域、不同季節的水體變化特征，采用累積量斜率變化率分析法探究人類活動和降水對水體面積變化的貢獻率。研究結論可以準確把控黃河流域水體在時間及空間上的分布和變化趨勢，對黃河流域復雜狀況的綜合管理、水利建設、河流規劃、水旱監測、農業發展以及水土保持都具有重要意義。

1 研究區概況和數據來源

1.1 研究區概況

黃河全長大約5464km，水面最高處和最低處相差可達約4480m，流域總面積約79.5萬km2，其中內流區約4.2萬km2。黃河流域跨度廣，地勢差異大。河口鎮分界了黃河上游與中游，桃花峪分界了黃河中游與下游。黃河流域水資源分布大致是由南向北逐漸減少的態勢。黃河支流眾多，但一級支流在黃河左岸和右岸的分布并不平均，右岸支流的數量明顯多于左岸支流的數量。黃河源頭地區水資源豐富，密布支流和湖泊，濕地種類多樣，年降水量充沛，來自高原永久冰川融化的水資源充足，黃河60%的降水量集中在6—9月，年內分配不均且年際變化大。在上世紀由于生態環境惡化，城鄉生產生活用水，各種自然原因與人文因素影響，黃河開始出現頻繁的季節性斷流，以1997年最為嚴重。

圖1 研究區概況圖

1.2 數據來源

本文采用Landsat-5 TM影像、Landsat 7 ETM+和Landsat 8 OLI數據提取1990—2020年間黃河地表水體，采用JRC(Joint Research Centre)水體數據集優化水體提取結果。JRC數據集包含1984—2021年地表水的位置和時間分布圖，并提供了有關這些水面范圍和變化的統計數據。為了檢驗水體提取結果精度，本文選用Sentinel-2 MSI遙感數據，在GEE平臺上進行去云掩膜處理以獲得高質量的影像。Sentinel-2衛星的重訪周期為10d成像時間更短，其遙感數據的空間分辨率為10m，圖像空間精度更高。

1901—2021年中國1km分辨率逐月降水量和氣溫數據集源于國家科技基礎條件平臺—國家地球系統科學數據中心(http：//www.geodata.cn)，見表1。

表1 數據來源

2 研究方法

2.1 基于大津算法的改進后歸一化差異水體指數的計算

為了提取黃河流域豐水期(7—9月)的水體，采用改進的歸一化差異水體指數(MNDWI)提取水體[14]。傳統歸一化差異水體指數(NDWI)主要利用綠光波段和近紅外波段進行計算，水體對近紅外波段的吸收強、植被對其反射率強，通過突出水體抑制植被來提取影像中水體信息。但NDWI提取的地物信息經常與建筑和土壤混合，導致水體面積值偏高。MNDWI利用短波紅外波段替代近紅外波段。該指數顯著改善了開闊水域特征，快速、準確地區分水體和非水體特征，在很大程度上降低或者消除建筑物的影響[15]。MNDWI是劃分開闊水域特征最廣泛、最有效的方法之一。其計算公式如下：

(1)

式中，Green—綠光波段；SWIR—短波紅外波段。大津算法(otsu)常用于計算機視覺和圖像處理中，通過返回單個強度的閾值將像素分為前景和背景兩類。本研究使用大津算法自動提取水體和非水體的分割閾值。當MNDWI值大于分割閾值時，將像素點劃分為有效水體；當水體頻率小于分割閾值時，將像素點劃分為非水體[16]。

2.2 基于閾值的多波段譜間關系法

對于黃河流域枯水期(12月—次年2月)內的水體，采用MNDWI指數法提取水體會出現大范圍的山體陰影誤判的現象。因此，為了降低陰影誤判提升整體精度，本文在枯水期則采用基于閾值的多波段譜間關系法[17]來進行水體提取。該方法可以增大水體和陰影灰度值的差異，從而將兩者很好地區分開來。水體的提取模型如下所示：

(TM2+TM3)-(TM4+TM5)>N(水體)

(2)

式中，N—閾值，閾值的選取可以通過直方圖來大致確定，然后通過反復實驗，直至得到最佳閾值。由于本文研究時間跨度長，涉及不同的數據來源，因此基于不同數據源與不同河段選擇的最佳閾值也不同(見表2)

表2 多波段譜間關系法閾值選擇

2.3 混淆矩陣計算

基于Sentinel-2 MSI遙感數據，選擇了2018年12月至2019年2月黃河下游影像數據與2017年7—9月黃河中游的影像數據，對所提取的水體結果進行精度驗證。通過生產者精度(PA)、用戶精度(UA)、總體精度(OA)和Kappa系數(Kc)計算混淆矩陣進行精度評價。

2.4 累積量斜率變化率比較法

2.4.1累積距平法

累積距平法是由曲線直觀判斷離散數據點變化趨勢的一種非線性統計方法：對于自變量x，其在某一時刻t的累計距平表示為：

(3)

該方法的核心是判斷離散數據對其均值的離散幅度，若累積距平值增大，表明離散數據大于其平均值，反之則小于其平均值[18]。如果曲線由上述兩個部分組成，則可確定變化趨勢的拐點。

1 生物標志物的概念 生物標志物是一種用于客觀測量和評價正常生物過程、致病過程或機體對治療干預的生物反應的指示物，可反映機體所做出的可被測定的改變[4]。由于其可以早期、快速、準確判斷疾病發生與發展情況，因此在新藥研發中發揮著重要作用。

2.4.2累積量斜率變化率比較法計算貢獻率[19]

水體面積累積距平曲線斜率在拐點前后兩個時期分別為Ka和Kb，降水量累積距平曲線的斜率在拐點前后兩個時期分別為Pa和Pb，則累積水體面積斜率變化率(Rw，單位：%)為：

Rw=100×(Kb/Ka-1)

(4)

累積降水量斜率變化率(Rp，單位：%)為：

Rp=100×(Pb/Pa-1)

(5)

降水量變化對水體面積變化的貢獻率(Cp，單位：%)可以表示為：

Cp=100×Rp/Rw

(6)

則人類活動對水體面積變化的貢獻率(CH，單位：%)可以表示為：

CH=1-Cp

(7)

3 結果與分析

3.1 水體提取精度驗證

與Landsat-5/7/8影像數據相比，Sentine-2影像具有更高的時間與空間分辨率。本文將Sentinel-2與Landsat影像所得結果進行對比，為了確保影像質量，本文選擇了2018年12月至2019年2月黃河下游影像數據與2017年7—9月黃河中游的影像數據(均為已通過大氣校正的地表反射率影像)，隨機選取水體樣本點數據與非水體樣本點數據，采用混淆矩陣的方式來驗證水體提取精度。基于大津算法的改進后歸一化差異水體指數的提取結果Kappa系數為0.88，基于閾值的多波段譜間關系法提取結果Kappa系數為0.86，一致性檢驗結果較好，可以用作后續研究。

3.2 黃河流域水體面積變化規律

3.2.1黃河流域水體面積隨時間變化規律

圖2 黃河流域1990—2020年水體面積變化

總體來看，1990—2020年，黃河流域豐水期內水體面積增加了1419.58km2，枯水期水體面積增加了704.34km2。其中豐水期水體面積最小值出現在1998年，為7295.49km2；1990—1998年，豐水期水體面積減少了1822.67km2，1998—2020年，豐水期水體面積增加了3242.26km2；枯水期水體面積的最小值出現在1997年，為4881.70km2；1990—1997年，流域內枯水期水體面積減少了1573.79km2，1997—2020年，枯水期水體面積增加了2278.14km2。

3.2.2黃河流域水體面積空間分布變化規律

黃河流域不同河段水體面積如圖3所示，豐水期黃河上游地區水域面積平均為5445.27km2，約占整個流域陸表水域面積的比例為62.13%；中游地區水域面積平均為2055.61km2，約占整個流域陸表水域面積的比例為23.45%；花園口以下的下游地區水域面積平均為1264.07km2，約占整個流域陸表水域面積的比例為14.42%。枯水期黃河上游地區水域面積平均為4053.375km2，約占整個流域陸表水域面積的比例為65.23%；中游地區水域面積平均為1286.66km2，約占整個流域陸表水域面積的比例為20.71%；花園口以下的下游地區水域面積平均873.59km2，約占整個流域陸表水域面積的比例為14.06%。黃河上游地區豐水期與枯水期以及黃河下游地區豐水期水體面積均呈現先減少后增加的趨勢，上游豐水期與枯水期水體面積最小值均出現在2001年，分別為4613.715km2與3282.69km2，下游地區豐水期水體面積最小值出現在1998年，面積為989.24km2。黃河中游地區枯水期以及黃河下游地區枯水期面積呈現持續性緩慢增加的趨勢，而黃河中游地區豐水期的水體面積則變化不明顯。

圖3 1990—2020年黃河流域不同河段水體面積變化

3.2.3黃河流域典型水體面積動態變化

為研究黃河流域不同類型的水體的變化趨勢及其在豐水期與枯水期內的差異，本文選擇了鄂陵湖、龍羊峽水庫、黃河龍門段、小浪底水庫、東平湖湖口和山陽水庫6個代表性水體。特征水體分布在上中下游的不同區域，水體的選擇包含了自然水體(大型自然湖泊、黃河干流河段、自然湖泊出湖口)與人工水體(山地大型水利工程、中下游平原大型水庫、中小型水庫)。

如圖4—5所示，上游的鄂陵湖水面呈現持續擴張狀態，2000年后在湖泊西南角枯水期與豐水期面積差異逐漸增大；中游龍門河段水體面積處于一個不斷波動的過程，其最顯著的特征是枯水期與豐水期的面積差異極大，枯水期面積明顯小于豐水期；下游東平湖的出湖口逐漸收縮，東平湖與干流交匯處卻略有增加趨勢，同時可以看出交匯處干流河段的枯水期面積與豐水期面積之間差異也較為明顯。上游大型水庫龍羊峽水庫水體面積也呈現持續增大狀態，水庫面積在枯水期與豐水期差異極小，主要體現在西南部與黃河交匯處；位于平原地區的小浪底水庫，2000年之后水面突增，與水庫1999年開始蓄水有關，小浪底水庫豐水期與枯水期的面積差異較龍羊峽水庫更大，1990年水庫未建之時河道枯水期面積也略小于豐水期面積；黃河下游中型水庫山陽水庫，水面呈現持續減少狀態，尤以2000年最為明顯，1990年山陽水庫豐水期與枯水期差異較大，2000年之后水庫枯水期與豐水期面積只有細微的差別。通過對黃河流域內不同特征水體的分析，整體來看，枯水期與豐水期的水體面積差異主要集中在黃河干流與湖泊等自然水體，尤以中游干流最為明顯，反而人工水體如水庫等，則差異并不十分明顯。

圖4 黃河流域典型水體面積動態變化

3.3 黃河流域水體面積變化的影響因素及貢獻率分析

已有研究結果表明，地表水體面積的時空分布主要受自然因素和人類活動的影響[22]，而自然因素又多分為降水與氣溫兩類。根據月均降水量值與月均溫數據，選擇7—9月平均值代表豐水期降水量與豐水期氣溫，12月—次年2月的平均值代表枯水期降水量與枯水期氣溫。1990—2020年黃河流域降水與氣溫的變化趨勢如圖5所示，可見在黃河流域，豐水期降水值遠遠高于枯水期，枯水期降水量波動較豐水期更大。

圖5 黃河流域不同時段典型水體面積差異

由圖6可知，黃河流域上游河段，中游豐水期與下游豐水期氣溫在1990—2020年間的變化并不顯著，中游枯水期與下游枯水期氣溫波動較大但無明顯變化趨勢，因此在計算貢獻率之前對黃河流域不同河段豐水期面積與枯水期面積與流域內氣溫做相關分析得，無論枯水期還是豐水期，上游、中游與下游的水體面積與氣溫的相關性均不高(|R|均小于0.4)。因此在對黃河流域水體面積變化的驅動因素分析時，重點考慮降水這一因素。

圖6 黃河流域氣溫及降水動態變化

本文采用累積量斜率變化率比較法來計算降水對黃河流域水體面積變化的貢獻率。累積量斜率變化率比較法可以應用在干旱—半干旱地區河流變化及其影響因素的定量評估中[20]。通過式(3)，得出1990—2020年降水量與水體面積的累積距平值，通過式(4)，式(5)，式(6)可得出降水量在豐水期與枯水期對黃河流域地表水體面積變化的貢獻率，通過式(7)得出人類活動對水體面積的貢獻率。

如圖7所示，在黃河流域內，無論上游中游還是下游，豐水期降水量對地表水體面積變化的貢獻率都比枯水期更大，降水量對水體面積的影響更加明顯，這源于豐水期降水量比枯水期更加豐沛，而降水恰好是流域內地表水體最主要的補給來源，在枯水期內則是人類活動的影響占比更大。從空間上看，豐水期內在黃河上游降水對地表水體面積的影響最為明顯，貢獻率為74.12%；在黃河中游與下游地區相對來說影響較??；豐水期降水量對地表水體面積的影響最小的地區是黃河中游地區，貢獻率為52.19%。枯水期內在黃河中游區域人類活動對地表水體面積的影響最大，貢獻率為76.47%；在黃河上游人類活動對地表水體面積的影響最小，貢獻率為37.13%。黃河流域內三大自然湖泊，鄂陵湖、扎陵湖與東平湖分別位于上游與下游，同時黃河流域上游地區還擁有龍羊峽、劉家峽等大型水庫，對于降水有很好的調節與涵養作用。中游地區，由于黃土高原占據了流域內中游地區極大的面積，是整個流域內水土流失最為嚴重的地區，缺少大型湖泊，天然降水不能夠很好地轉變為地表水，因此中游地區地表水面積受人類活動影響更大。這些影響包括正向影響如淤地壩的建設，陜西、山西兩省的淤地壩數量占流域內淤地壩總數的89.70%；負向影響如生產生活用水與工業用水等，山西省煤炭開采極大程度影響了地表水體面積。下游地區河網分布較為稠密，城市內部水體多且黃汶區內有較大數量的中小型水庫，因此相比中游地區，自然因素的影響更大，但由于下游地區城市發展速度較中上游更為迅速，近年來人類活動強度與人口密度也迅速增大，導致下游水體變化受人類活動影響大。因此，在對黃河流域進行水資源保護和整體規劃時，可以重點考慮人類活動的影響，如耕地用水，城市建設用水與大型水利工程的建設等。

圖7 黃河流域不同河段降水與人類活動對水體面積變化貢獻率

4 結論

本文以黃河流域為研究區，將全年內時間段分為豐水期(7—9月)與枯水期(12月—次年2月)，采用MNDWI指數(利用大津算法自動分割閾值)和多波段譜間關系法，提取并分析了近30a黃河流域不同河段豐水期與枯水期的地表水體面積變化特征，采用累積量斜率變化分析法識別了研究區地表水體面積變化的主要影響因素。得到如下結論。

(1)研究區內水體主要分布在黃河流域上游，占比超過60%，下游水體面積最小，占比不到15%。豐水期水體面積最小值出現在1998年，為7295.489km2，枯水期的最小值出現在1997年，為4881.695km2，同一年內豐水期面積普遍大于枯水期面積。

(2)黃河流域內枯水期與豐水期的水體面積差異主要體現在黃河干流與湖泊等自然水體，尤以中游干流最為明顯；人工水體(如水庫)豐水期與枯水期的差異并不十分明顯。

(3)同一地區枯水期人類活動影響較豐水期更明顯，豐水期自然因素貢獻率更大。中游地區人類活動為地表水體面積變化的主要影響因素(76.47%)，上游地區則是降水為主導因素(74.12%)，下游地區人類活動與降水對地表水體面積貢獻率差異不大。