朱凝嫣，肖作林，夏鴻萱

(重慶師范大學(xué)地理信息系統(tǒng)應(yīng)用研究重慶高校市級(jí)重點(diǎn)實(shí)驗(yàn)室，重慶 401331)

0 引言

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

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

目前對(duì)黃河流域長(zhǎng)時(shí)序水體變化的研究有一定的開展[9-11][21]，但多集中于某一流域或某一地區(qū)，如黃土高原地區(qū)[12]，黃河灘區(qū)[13]等，且未能揭示黃河流域豐水期和枯水期的水體面積差異特征。鑒于此，本文以黃河流域?yàn)檠芯繀^(qū)域，基于Google earth engine平臺(tái)，針對(duì)不同河段、不同季節(jié)遴選最優(yōu)水體識(shí)別方法，分析黃河流域1990—2020年間黃河流域不同河段、不同典型流域、不同季節(jié)的水體變化特征，采用累積量斜率變化率分析法探究人類活動(dòng)和降水對(duì)水體面積變化的貢獻(xiàn)率。研究結(jié)論可以準(zhǔn)確把控黃河流域水體在時(shí)間及空間上的分布和變化趨勢(shì)，對(duì)黃河流域復(fù)雜狀況的綜合管理、水利建設(shè)、河流規(guī)劃、水旱監(jiān)測(cè)、農(nóng)業(yè)發(fā)展以及水土保持都具有重要意義。

1 研究區(qū)概況和數(shù)據(jù)來源

1.1 研究區(qū)概況

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

圖1 研究區(qū)概況圖

1.2 數(shù)據(jù)來源

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

1901—2021年中國1km分辨率逐月降水量和氣溫?cái)?shù)據(jù)集源于國家科技基礎(chǔ)條件平臺(tái)—國家地球系統(tǒng)科學(xué)數(shù)據(jù)中心(http：//www.geodata.cn)，見表1。

表1 數(shù)據(jù)來源

2 研究方法

2.1 基于大津算法的改進(jìn)后歸一化差異水體指數(shù)的計(jì)算

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

(1)

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

2.2 基于閾值的多波段譜間關(guān)系法

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

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

(2)

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

表2 多波段譜間關(guān)系法閾值選擇

2.3 混淆矩陣計(jì)算

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

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

2.4.1累積距平法

累積距平法是由曲線直觀判斷離散數(shù)據(jù)點(diǎn)變化趨勢(shì)的一種非線性統(tǒng)計(jì)方法：對(duì)于自變量x，其在某一時(shí)刻t的累計(jì)距平表示為：

(3)

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

1 生物標(biāo)志物的概念 生物標(biāo)志物是一種用于客觀測(cè)量和評(píng)價(jià)正常生物過程、致病過程或機(jī)體對(duì)治療干預(yù)的生物反應(yīng)的指示物，可反映機(jī)體所做出的可被測(cè)定的改變[4]。由于其可以早期、快速、準(zhǔn)確判斷疾病發(fā)生與發(fā)展情況，因此在新藥研發(fā)中發(fā)揮著重要作用。

2.4.2累積量斜率變化率比較法計(jì)算貢獻(xiàn)率[19]

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

Rw=100×(Kb/Ka-1)

(4)

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

Rp=100×(Pb/Pa-1)

(5)

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

Cp=100×Rp/Rw

(6)

則人類活動(dòng)對(duì)水體面積變化的貢獻(xiàn)率(CH，單位：%)可以表示為：

CH=1-Cp

(7)

3 結(jié)果與分析

3.1 水體提取精度驗(yàn)證

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

3.2 黃河流域水體面積變化規(guī)律

3.2.1黃河流域水體面積隨時(shí)間變化規(guī)律

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

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

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

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

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

3.2.3黃河流域典型水體面積動(dòng)態(tài)變化

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

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

圖4 黃河流域典型水體面積動(dòng)態(tài)變化

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

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

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

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

圖6 黃河流域氣溫及降水動(dòng)態(tài)變化

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

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

圖7 黃河流域不同河段降水與人類活動(dòng)對(duì)水體面積變化貢獻(xiàn)率

4 結(jié)論

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

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

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

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