基于多時(shí)相主被動(dòng)遙感的漓江水面監(jiān)測與水質(zhì)參數(shù)反演(2016-2020年)*

劉 曼，付波霖，何宏昌，解淑毓，仇繼圣，孫習(xí)東，龔燁云，勞植楠，左萍萍

(桂林理工大學(xué)測繪地理信息學(xué)院，桂林 541006)

河流水文信息動(dòng)態(tài)監(jiān)測是分析流域水體演變規(guī)律的基礎(chǔ)，而水質(zhì)信息的采集則是評(píng)價(jià)流域水體質(zhì)量優(yōu)劣程度、衡量水體初級(jí)生產(chǎn)力的前提，及時(shí)準(zhǔn)確地監(jiān)測、獲取河流水文和水質(zhì)信息對河流水資源開發(fā)管理、水環(huán)境保護(hù)、水生態(tài)修復(fù)以及防洪減災(zāi)等方面具有重要意義[1-2]. 漓江是桂林風(fēng)景的靈魂，是中國唯一一條入選的全球最美河流[3]. 然而，近年來漓江水域水量暴漲、干涸等現(xiàn)象突出，影響了桂林人民的正常生活. 快速、準(zhǔn)確地掌握漓江水體信息，并依據(jù)流域時(shí)空演變特點(diǎn)研究相應(yīng)對策，對漓江流域進(jìn)行合理開發(fā)、規(guī)劃具有重大意義.

目前，遙感技術(shù)已被用于河流水文、水質(zhì)信息的監(jiān)測，相較于傳統(tǒng)的地面監(jiān)測方法，其在宏觀性、時(shí)效性、便利性上具有較大優(yōu)勢. Yang等[4]和王大釗等[5]使用歸一化差異水體指數(shù)(NDWI)、改進(jìn)歸一化水體指數(shù)(MNDWI)與自動(dòng)提取水體指數(shù)(AWEI)等方法，基于Sentinel-2和Landsat 8影像提取了城市地表水分布信息，結(jié)果表明Sentinel-2 影像提取的水體信息更為準(zhǔn)確，Kappa系數(shù)達(dá)0.88以上；段秋亞等[6]基于高分一號(hào)(GF-1)數(shù)據(jù)，采用NDWI閾值法，提取了 2 塊不同尺度和不同復(fù)雜度的代表性水域，與人工解譯的水體信息相比較，閾值法提取精度分別達(dá)到 97.59% 與99.15%；Yao等[7]利用城市水自動(dòng)提取方法和資源三號(hào)(ZY-3)影像提取了青島多個(gè)地區(qū)的地表水信息，與NDWI指數(shù)提取結(jié)果相比，該方法Kappa系數(shù)提高了7%～32%. 因光學(xué)影像易受到天氣的影響，部分學(xué)者則利用微波遙感全天時(shí)、全天候監(jiān)測以及后向散射系數(shù)對水體較敏感等特點(diǎn)來進(jìn)行水體提取研究. Hsu等[8]基于臺(tái)灣地區(qū)四景Sentinel-1雷達(dá)影像，利用后向散射系數(shù)提取了河流水面寬度，與無人機(jī)影像相比，利用雷達(dá)影像測量的結(jié)果易受陰影的干擾；Liang等[9]基于Sentinel-1合成孔徑雷達(dá)圖像，采用局部閾值化方法提取了Louisiana部分洪水淹沒地區(qū)的水體輪廓，與傳統(tǒng)閾值法相比，局部閾值法在區(qū)別水體與非水體方面精度更高，用戶精度與生產(chǎn)者精度提高了4%～13%；谷鑫志等[10]提出了水體自動(dòng)化提取算法，利用兩景不同成像模式的GF-3 SAR影像對湖南省東北部水體進(jìn)行信息提取，結(jié)果表明該算法提取精度達(dá)到了85%以上，精度優(yōu)于基于Landsat 8與GF-1影像的NDWI、MNDWI、AWEI指數(shù). 葉綠素濃度遙感反演研究對水體的管理起著積極作用，為驗(yàn)證遙感影像對流域水質(zhì)信息的監(jiān)測效果，Ansper等[11]和Toming等[12]利用Sentinel-2數(shù)據(jù)反演了Estonian湖泊葉綠素a(Chl.a)濃度，證明了Sentinel數(shù)據(jù)對水體中葉綠素的估算和湖泊時(shí)空動(dòng)態(tài)的跟蹤的適用性；Nazeer等[13]為監(jiān)測香港沿海復(fù)雜水域的Chl.a濃度，使用Landsat TM 、ETM+影像中藍(lán)、綠、紅和近紅外4個(gè)波段數(shù)據(jù)及實(shí)測葉綠素?cái)?shù)據(jù)建立了回歸模型，相關(guān)系數(shù)達(dá)0.89；Filipponi[14]利用C2RCC算法從Sentinel-2數(shù)據(jù)中提取意大利波河水體信息，并分析了水體中Chl.a和總懸浮物質(zhì)(TSM)濃度變化，證明了利用衛(wèi)星多光譜數(shù)據(jù)監(jiān)測河流水質(zhì)變化的可行性. 綜上所述，以上研究主要集中在比較不同遙感水體指數(shù)提取水體的精度，少有學(xué)者利用長時(shí)間序列的主被動(dòng)遙感影像提取水體信息和反演水體Chl.a、TSM濃度等水質(zhì)參數(shù)，定量分析河流水文和水質(zhì)信息的年內(nèi)動(dòng)態(tài)變化規(guī)律.

因此，本文以漓江流域?yàn)檠芯繀^(qū)，以2016-2020年Landsat 8 OLI、GF-1和Sentinel-2A多光譜數(shù)據(jù)和2018年逐月Sentinel-1A IW模式的SAR影像為主要數(shù)據(jù)源，利用構(gòu)建的多種主被動(dòng)遙感水體指數(shù)逐月提取漓江水面信息，反演漓江逐月的Chl.a和TSM濃度，構(gòu)建了漓江流域的278個(gè)基本評(píng)價(jià)單元，并在此單元內(nèi)基于河岸線發(fā)育率、水面變化區(qū)域差異值和水體信息變化動(dòng)態(tài)度定量分析漓江水面、Chl.a和TSM濃度的年內(nèi)、枯水期和汛期、上、中和下游的時(shí)空演化特征，為漓江的水環(huán)境保護(hù)、流域合理開發(fā)與保護(hù)提供技術(shù)支撐和理論依據(jù).

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

1.1 研究區(qū)概況

漓江流域(24°18′～25°41′N，109°45′～110°40′E)位于廣西壯族自治區(qū)桂林市境內(nèi)，漓江發(fā)源于興安縣華江鄉(xiāng)貓兒山老山界東南側(cè)，自北向南流經(jīng)興安縣、靈川縣、桂林市區(qū)、陽朔縣、荔浦縣、恭城縣、平樂縣等地，全長164 km，屬珠江水系西江的支流[15]. 桂林位于低緯度區(qū)域，屬亞熱帶濕潤季風(fēng)氣候，雨季集中在夏季期間. 漓江年徑流量大，但分布不均勻，水位全年變化顯著，3-8月為汛期，9月至次年2月為枯水期，且兩期差異顯著[16]. 枯水期河床外露，嚴(yán)重時(shí)船只無法通行；汛期水量暴漲，可淹沒兩岸河床和灘涂. 桂林山水甲天下，而漓江則是桂林山水的靈魂所在，漓江部分河段有著喀斯特地貌，是廣西東北部喀斯特地貌發(fā)育的最典型地區(qū)，象鼻山、冠巖、黃布倒影等地較為出名，1982年國務(wù)院將漓江流域列入“首批國家重點(diǎn)風(fēng)景名勝區(qū)”之一[17]，桂林旅游業(yè)得以快速發(fā)展. 近年來漓江水量暴漲、干涸現(xiàn)象突出，因受水位漲幅的劇烈影響，桂林旅游業(yè)呈現(xiàn)明顯的淡旺季，桂林人民的正常生活也受到嚴(yán)重影響，汛期時(shí)因洪水造成的經(jīng)濟(jì)損失尤為嚴(yán)重. 漓江流域位置如圖1所示.

圖1 漓江流域位置Fig.1 Location of Lijiang River Basin

1.2 數(shù)據(jù)源及數(shù)據(jù)預(yù)處理

為探究漓江水體年內(nèi)動(dòng)態(tài)變化特征，本文獲取了漓江流域全年逐月的主被動(dòng)遙感影像. 光學(xué)遙感影像主要包括：①2016-2020年12個(gè)月的Sentinel-2A多光譜數(shù)據(jù)，下載網(wǎng)址為https://scihub.copernicus.eu/；②2016-2018年11個(gè)月Landsat 8 OLI和GF-1 WFV多光譜數(shù)據(jù)，其中Landsat 8 OLI數(shù)據(jù)是從美國地質(zhì)調(diào)查局USGS網(wǎng)站(https://earthexplorer.usgs.gov/)下載獲取，因部分月份出現(xiàn)Landsat 8 OLI影像研究區(qū)覆蓋不全、云覆蓋率較高等問題，本文將從中國資源衛(wèi)星應(yīng)用中心網(wǎng)站(http://www.cresda.com)下載的對應(yīng)月份的GF-1 WFV多光譜數(shù)據(jù)作為補(bǔ)充數(shù)據(jù)，而在2016-2018年研究區(qū)3月的Landsat 8影像和GF-1影像均被較多云遮擋，無法進(jìn)行漓江水體的提取. 本文利用Sen2Cor v2.8插件對Sentinel-2A影像進(jìn)行輻射校正、大氣校正處理，再用SNAP軟件將多光譜影像重采樣為15 m的空間分辨率，并導(dǎo)入到ENVI 5.4軟件中對數(shù)據(jù)進(jìn)行波段合成、影像拼接和裁剪，生成研究區(qū)影像；而基于Sentinel-2A影像對Chl.a及TSM濃度的反演，則直接使用SNAP軟件對原數(shù)據(jù)重采樣，并進(jìn)行C2RCC大氣校正，將Salinity參數(shù)設(shè)置為0.0001；運(yùn)用ENVI 5.4軟件分別對Landsat 8 OLI多光譜影像進(jìn)行輻射定標(biāo)和FLAASH大氣校正，對全色波段進(jìn)行輻射校正，并利用Gram-Schmidt Spectral Sharpening算法將校正后的全色與多光譜進(jìn)行圖像融合，生成空間分辨率為15 m的多光譜影像，最后裁剪得到漓江流域影像數(shù)據(jù)；同時(shí)，利用ENVI 5.4軟件對GF-1 WFV影像進(jìn)行輻射定標(biāo)、大氣校正、幾何校正、影像拼接與裁剪.

本文在歐空局遙感數(shù)據(jù)中心網(wǎng)站(http://www.cnblogs.com/)下載了漓江流域2018年12個(gè)月IW模式SLC格式的Sentinel-1A影像. 運(yùn)用SNAP 軟件處理Sentinel-1A 影像，先后經(jīng)過熱噪聲去除、輻射定標(biāo)、Burst去除、Refined Lee濾波處理、多視處理、斜距影像轉(zhuǎn)換為地距影像、地形糾正、地理編碼和影像裁剪，導(dǎo)出空間分辨率為15 m為ENVI格式的后向散射系數(shù)數(shù)據(jù). 為了減少相干噪聲對準(zhǔn)確提取和解譯研究區(qū)域水體的影響，本文分別將Refined Lee濾波窗口設(shè)置為7×7，多視處理視數(shù)為4. 基于 30 m SRTM DEM作為參考DEM 對Sentinel-1A影像進(jìn)行地理編碼，生成投影坐標(biāo)系為 WGS_1984_UTM_Zone_49N. 研究區(qū)的主被動(dòng)遙感數(shù)據(jù)集以Sentinel-2A多光譜影像為基準(zhǔn)影像，其他數(shù)據(jù)在ENVI 5.4 軟件中進(jìn)行相對配準(zhǔn). 數(shù)據(jù)源具體信息如附表Ⅰ所示.

水體Chl.a和TSM濃度實(shí)測數(shù)據(jù)來自于廣西壯族自治區(qū)桂林市水文水資源局提供的2018年逐月數(shù)據(jù)，用于驗(yàn)證基于C2RCC算法和Sentinel-2A影像反演的逐月Chl.a和TSM濃度. 漓江水面寬度數(shù)據(jù)驗(yàn)證數(shù)據(jù)基于Google Earth中高空間分率遙感影像進(jìn)行解譯和矢量化提取. 其中，在解譯提取漓江水面寬度數(shù)據(jù)時(shí)，將Google Earth中影像調(diào)整為與本文所用遙感影像對應(yīng)的月份，并從漓江上、中和下游各選取5個(gè)漓江水面橫截面計(jì)算出平均值作為最終漓江水面寬度的驗(yàn)證數(shù)據(jù).

2 研究方法

為系統(tǒng)研究漓江水文和水質(zhì)信息年內(nèi)動(dòng)態(tài)變化特征，本文基于水利部《河流健康評(píng)估指標(biāo)、標(biāo)準(zhǔn)與方法 V1.0》，以500 m 為步長等分主航道中心線且生成覆蓋整個(gè)研究區(qū)的 278個(gè)矩形緩沖區(qū)，格網(wǎng)化漓江水體建立278個(gè)基本評(píng)價(jià)單元. 基于基本評(píng)價(jià)單元分別計(jì)算河岸線發(fā)育系數(shù)(SDI)、水面變化區(qū)域差異值(WDr)及水體信息變化動(dòng)態(tài)度(K)3種指標(biāo)參數(shù)，來定量分析漓江水面、Chl.a和TSM濃度的年內(nèi)、枯水期和汛期上、中和下游的變化規(guī)律. 本文采用閾值分割算法對漓江遙感水體指數(shù)進(jìn)行分割運(yùn)算，提取逐月的水體信息，分別計(jì)算了4類光學(xué)遙感水體指數(shù)和雷達(dá)后向散射系數(shù)：Landsat 8 OLI和GF-1多光譜影像分別計(jì)算了歸一化水體指數(shù)(NDWI)、改進(jìn)型歸一化水體指數(shù)(MNDWI)、增強(qiáng)型水體指數(shù)(EWI)，Sentinel-2A多光譜影像計(jì)算了歸一化差值池指數(shù)(NDPI)，Sentinel-1A SAR影像計(jì)算了后向散射系數(shù)(S). 綜合利用多種光學(xué)遙感水體指數(shù)和雷達(dá)后向散射系數(shù)進(jìn)一步構(gòu)建了主被動(dòng)遙感加權(quán)指數(shù)(JQ)，增強(qiáng)水體信息的同時(shí)抑制其他非水體信息. 本文采用二類水體區(qū)域性近岸海域水色算法(C2RCC)、最大葉綠素指數(shù)(MCI)、雙波段比值法(DoubleR)及Chl.a反射峰強(qiáng)度(ρchl)4種方法，反演計(jì)算漓江水體Chl.a和TSM濃度，進(jìn)而探究其年內(nèi)動(dòng)態(tài)變化規(guī)律，具體技術(shù)流程如圖2所示.

圖2 研究方法和技術(shù)路線Fig.2 Research method and technical route

2.1 基于主被動(dòng)遙感的漓江水體提取

由于NDWI水體指數(shù)對植被與水體有較好的區(qū)分能力，MNDWI水體指數(shù)適用于增強(qiáng)和提取背景為建成區(qū)的水域信息，EWI水體指數(shù)則結(jié)合了以上兩種指數(shù)的優(yōu)點(diǎn)，較好地對植被、土壤、建筑區(qū)進(jìn)行反射抑制，增強(qiáng)了水體亮度，NDPI水體指數(shù)不僅可以提取小型池塘和低于0.01 hm2的水體，還可以區(qū)分水體與周圍環(huán)境中的植被信息[18-21]. 因此，本文分別基于NDWI、MNDWI、EWI、NDPI指數(shù)和后向散射系數(shù)(S)數(shù)據(jù)進(jìn)行灰度分割，確定佳閾值分割提取漓江水體信息. 具體計(jì)算公式為：

歸一化水體指數(shù)(NDWI)：

(1)

改進(jìn)型歸一化水體指數(shù)(MNDWI)：

(2)

增強(qiáng)型水體指數(shù)(EWI)：

(3)

歸一化差值池指數(shù)(NDPI)：

(4)

以上各式計(jì)算結(jié)果均浮點(diǎn)型數(shù)據(jù)，式中，BGreen為多光譜中的綠波波段，BNIR為近紅外波段，BSWIR表示短波紅外波段.

為了提高水體提取精度，融合了3種指數(shù)和后向散射系數(shù)S的水體提取優(yōu)勢，本文在對比NDWI、MNDWI、EWI指數(shù)與后向散射系數(shù)S提取漓江水體效果的基礎(chǔ)上，根據(jù)水體提取精度確定加權(quán)系數(shù)，并利用ENVI 5.4軟件中波段運(yùn)算工具進(jìn)行加權(quán)求和運(yùn)算，得出主被動(dòng)遙感加權(quán)水體指數(shù)(JQ)影像，再對其灰度進(jìn)行分割、確定最佳閾值提取漓江流域逐月水體信息，具體計(jì)算公式為：

JQ=a·BNDWI+b·BMNDWI+c·BEWI+d·BS

(5)

表1 主被動(dòng)遙感加權(quán)水體指數(shù)中各組成 部分逐月權(quán)重系數(shù)

式中，a、b、c、d為權(quán)重系數(shù)，依據(jù)試錯(cuò)法進(jìn)行反復(fù)試驗(yàn)所設(shè)定，首先通過閾值分割確定NDWI、MNDWI、EWI與后向散射系數(shù)4種指數(shù)逐月水面提取結(jié)果的最適閾值范圍，再依據(jù)不同月份各類指數(shù)分別對漓江水體的提取效果進(jìn)行權(quán)重設(shè)定，逐月權(quán)重系數(shù)值如表1所示，BNDWI為研究區(qū)NDWI指數(shù)數(shù)據(jù)，BMNDWI為研究區(qū)MNDWI指數(shù)數(shù)據(jù)，BEWI為研究區(qū)的EWI指數(shù)數(shù)據(jù)，BS為研究區(qū)后向散射系數(shù)數(shù)據(jù).

Chl.a是水生植物的重要組成成分，懸浮顆粒物是水體中營養(yǎng)物質(zhì)和污染物的重要載體，兩者皆是反映水體質(zhì)量優(yōu)劣程度的重要參數(shù). 基于神經(jīng)網(wǎng)絡(luò)技術(shù)的C2RCC算法可對多種數(shù)據(jù)進(jìn)行大氣校正，該算法可輸出校正后的遙感反射率、固有光學(xué)量、Chl.a濃度、TSM濃度和黃色物質(zhì)等多種產(chǎn)品數(shù)據(jù)，本文則基于此算法反演漓江水體Chl.a及TSM濃度[22-23]；Sentinel-2A中MSI傳感器有3個(gè)波段可獲取植被在近紅外“紅邊”區(qū)間的光譜特征信息，對監(jiān)測植被長勢與健康、水體Chl.a濃度、浮游植物分布等水生態(tài)信息靈敏有效[24]，運(yùn)用紅邊波段的MCI與DoubleR兩種模型可較好地檢測內(nèi)陸、沿海和海洋水域的表面水華和近地表植被信息[25-26]；而可見光區(qū)間內(nèi)，綠光中心波長處Chl.a的光譜反射率，高于相鄰兩側(cè)藍(lán)光、紅光波長的反射率，而形成典型的Chl.a反射峰，該反射峰強(qiáng)度(ρchl)與Chl.a濃度呈正相關(guān)[27-29]. 因此，本文基于Sentinel-2A影像獲取Chl.a、TSM、MCI、DoubleR及ρchl5種指數(shù)，以此來探究漓江年內(nèi)水質(zhì)狀況，具體計(jì)算公式為：

(6)

TSM=1.72IOPbp+3.1IOPbw

(7)

MCI=BVRE-BRed-0.53(BVRE-BRed)

(8)

(9)

(10)

式中，IOPapig表示浮游植物色素吸收特性，IOPbp表示海水中一般顆粒物的后向散射，IOPbw表示海水鈣質(zhì)白色粒子的后向散射，BVRE表示光譜中的紅邊波段，BRed表示光譜中的紅波波段，BBlue表示光譜中的藍(lán)波波段.

2.2 基于基本評(píng)價(jià)單元的漓江水文和水質(zhì)信息定量分析

SDI反映岸線的不規(guī)則程度，刻畫水面寬度差異，以此來研究漓江水體信息演變特征；WDr可定量表征基本評(píng)價(jià)單元內(nèi)水域變化區(qū)域差異；土地利用動(dòng)態(tài)度則適用于分析一定時(shí)間范圍內(nèi)某種土地利用類型的數(shù)量變化情況[30-32]，本文將土地利用動(dòng)態(tài)度用于描述漓江水體動(dòng)態(tài)變化分析，即計(jì)算水體信息變化動(dòng)態(tài)度(K)，研究汛期、枯水期內(nèi)漓江水文與水質(zhì)信息動(dòng)態(tài)變化程度. 因此，本文在278個(gè)基本評(píng)價(jià)單元內(nèi)，采用SDI、WDr和K3個(gè)指標(biāo)來定量分析漓江水體水文和水質(zhì)的年內(nèi)演變特征，具體計(jì)算公式為：

(11)

(12)

(13)

式中，S為水面周長；A為水域面積；Ua、Ub分別為基期、現(xiàn)期的數(shù)值；T為研究時(shí)段長，當(dāng)T的時(shí)段設(shè)定為月時(shí)，K值就是該研究區(qū)某種水文信息的月變化率.

均方根誤差(RMSE)[33]可反映觀測值的可信度、精確度，本文利用均方根誤差進(jìn)行精度驗(yàn)證，計(jì)算公式為：

(14)

式中，xi表示觀測值，yi表示真值，n表示觀測次數(shù).

3 研究結(jié)果

3.1 水體提取、葉綠素a和懸浮物質(zhì)濃度的反演精度驗(yàn)證

為了探究NDWI、MNDWI、EWI、JQ、NDPI與后向散射系數(shù)S6種遙感指數(shù)提取漓江水體效果，本文選取漓江中、下游河段進(jìn)行對比分析. 由圖3可知，MNDWI指數(shù)能較好地提取出處于漓江中游的桂林市城區(qū)漓江水體，但在下游兩岸山峰較多的區(qū)域，山體陰影和水體的區(qū)分不明顯，MNDWI指數(shù)會(huì)把部分山體陰影錯(cuò)誤歸類為漓江水體，導(dǎo)致水域面積偏大. 相對于MNDWI指數(shù)而言，NDWI指數(shù)能較好劃分水體和山體陰影，但在桂林城區(qū)段的水體提取效果不佳，臨近建筑與水體難以通過閾值范圍進(jìn)行區(qū)分.EWI指數(shù)在山區(qū)的提取效果優(yōu)于MNDWI指數(shù)，水體和山體陰影的區(qū)分比較明顯，對城區(qū)的提取能力比NDWI指數(shù)強(qiáng)，能提取出較多的水體部分. 后向散射系數(shù)S完整提取整條漓江的能力比較強(qiáng)，可提取出河流水面寬度較小區(qū)域，但在下游會(huì)把部分山體陰影也劃入漓江，整個(gè)范圍內(nèi)錯(cuò)誤劃分閾值的小圖斑較多. 與其他指數(shù)相比，JQ與NDPI指數(shù)較好地區(qū)分出了桂林城區(qū)的水體，下游被錯(cuò)誤地劃分為水體的山體陰影面積均有所降低，細(xì)小河段以及河流連貫性也有所提高.

圖3 基于Landsat 8、GF-1、Sentinel-1A與Sentinel-2A的6種指數(shù)最佳閾值提取漓江水體結(jié)果對比分析Fig.3 Comparison and analysis of the results of Lijiang River water extraction based on six index optimal thresholds of Landsat 8, GF-1, Sentinel-1A and Sentinel-2A

為驗(yàn)證主被動(dòng)遙感加權(quán)指數(shù)JQ與NDPI兩種水體指數(shù)的提取精度，選用汛期4月和8月與枯水期1月和12月水面提取結(jié)果進(jìn)行對比分析，由圖4中JQ指數(shù)與NDPI指數(shù)提取的水面對比結(jié)果可知，基于NDPI指數(shù)提取的漓江水體有較好的完整性、連貫性；表2為4個(gè)月水面寬度對比結(jié)果，其中JQ指數(shù)的水面寬度選取與實(shí)測水面寬度選取方式一致，NDPI指數(shù)的水面寬度則基于格網(wǎng)取均值；表3中NDPI指數(shù)提取水面寬度的RMSE在2.61～7.04 m之間，而JQ指數(shù)提取水面寬度的RMSE在6.05～9.79 m之間，對比兩種數(shù)據(jù)提取誤差范圍，NDPI指數(shù)提取的數(shù)據(jù)整體誤差明顯小于JQ指數(shù). 以上對比說明基于Sentinel-2A的NDPI指數(shù)提取水體信息有更高的精確度和可信度，下文將基于NDPI指數(shù)提取的水文信息分析漓江水體年內(nèi)變化特征.

圖4 JQ與NDPI指數(shù)水面提取結(jié)果對比Fig.4 Comparison of water surface extraction results of JQ and NDPI indexes

表2 JQ與NDPI指數(shù)提取水面寬度精度對比

表3 JQ與NDPI指數(shù)提取水面寬度的RMSE數(shù)值差異

利用36個(gè)Chl.a與TSM濃度站點(diǎn)監(jiān)測數(shù)據(jù)，通過分析變化趨勢、計(jì)算相關(guān)性來對比驗(yàn)證基于C2RCC算法的Chl.a、MCI、DoubleR及ρchl4種反演模型精度，并在每月上、中、下游反演數(shù)據(jù)中隨機(jī)選取10個(gè)樣點(diǎn)進(jìn)行RMSE計(jì)算. 圖5a～c中基于C2RCC算法的反演數(shù)據(jù)與實(shí)測Chl.a濃度變化趨勢基本一致，使用SPSS統(tǒng)計(jì)分析軟件進(jìn)行相關(guān)性計(jì)算，在95%置信區(qū)間內(nèi)R2為0.980，因每月河流徑流量和水體流動(dòng)速度不同，上、中、下游各部分Chl.a與TSM濃度不均，且河流兩岸旅游開發(fā)程度與人類活動(dòng)影響程度不同，所以通過隨機(jī)選取的樣點(diǎn)計(jì)算出的每月RMSE差異較大，結(jié)合表4可知1-12月Chl.a濃度的RMSE最小為0.18 mg/m3，最大為7.88 mg/m3；而圖5d～f中MCI、DoubleR、ρchl與實(shí)測數(shù)據(jù)的歸一化趨勢可以看出，DoubleR的變化趨勢吻合度較好，R2為0.906，MCI次之，R2為0.883，ρchl則不適于反映漓江Chl.a濃度變化，R2為-0.110；圖5g～i中TSM濃度實(shí)測數(shù)據(jù)與反演數(shù)據(jù)變化趨勢吻合度較好，R2為0.993，但上游12月與中、下游2月誤差較大，1-12月中RMSE最小為0.17 g/m3，最大為12.55 g/m3，下文將基于C2RCC算法的反演數(shù)據(jù)分析漓江Chl.a與TSM濃度的年內(nèi)動(dòng)態(tài)變化規(guī)律.

圖5 Chl.a、TSM濃度的遙感反演數(shù)據(jù)與實(shí)測數(shù)據(jù)的精度驗(yàn)證與趨勢分析Fig.5 Accuracy verification and trend analysis of remote sensing inversion data and measured data of Chl.a and TSM concentrations

表4 Chl.a與TSM 濃度的RMSE差異

3.2 漓江流域水文與水質(zhì)信息年內(nèi)變化定量分析

利用從278個(gè)基本評(píng)價(jià)單元統(tǒng)計(jì)出的河流水面寬度、水域面積、Chl.a與TSM濃度，計(jì)算每月均值與標(biāo)準(zhǔn)差，并利用直方圖(圖6)進(jìn)行趨勢分析. 其中，均值表示一組數(shù)據(jù)的平均水平，反映數(shù)據(jù)集的集中趨勢；標(biāo)準(zhǔn)差反映了數(shù)據(jù)集的離散程度；誤差棒反應(yīng)所測數(shù)據(jù)不確定度的大小.

由圖6a和b可看出漓江1-12月水面寬度與水域面積的整體變化趨勢基本相同，1-8月呈波動(dòng)增長，8-12月呈波動(dòng)型下降. 其中，3-8月上升趨勢較明顯，上升幅度也較大，3月開始進(jìn)入雨季，降雨量增加，水面寬度、水域面積分別由3月的112.56 m、0.057 km2上升至8月的175.32 m、0.088 km2，最大值約為最小值的1.56倍；漓江降水主要集中在夏季，因此5-8月水面寬度、水域面積略高于其他月份. 圖6c與d中Chl.a與TSM濃度變化趨勢基本一致，2月富營養(yǎng)化程度最深，4-11月程度較低、水質(zhì)較好，Chl.a與TSM濃度最大值分別為最小值的38與37倍；兩者變化趨勢與水面寬度、水域面積變化趨勢有一定的負(fù)相關(guān)，即水面寬度、水域面積大的月份，水生植物、藻類等物質(zhì)少，富營養(yǎng)化程度低，水質(zhì)好. 為詳細(xì)觀察漓江水體Chl.a濃度的季節(jié)變化情況，選用離散度較低的1、4、7和10月進(jìn)行分析，由圖7可知每月Chl.a濃度最大值皆位于下游興坪鎮(zhèn)地區(qū)，該地區(qū)為陽朔縣旅游重鎮(zhèn)，旅游開發(fā)區(qū)較多，漓江水體富營養(yǎng)化最嚴(yán)重，4個(gè)月中Chl.a濃度最大值在7月，為39.27 mg/m3，而上、中游水體污染程度較高處均位于興安縣、靈川縣與桂林市城鎮(zhèn)居民區(qū)，其中1月上游水體污染區(qū)域較多，上游地區(qū)興安縣與靈川縣居民區(qū)聚集、人口密度大，農(nóng)田、廠礦分布集中，生活污水、工業(yè)廢水等使得漓江水質(zhì)變差.

圖6 漓江1-12月水文與水質(zhì)信息變化趨勢Fig.6 Change trend of hydrological and water quality information of Lijiang River from January to December

圖7 漓江水體1、4、7和10月各個(gè)河段Chl.a濃度對比分析Fig.7 Comparative analysis of Chl.a concentration in Lijiang River in January, April, July and October

為進(jìn)一步探究漓江1-12月水體變化規(guī)律，本文將每月數(shù)據(jù)分為上、中和下游3部分，表5為水面寬度、水域面積、Chl.a與TSM濃度的全年上、中和下游均值與占比情況，圖8描述了漓江逐月上、中、下游水體信息的動(dòng)態(tài)變化趨勢.

表5 漓江1-12月水文與水質(zhì)信息上、中、下游均值與占比

根據(jù)表5可知，上、中、下游水面寬度、水域面積與Chl.a濃度變化較小，TSM濃度變化最大，且為中游>下游>上游. 結(jié)合圖8可知，①上、中、下游水面寬度變化趨勢相似，但存在差異. 1-3月水面寬度呈下降趨勢，且2-3月下降顯著，該期間氣溫回升但降水少，水面寬度明顯偏低；3-8月雨量增多，水面寬度變大，且中、下游3-5月上升幅度最大，由114 m上升至178 m；8-12月水面寬度雖有波動(dòng)但下降趨勢明顯，上、中、下游降幅分別為32、54和38 m，中游降幅最為顯著；上、中游水面寬度在8月達(dá)到最大值，下游最大值出現(xiàn)在5月，上、中、下游最小值均出現(xiàn)在3月(圖8a)，說明3月漓江流域徑流量最小. ②水域面積與水面寬度總體變化趨勢基本一致，說明水域面積與水面寬度有一定正相關(guān)性，即降水較多的月份，水面寬度增大，水域面積也增大(圖8b). ③中、下游1-12月Chl.a濃度變化趨勢相同，1-3月變化幅度較大，2月達(dá)全年最大值，說明該月份中、下游水體中水生植物最多，富營養(yǎng)化程度嚴(yán)重，而上游1-3月Chl.a濃度呈下降趨勢；漓江水體4-10月Chl.a濃度均低于1.5 mg/m3，該時(shí)期降水較多，水體流動(dòng)性大，水生植物的生長受到抑制，水質(zhì)較好，而10-12月降水少，水生植物覆蓋率隨之上升，Chl.a濃度略有增加(圖8c). ④TSM濃度變化趨勢與Chl.a濃度基本相同，但上游1-3月變化趨勢稍有差異，2月TSM濃度較低，該期間上游降水較多，水體流動(dòng)性大，懸浮物覆蓋率降低，該趨勢與中、下游變化趨勢相反(圖8d). 由于1月數(shù)據(jù)離散度小，且基于基本評(píng)價(jià)單元提取的數(shù)據(jù)較完整，因此選取1月漓江流域278個(gè)基本評(píng)價(jià)單元的數(shù)據(jù)分析水文與水質(zhì)信息的變化情況(圖8e～h)，其中，水面寬度與水域面積變化趨勢一致，上、中游水面寬度與水域面積略大于下游，兩者下游變化趨勢皆呈上升狀態(tài)，但整條河流水面寬度與水域面積變化幅度較大，最大值分別為最小值的15.8和7.6倍；Chl.a與TSM濃度整體呈下降趨勢，上游下降幅度最為顯著，Chl.a濃度由第1個(gè)基本評(píng)價(jià)單元的11.99 mg/m3下降至第71個(gè)基本評(píng)價(jià)單元的2.63 mg/m3，TSM濃度由第1個(gè)基本評(píng)價(jià)單元的24.32 g/m3下降至第71個(gè)基本評(píng)價(jià)單元的0.59 g/m3，上游地區(qū)村莊較多，居民環(huán)保意識(shí)不足，致使河流污染嚴(yán)重；Chl.a與TSM濃度在中游第115個(gè)基本評(píng)價(jià)單元處有一個(gè)峰值，該區(qū)域的水域面積與水面寬度較小，表明1月桂林市該區(qū)域降水少、水生植物與懸浮物質(zhì)較多.

圖8 漓江水文與水質(zhì)信息變化趨勢Fig.8 Change trend of hydrological and water quality information of Lijiang River

3.3 漓江流域水文與水質(zhì)信息變化差異的定量分析

水面變化區(qū)域差異值可以比較直接地觀察河流水面變化的情況，WDr值越大說明變化程度越大. 以全年數(shù)據(jù)均值作為Ua，每月均值作為Ub，計(jì)算河流信息變化差異值，并以折線圖來刻畫每個(gè)月份的變化情況.

由圖9可以看出，漓江水面寬度、水域面積、Chl.a和TSM濃度變化與圖8中變化趨勢呈正相關(guān). 圖9a和b中水面寬度WDr與水域面積WDr變化趨勢相同，4條差異變化折線在8-10月都是正值，即為正向變化，3、4和12月都為負(fù)向變化，3月變化差異最大，水面寬度與水域面積流域WDr分別為-0.24與-0.23. 漓江上游段在3-7、12月為負(fù)值，即相對全年來說，該期間降水較少，水面寬度和水域面積處于減小狀態(tài)，而1-2、8-11月降水較多. 其中，3月干旱嚴(yán)重，負(fù)向變化差異最大，水面寬度差異為-0.28，水域面積差異為-0.30. 漓江水面寬度與水域面積的流域變化差異與中、下游部分的變化趨勢大體一致，且中、下游5-10月雨量大，水面寬度、水域面積較大. 中、下游正向變化差異最大的月份為5和8月，負(fù)向變化最大的仍為3月. 漓江3-12月各河段Chl.a濃度變化趨勢相同，4-11月WDr皆為負(fù)值，說明該時(shí)期Chl.a低于全年平均濃度，而1-3月上游Chl.a濃度變化趨勢與中、下游部分稍有差異，該時(shí)期上游WDr值雖呈下降趨勢，但Chl.a均高于全年濃度(圖9c)；漓江中、下游Chl.a濃度變化與整條流域變化趨勢相同，2月正向變化差異最大，總流域WDr為4.16，即與全年相比2月Chl.a濃度增幅最大，同時(shí)說明該月份水生植物覆蓋率最高. 漓江TSM與Chl.a濃度變化趨勢大體一致，兩者變化趨勢與水域面積、水面寬度的變化趨勢呈一定的負(fù)相關(guān)；上游1-3月TSM濃度與中、下游變化趨勢相反，2月TSM濃度增幅最小(圖9d). 1月漓江水面寬度WDr與水域面積WDr整體變化趨勢較穩(wěn)定，但仍有部分區(qū)域變化幅度較大；水面寬度與水域面積下游區(qū)域WDr負(fù)值較多，兩者負(fù)值基本評(píng)價(jià)單元數(shù)量分別占總負(fù)值數(shù)量的45.14%與46.56%，即與全年相比，1月下游地區(qū)降水較少，水面寬度、水域面積多處于減小狀態(tài)(圖9e，f). 而漓江水體Chl.a的WDr與TSM的WDr皆呈下降趨勢，其中上游Chl.a的WDr均為正值，TSM的WDr僅有第26個(gè)基本評(píng)價(jià)單元為負(fù)值，上游興安縣與靈川縣內(nèi)漓江水質(zhì)較差；中游42%格網(wǎng)Chl.a的WDr為負(fù)值，85%格網(wǎng)TSM的WDr為負(fù)值，其中第113個(gè)格網(wǎng)正向變化最大，WDr達(dá)到了1.96；下游81.48%的格網(wǎng)Chl.a的WDr為負(fù)值，TSM的WDr均為負(fù)值(圖9g，h)，說明1月漓江上游水面水生植物與藻類等懸浮物質(zhì)覆蓋率最大，富營養(yǎng)化程度最深，下游水體水面寬度與水域面積雖呈減少趨勢，但水質(zhì)較好.

圖9 漓江水文與水質(zhì)信息變化區(qū)域差異Fig.9 Regional differences of hydrological and water quality information changes of Lijiang River

3.4 漓江流域河岸線發(fā)育系數(shù)定量分析

河岸線發(fā)育系數(shù)可表征岸線的不規(guī)則程度，也可表征水面變化程度. 漓江河流SDI與水面周長變化呈正相關(guān)，而與水域面積變化呈負(fù)相關(guān)(圖10). 由圖10a和b可知，1-12月水面周長與SDI皆有下降趨勢，但岸線發(fā)育系數(shù)變化幅度較大，可分成3個(gè)階段：小幅下降(1-3月)、急劇下降(3-5月)和穩(wěn)定(5-12月). 1-3月雨量較少，流域水量低，水域面積小，河岸線較為曲折，河流水面周長較大，但仍有下降趨勢. 3-5月降水增多，河流面積迅速增加，該期間河岸線變的平緩，河流水面周長降低. 5-12月SDI值在1.6～1.7之間波動(dòng)，岸線發(fā)育率較低，該階段河流水面周長雖有下降趨勢，但由于雨期的來臨，水量增多，水域面積明顯大于1-4月，岸線不規(guī)則程度也較低. 圖10c表征了漓江1月278個(gè)基本評(píng)價(jià)單元的SDI變化情況，整體來看該月份岸線發(fā)育率呈下降趨勢，但中游起伏較大，其中，第106個(gè)基本評(píng)價(jià)單元的SDI值達(dá)到5.89，桂林市內(nèi)該水域面積數(shù)值較小，水面周長達(dá)到最大值，說明此基本評(píng)價(jià)單元中的水量少，岸線較為曲折.

圖10 漓江河流水面周長和岸線發(fā)育系數(shù)的變化趨勢Fig.10 Variation trend of river circumference and coastline development coefficient of Lijiang River

表6 基于汛期和枯水期的水文與 水質(zhì)信息動(dòng)態(tài)度變化

3.5 漓江流域水文與水質(zhì)信息汛期、枯水期變化定量分析

為探究漓江汛期和枯水期水體信息動(dòng)態(tài)變化規(guī)律，本文以枯水期均值作為動(dòng)態(tài)度公式中的Ua，汛期均值作為Ub，來計(jì)算汛期相對于枯水期的變化情況. 由表6可知，水面寬度與水域面積的總體流域動(dòng)態(tài)度都是正值，Chl.a與TSM的總河流動(dòng)態(tài)度都是負(fù)值，即汛期時(shí)降雨多，漓江水面寬度、水域面積都大于枯水期，而水體中的Chl.a與TSM濃度低于枯水期，其中Chl.a與TSM濃度的變化較大，動(dòng)態(tài)度均小于-7%，下游水面寬度與水域面積動(dòng)態(tài)度分別為2.24%與2.41%，說明水中植物與TSM濃度變化最大；上游水面寬度與水域面積為負(fù)值，即上游汛期水面寬度與水域面積小于枯水期，說明上游汛期降水少于枯水期；而中、下游降水多于枯水期，且下游降水最多、動(dòng)態(tài)變化最大；汛期時(shí)漓江水質(zhì)好，富營養(yǎng)化程度低，整條河流Chl.a與TSM濃度均小于枯水期，其中下游變化最大，動(dòng)態(tài)度為-13.09%.

本文基于278個(gè)基本評(píng)價(jià)單元進(jìn)一步分析各河段漓江汛期和枯水期水體信息動(dòng)態(tài)變化，具體見圖11. 由圖11a可知，汛期大部分地區(qū)的水面寬度在100～250 m范圍內(nèi)波動(dòng)，但上游部分地區(qū)變化趨勢較大，其中最大值為靈川縣第31個(gè)格網(wǎng)，水面寬度達(dá)到了338.31 m，最小值為第5個(gè)格網(wǎng)位于興安縣境內(nèi)，數(shù)值為40.97 m. 枯水期上游水面寬度波動(dòng)較大，中游水面寬度有下降趨勢，下游略呈上升狀態(tài)(圖11b). 汛期與枯水期水域面積與水面寬度變化趨勢基本一致(圖11c，d)，但中游第105個(gè)基本評(píng)價(jià)單元水域面積較大，汛期與枯水期水域面積分別為0.14和0.17 km2，該地區(qū)全年降水多、水量多、水面寬. 汛期與枯水期Chl.a濃度整體呈“雙谷型”變化(圖11e，f)，中游第100～110個(gè)基本評(píng)價(jià)單元處濃度較高，該地區(qū)村莊聚集，較多生活污水的排放導(dǎo)致河流富營養(yǎng)化程度加深. 除第26個(gè)基本評(píng)價(jià)單元TSM濃度為16.85 g/m3，達(dá)整條河流的最大值外，汛期TSM濃度均在0～3 g/m3之間(圖11g)，該水域的水面寬度、水域面積與Chl.a濃度也較大，說明汛期該地區(qū)水量多，但水質(zhì)差. 枯水期TSM與Chl.a濃度呈正相關(guān)(圖11h)，上游興安縣第1～10個(gè)、中游桂林市區(qū)第10～110個(gè)基本評(píng)價(jià)單元所在水域水面富營養(yǎng)化程度高，其中第107個(gè)基本評(píng)價(jià)單元處TSM濃度為15.03 g/m3，達(dá)整條河流的最大值，該水域的水面寬度、水域面積與Chl.a濃度也較大，說明該處水體污染程度最嚴(yán)重.

圖11 漓江汛期和枯水期水文與水質(zhì)信息變化趨勢Fig.11 Change trend of hydrological and water quality information in flood and dry seasons of Lijiang River

4 結(jié)論與討論

本文以2016-2020年Landsat 8 OLI、GF-1和Sentinel-2A多光譜數(shù)據(jù)和2018年逐月Sentinel-1A IW模式的SAR影像為主要數(shù)據(jù)源，利用構(gòu)建的多種主被動(dòng)遙感水體指數(shù)提取了漓江逐月的水域信息，選用4種模型反演了漓江逐月的Chl.a與TSM濃度，并基于基本評(píng)價(jià)單元采用多種水體指標(biāo)定量分析漓江水文、水質(zhì)的年內(nèi)、枯水期和汛期上、中和下游的時(shí)空演化特征，結(jié)果表明：(1)主被動(dòng)遙感加權(quán)指數(shù)JQ與NDPI指數(shù)的提取效果優(yōu)于NDWI、MNDWI、EWI與后向散射系數(shù)，而NDPI指數(shù)與JQ指數(shù)相比提取精度更高、可信度更強(qiáng)，水面寬度RMSE在2.61～7.04 m之間. (2)基于C2RCC算法反演的數(shù)據(jù)可較好地揭示漓江Chl.a與TSM濃度的年內(nèi)動(dòng)態(tài)變化規(guī)律，Chl.a濃度的RMSE處于0.18～7.88 mg/m3之間，R2為0.980，TSM濃度的RMSE為0.17～12.55 g/m3，R2為0.993，而DoubleR可較好地反映Chl.a濃度的變化趨勢，R2為0.906，MCI的R2為0.883，ρchl的R2為-0.110，則不適用于反映水體Chl.a濃度變化. (3)1月水面寬度與水域面積變化較穩(wěn)定，大部分基本評(píng)價(jià)單元所在水域的水面寬度和水域面積分別處于60～240 m與0.03～0.12 km2區(qū)間內(nèi)，而Chl.a與TSM濃度明顯呈下降趨勢，其中上游降幅最大，由第1個(gè)基本評(píng)價(jià)單元至第71個(gè)基本評(píng)價(jià)單元Chl.a與TSM濃度降幅分別為9.36 mg/m3與23.73 g/m3，中、下游變化則較穩(wěn)定，實(shí)測數(shù)據(jù)則依靠站點(diǎn)監(jiān)測，所得數(shù)據(jù)較分散，無法進(jìn)行連續(xù)性分析. (4)漓江1-12月平均水面寬度在100～175 m范圍內(nèi)，水域面積處于0.05～0.09 km2之間，其中 5-10月降水較多，水面寬度在150 m以上，水域面積在0.07 km2以上，水面寬度和水域面積最大值出現(xiàn)在8月；而5--10月水體富營養(yǎng)化程度低，水質(zhì)較好，Chl.a濃度低于0.6 mg/m3，TSM濃度低于1.3 g/m3，漓江2月水質(zhì)最差，Chl.a濃度為9.82 mg/m3，TSM濃度為15.32 g/m3，水體富營養(yǎng)化程度較高地區(qū)主要集中于上、中游的興安縣、靈川縣等城鎮(zhèn)居民區(qū)，以及下游旅游開發(fā)較多的興坪鎮(zhèn). ②漓江水體水面寬度WDr與水域面積WDr在8-10月都大于0，即為正向變化，3、4、12月全為負(fù)向變化，其中3月負(fù)向變化最大，水面寬度與水域面積的WDr分別為-0.24與-0.23，該月份降水最少；漓江流域4-11月Chl.a的WDr均小于0，該時(shí)期Chl.a濃度低于全年平均值，TSM濃度的WDr值也在0附近波動(dòng)，該時(shí)期漓江水質(zhì)較好，富營養(yǎng)化程度低. 漓江1-12月水面周長與SDI呈下降趨勢，4-5月降幅較大，其中上游降幅為0.57，下游降幅為0.27，5-12月SDI則基本穩(wěn)定在1.6～1.8之間，該期間水面變化程度較小. 漓江水體水面寬度、水域面積與Chl.a濃度在上、中、下游總體變化較小，TSM濃度變化最大，且為中游>下游>上游；汛期時(shí)降水多、水體流動(dòng)性強(qiáng)，大部分地區(qū)平均水面寬度在100～250 m范圍內(nèi)，但上游地區(qū)波動(dòng)較大，TSM濃度則基本處于0～3 g/m3之間，水體富營養(yǎng)化程低，水質(zhì)較好.

本文利用Landsat 8 OLI、GF-1、Sentinel-1A與Sentinel-2A 4種衛(wèi)星影像提取漓江水面數(shù)據(jù)，并選用C2RCC算法、MCI、DoubleR與ρchl進(jìn)行漓江水質(zhì)反演，研究結(jié)果與漓江實(shí)測數(shù)據(jù)較為吻合，但還存在一些不足，將會(huì)在以后工作研究中逐一完善：1)因部分月份中影像受云霧和陰雨天氣影響，部分星地?cái)?shù)據(jù)存在不同步現(xiàn)象. 2)在C2RCC處理過程中，所有參數(shù)均依據(jù)美國宇航局生物光學(xué)海洋算法數(shù)據(jù)集的訓(xùn)練結(jié)果進(jìn)行設(shè)置，但為進(jìn)一步提高實(shí)驗(yàn)結(jié)果的準(zhǔn)確性，應(yīng)盡可能根據(jù)全球各地不同近岸水體特點(diǎn)，通過當(dāng)?shù)厮|(zhì)實(shí)測數(shù)據(jù)來優(yōu)化調(diào)整為更有區(qū)域適應(yīng)性的水質(zhì)參數(shù). 3)主被動(dòng)遙感加權(quán)指數(shù)JQ中各指數(shù)的權(quán)重是通過試錯(cuò)法反復(fù)試驗(yàn)所設(shè)定的，存在主觀賦值所帶來的結(jié)果不確定性等問題，同時(shí)還應(yīng)當(dāng)考慮不同河段地物類型的特點(diǎn)，進(jìn)行分段差別化處理.

5 附錄

附表Ⅰ見電子版(DOI: 10.18307/2021.0306).