吳 杰,董小濤,張 珂,3,4,李 曦,吳 南

(1.河海大學(xué)水文水資源學(xué)院,江蘇 南京 210098; 2.水利部綜合事業(yè)局,北京 100032;3.長江保護(hù)與綠色發(fā)展研究院,江蘇 南京 210098; 4.中國氣象局水文氣象重點開放實驗室,江蘇 南京 210098)

消落帶作為一種特殊的濕地生態(tài)系統(tǒng),指由于季節(jié)或人為控制因素導(dǎo)致水位消漲而在水庫周邊形成的陸地和水生生態(tài)系統(tǒng)功能相互作用的巨大隔離帶,是庫區(qū)生態(tài)保護(hù)的重要對象[1-3]。水位波動是對湖泊生態(tài)系統(tǒng)影響最為重要的因素[4],也是形成消落帶的主要原因。與自然形成的河岸帶不同,水庫消落帶的水位變化主要受水庫實際調(diào)度影響,水庫調(diào)蓄造成的庫區(qū)水位周期性波動變化會對消落帶的生態(tài)結(jié)構(gòu)造成重大影響。目前對消落區(qū)的研究主要集中在土地利用變化[5-6]、消落帶分區(qū)等方面,對消落帶時空變化特征分析的研究較少。

淹沒頻率作為一項重要的水文特征,反映了空間位置上水體淹沒時間的長短及次數(shù),是影響濕地生態(tài)系統(tǒng)的重要因素[7],同時會影響消落帶的生態(tài)結(jié)構(gòu)和生態(tài)類型劃分。張麗麗等[8]通過計算鄱陽湖自然保護(hù)區(qū)濕地的水深和淹水頻率,分析了植被群落對這2種水文特征的耐受性和敏感性。許超等[9]基于MODIS遙感數(shù)據(jù)計算淹沒頻率,分析其變化;谷娟等[10]基于混合像元分解模型分析了鄱陽湖的淹沒頻率變化,探究時間序列上濕地植被特征的空間響應(yīng)。水位的季節(jié)性波動變化導(dǎo)致水庫周邊形成不同的淹沒頻率,濕地中的植被也受到不同程度的水淹影響,對水位變化不適應(yīng)的植被大量消亡,造成該濕地區(qū)域植被群落減少,環(huán)境的穩(wěn)定性受到破壞,水土流失、水質(zhì)惡化等損害生態(tài)系統(tǒng)的問題層出不窮。因此,綜合分析消落帶淹水范圍及淹水頻率的時空變化特征,加強(qiáng)對消落帶水體分布的高精度、高頻次監(jiān)測具有重要意義。

遙感信息技術(shù)作為一種可以高效獲取、分析及處理空間地理信息的技術(shù)手段,已廣泛應(yīng)用于水文領(lǐng)域的趨勢變化分析及宏觀監(jiān)測[11-15]等領(lǐng)域。高分一號衛(wèi)星由我國自主研發(fā)并于2013年4月發(fā)射成功,有效載荷為2臺2m分辨率全色、1臺8m分辨率多光譜相機(jī)(PMS)和4臺16m分辨率多光譜寬幅相機(jī)(WFV),具有高空間分辨率、多相機(jī)等特點,滿足了我國對高時空分辨率遙感數(shù)據(jù)的應(yīng)用需求[16]。目前常見的遙感影像水體信息提取方法包括單波段閾值法[17]、多波段譜間關(guān)系法[18]、水體指數(shù)法及分類法[19- 20]等。

為研究消落帶生態(tài)類型劃分,本文提出一種線性插值模擬計算水體淹沒頻率的方法,并結(jié)合遙感影像水體信息提取[21-23]和幾何校正[24-25]技術(shù),分析了三河口水庫蓄水階段的水體時空變化特征,并根據(jù)模擬水體的精度反映淹沒頻率結(jié)果的可靠性,以期提高淹沒頻率計算的可靠性。

1 研究方法

1.1 遙感數(shù)據(jù)處理及精細(xì)配準(zhǔn)

本文采用的遙感數(shù)據(jù)為中國資源衛(wèi)星網(wǎng)站發(fā)布的高分一號衛(wèi)星WFV16 m分辨率遙感影像。在進(jìn)行后續(xù)數(shù)據(jù)分析及應(yīng)用前,利用ENVI二次開發(fā)工具IDL批量預(yù)處理不同日期遙感影像,步驟包括輻射定標(biāo)、大氣校正[26]、正射校正等。為了消除因傳感器自身因素、大氣條件、太陽位置和角度等引起的遙感影像空間位置誤差,選取固定的地面建筑物作為參照物進(jìn)行精細(xì)配準(zhǔn),使配準(zhǔn)后的影像空間位置差異在誤差允許范圍內(nèi)[24]。根據(jù)研究區(qū)矢量文件對預(yù)處理后的遙感影像進(jìn)行裁剪得到用于水體提取的源數(shù)據(jù)。

1.2 最大類間方差迭代水體提取法

最大類間方差法又稱為大津法(Otsu),由學(xué)者Otsu[27]于1979年提出,是一種簡單可行、應(yīng)用廣泛的圖像分割閾值算法。該算法可以找到某一灰度值,即閾值,將圖像分為前景和背景。本文根據(jù)歸一化差分水體指數(shù)(normalized difference water index,NDWI)將圖像分為水體和非水體。閾值為使前景和背景的類間方差最大的灰度值,類間方差越大說明前景與背景灰度值的離散程度越大,同時也說明類內(nèi)方差越小,即前景和背景內(nèi)部灰度值離散程度越小,故此時的閾值選取也更加合理。

基于改進(jìn)的最大類間方差迭代水體提取法[11],通過迭代不同膨脹算子確定最優(yōu)提取結(jié)果。首先通過經(jīng)驗閾值和邊緣檢測法初步篩選水體信息并利用連通域標(biāo)記法去除其他無關(guān)水體,然后迭代多個膨脹算子并基于形態(tài)學(xué)膨脹法建立庫區(qū)水體的緩沖區(qū),再對緩沖區(qū)內(nèi)區(qū)域應(yīng)用最大類間方差法,根據(jù)判斷準(zhǔn)則自適應(yīng)地確定水體與非水體的分割閾值。自適應(yīng)最優(yōu)閾值選擇公式為

Smin=min{|Si-Si-1|} (i=1,2,3)

(1)

式中:Smin為最小的水體提取面積差值;Si、Si-1分別為第i次、第i-1次迭代得到的水體面積。

1.3 連通域標(biāo)記法

連通域指的是圖像中由相同并且相鄰的像素值組成的圖像區(qū)域。不同的相鄰像素點形成不同的區(qū)域,這些不同區(qū)域代表了不同的連通區(qū)域,像素的鄰接關(guān)系決定了像素與周邊的連通關(guān)系,通常分為四鄰接和八鄰接關(guān)系[28]。本文圖像為二值圖像,最小單位為像素,采用四鄰接關(guān)系,即考慮像素的上、下、左、右等4個方向的連接狀態(tài),通過遍歷圖像,記下每一行或列中團(tuán)和標(biāo)記的等價對,并通過等價對對原來的圖像進(jìn)行重新標(biāo)記,步驟分為橫向合并、縱向合并和逆向合并。

通過連通域區(qū)域標(biāo)記算法對各連通水體進(jìn)行標(biāo)記,對研究區(qū)內(nèi)不同水體進(jìn)行區(qū)分,選擇最大連通域作為庫區(qū)水體范圍,可以有效去除圖像中存在的其他水體和噪音。

1.4 水庫淹沒頻率計算方法

水庫淹沒頻率是指某個位置上被水淹的時間與總研究時間的比值。受遙感影像時間分辨率、云體干擾、數(shù)據(jù)質(zhì)量等因素干擾,無法獲得日尺度的水體時空變化特征。本文選擇影像在時間尺度上保證相鄰兩次日期高分一號衛(wèi)星數(shù)據(jù)的水庫水體時空變化呈單調(diào)變化趨勢,對缺少數(shù)據(jù)的日期,基于線性插值模擬結(jié)果構(gòu)建淹沒頻率計算公式:

(2)

式中:P為對應(yīng)像元在研究時間內(nèi)的淹沒頻率;Nj為第j+1期影像該位置的水體劃分,水體與非水體分別對應(yīng)為1、0;ΔT為兩期影像間隔時間;T為總研究時間。

2 實例驗證

2.1 工程概況

引漢濟(jì)渭工程為陜西省內(nèi)重大跨流域調(diào)水工程,三河口水利樞紐是引漢濟(jì)渭工程的調(diào)蓄中樞,目前已完成水庫初期蓄水階段,消落帶正在逐漸形成。三河口水利樞紐位于佛坪縣與寧陜縣境交界、漢江一級支流子午河中游峽谷段,其壩址位于大河壩鄉(xiāng)三河口村下游約2km處。三河口水利樞紐主要由攔河大壩、泄洪放空系統(tǒng)、供水系統(tǒng)和連接洞等水工建筑物組成。

三河口水庫回水淹沒范圍涉及陜西省寧陜縣,庫區(qū)上游左岸有椒溪河、蒲河和汶水河等支流匯入。該水庫庫容較小,是整個引漢濟(jì)渭工程中具有較大水量調(diào)節(jié)能力的核心項目,承擔(dān)供水、調(diào)蓄、發(fā)電以及滿足下游綜合用水要求等多項任務(wù)。

根據(jù)三河口水庫的蓄泄運行規(guī)律,水庫水位在正常蓄水位643m至庫正常運行死水位558.0m之間變動。圖1展示了三河口水庫最大淹沒范圍1000m緩沖區(qū)。

圖1 研究區(qū)示意圖Fig.1 Schematic diagram of the study area

2.2 結(jié)果與分析

2.2.1 庫區(qū)水位變化趨勢

選取2020年1月22日到2021年11月8日共25幅預(yù)處理后的高分一號衛(wèi)星遙感影像,應(yīng)用Otsu法選取閾值得到水體提取結(jié)果。自2020年1月22日至2020年2月29日水庫始終處于最低水位,2021年9月8日后達(dá)到最大水位,漲水階段共歷時約18月,預(yù)計消落帶總面積為13.95km2。對遙感影像進(jìn)行目視解譯,得到三河口水庫水體面積的變化階段:2020年2月以前,水庫維持在死水位不變;2020年2月至2021年2月,水庫水位第一次上漲后維持固定水位,漲落周期為11月;2021年2—6月,水庫水位在大幅上漲后維持固定水位,漲落周期為4月;2021年6—10月水庫水位繼續(xù)上漲至正常蓄水位,之后水位維持在正常蓄水位。

圖2為時間序列上Otsu法提取水體的變化趨勢,線性擬合R2達(dá)到0.58468,總體呈現(xiàn)穩(wěn)步上升趨勢。20201129景影像水體提取結(jié)果較相鄰日期波動較大,結(jié)合目視解譯遙感影像發(fā)現(xiàn)將大壩周圍建筑及山體陰影誤分為水體導(dǎo)致水體提取結(jié)果較實際結(jié)果偏大;20210117、20210605景影像水體提取結(jié)果較周圍日期向下波動則是由于受云層影響水體上游末端細(xì)小水體未充分提取,導(dǎo)致水體提取結(jié)果較實際結(jié)果偏小。三河口水庫周邊存在部分建筑和山體陰影,而NDWI指數(shù)對山體陰影、建筑物和水體的區(qū)分能力較弱,存在將其誤分為水體導(dǎo)致部分日期水體面積偏大的情況,呈現(xiàn)為圖2中蓄水過程線的波動。

圖2 三河口庫區(qū)水體面積時間序列Fig.2 Time series of water body area in Sanhekou Reservoir

由圖3可知,水體邊界整體沿水系擴(kuò)展,水體面積增加,水位由558m上漲至643m,落差達(dá)85m,符合水庫實際調(diào)度情況。其中,20210219景影像之前庫區(qū)水體面積變化較小,存在將壩址、山體陰影錯分的現(xiàn)象;20210219景影像后隨著水庫大量蓄水,水位大幅上升,水體面積增加,水體與周邊噪音的區(qū)分程度較好,對水體的誤提取情況明顯減少。

圖3 2020—2022年三河口庫區(qū)水體提取結(jié)果Fig.3 Results of extracted water body area of Sanhekou Reservoir from 2020 to 2022

2.2.2 水庫空間淹沒頻率特征

對提取的2020—2022年三河口水庫水體邊界,在進(jìn)行大氣校正等預(yù)處理手段后,通過精細(xì)幾何校正使影像之間地物的誤差控制在誤差允許范圍內(nèi)。將處理后的影像進(jìn)行疊置分析并按時間序列進(jìn)行插值處理,計算年淹沒頻率。

圖4(a)為得到的空間上16m分辨率每個像元的淹沒頻率,反映了庫區(qū)空間淹沒狀態(tài),淹沒程度沿三河口水庫三條河流上游向庫區(qū)逐漸增加,較真實地反映了庫區(qū)淹沒狀態(tài)。但永久淹沒區(qū)較死水位對應(yīng)的水面面積偏小,分析誤差原因主要為不同遙感影像間空間位置差異。通過人工修正按淹沒頻率(0, 33.3%]、(33.3%, 66.6%]、(66.6%, 100%)、100%將三河口水庫水位變動區(qū)域劃分為4個淹沒等級(圖4(b)):經(jīng)常出露、半淹半露、經(jīng)常淹沒和永久淹沒。可以看出壩址周圍建筑及山體陰影的誤分影響了淹沒頻率的計算。

圖4 三河口水庫1000m緩沖區(qū)內(nèi)的淹沒頻率特征及淹沒頻率分區(qū)Fig.4 Characteristics and zoning of inundation frequency with 1000m buffer zone in Sanhekou Reservoir

從表1可以看出淹沒頻率面積占比隨淹沒程度的增加逐漸遞減,空間上表現(xiàn)為由內(nèi)向外遞減的趨勢。由于研究期限內(nèi)水庫正處于初期蓄水階段,淹沒頻率區(qū)間為(0, 33.3%]的區(qū)域占比最大,面積達(dá)到11.72km2,占淹沒區(qū)面積的81.84%,是淹沒區(qū)的主體部分;(33.3%, 66.6%]、(66.6%, 100%)、100%的區(qū)域面積分別為1.55、0.68、0.37km2,占淹沒區(qū)面積的10.82%、4.75%、2.58%。

表1 緩沖區(qū)淹沒頻率統(tǒng)計

2.2.3 水庫淹沒頻率可靠性分析

選擇存在遙感影像的日期進(jìn)行驗證,基于淹沒頻率空間特征結(jié)果,模擬得到水體分布結(jié)果。淹沒頻率結(jié)果的可靠性通過模擬得到的水體分布結(jié)果的精度來反映。用于驗證的遙感影像需要滿足未參與頻率計算且處于水庫調(diào)蓄階段內(nèi)。

由圖5可知,20210312與20210323、20210708、20210717水體提取結(jié)果基本一致,說明它們處于水庫上漲的同一階段。結(jié)合人工目視解譯實際遙感影像,分別選取100個水體和非水體像元構(gòu)建混淆矩陣計算總體精度和Kappa系數(shù),結(jié)果見表2。

表2 水體模擬結(jié)果精度

圖5 不同日期模擬水體分布Fig.5 Simulated water body distribution

由表2可知, 6次水體模擬結(jié)果的總體精度和Kappa系數(shù)均值分別達(dá)到86.5%和73%,結(jié)果較為滿意,說明由本文方法得到的淹沒頻率特征符合該地區(qū)的水庫調(diào)蓄情況,具有較強(qiáng)的可靠性。

3 結(jié) 語

本文采用線性插值法計算了不同時段的水淹狀態(tài),模擬了蓄水時段的水體變化,結(jié)合高分辨率遙感數(shù)據(jù),得到庫區(qū)淹沒頻率的時空分布特征。根據(jù)淹沒程度不同分成經(jīng)常出露、半淹半露、經(jīng)常淹沒、長久淹沒4種淹沒區(qū),空間上呈現(xiàn)為由內(nèi)向外遞減,面積占比分別為81.84%、10.82%、4.75%、2.58%,在水庫進(jìn)入正常調(diào)蓄階段后,為消落帶形成后的生態(tài)類型劃分提供了借鑒思路。本文選擇未參與計算且存在遙感數(shù)據(jù)的日期對水體模擬結(jié)果進(jìn)行精度驗證,總體精度最大值為90%,最小值為82%,均值為86.5%,Kappa系數(shù)最大值為80%,最小值為64%,均值為73%,說明淹沒頻率時空分布特征可靠性較好。總體上該頻率計算方法簡單有效,在水庫調(diào)蓄情況下較為適用,精細(xì)化了淹沒頻率的計算。但強(qiáng)降雨等自然因素的干擾會影響其適用性;水體提取結(jié)果的準(zhǔn)確性也會影響頻率計算的結(jié)果;在小型水庫由于對幾何校正的要求變高也會導(dǎo)致精度的降低,對此在未來的研究中仍需進(jìn)一步深入探討。