張 珂，吳 南，徐國鑫，劉林鑫1，，范亞洲1，，張企諾，周佳奇1，，劉 挺

(1.河海大學(xué)水文水資源與水利工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210098；2.河海大學(xué)水文水資源學(xué)院，江蘇 南京 210098； 3.長江保護(hù)與綠色發(fā)展研究院，江蘇 南京 210098； 4.中國氣象局-河海大學(xué)水文氣象研究聯(lián)合實(shí)驗(yàn)室，江蘇 南京 210098； 5.陜西省引漢濟(jì)渭工程建設(shè)有限公司，陜西 西安 710010)

隨著空間技術(shù)的發(fā)展，通過遙感影像提取水體在洪澇災(zāi)害、水資源監(jiān)測以及水環(huán)境保護(hù)等領(lǐng)域得到廣泛應(yīng)用[1-3]，水體提取已經(jīng)成為遙感技術(shù)研究領(lǐng)域的重要分支之一[4]。目前常見的水體提取方法有單波段閾值法、譜間關(guān)系法、水體指數(shù)法、決策樹法、密度分割法以及面向?qū)ο蠓ǖ萚5-9]，其中水體指數(shù)法是一種方便高效的水體識(shí)別法，應(yīng)用也比較成熟[10]。由于衛(wèi)星的分辨率不同，水體提取結(jié)果也有所差異，鄧勁松等[11]基于SPOT高分辨率影像提出一種波段間的提取方法，有效地將居民地和水體分開；周成虎等[12]針對(duì)AVHRR影像提出了自動(dòng)提取識(shí)別水體的描述模型，并在太湖、淮河等多個(gè)地區(qū)得到應(yīng)用；王培培[13]利用ETM遙感影像數(shù)據(jù)，通過光譜特征分析和歸一化植被指數(shù)識(shí)別水體，實(shí)現(xiàn)了水體的自動(dòng)提取。以上研究所使用的數(shù)據(jù)源均來自國外遙感衛(wèi)星，利用我國高分衛(wèi)星進(jìn)行水體提取的研究和算法還不多見。

高分一號(hào)衛(wèi)星(GF-1)于2013年4月26日發(fā)射成功，是我國第一顆考核壽命大于5 a的衛(wèi)星，衛(wèi)星除裝載有兩臺(tái)2 m分辨率全色和8 m分辨率多光譜相機(jī)外，還搭載了4臺(tái)16 m分辨率多光譜相機(jī)。6臺(tái)相機(jī)的配合使得我國首次實(shí)現(xiàn)了單顆衛(wèi)星高分辨率與大幅寬技術(shù)的有效結(jié)合[14]，大大提高了我國高分辨率影像數(shù)據(jù)的自給率[15]，將為我國的科技發(fā)展帶來不可估量的影響[16]。

受天氣影響，遙感影像在地表水體提取過程中易受到云層的影響。遙感影像上受云層影響的區(qū)域覆蓋了大面積的地表信息，使得圖像信息不完善，給影像數(shù)據(jù)的處理和應(yīng)用造成了極大的障礙[17-18]。如何消除遙感影像上云層的影響，至今仍是研究的熱門話題。目前已有不少國內(nèi)外相關(guān)專家和學(xué)者對(duì)此展開研究，研究結(jié)果表明：同態(tài)濾波法適用于大面積有薄云存在的遙感影像[19-21]，對(duì)厚云區(qū)域基本不產(chǎn)生效果；多源數(shù)據(jù)融合法使用不同傳感器的遙感影像，將無云部分替代其他影像的云干擾部分[22]；時(shí)間序列傅立葉分析法簡化了壓縮時(shí)間序列遙感數(shù)據(jù)，對(duì)遙感影像中的云進(jìn)行檢測、去除和替代處理[23]。除此之外，將遙感影像與地形信息有效結(jié)合的研究也有很多[24-25]，王鑫蕊等[26]利用湖泊邊界位于同一高度的特點(diǎn)，提出了一種利用湖泊的已知部分邊界結(jié)合DEM數(shù)據(jù)，還原完整邊界的自動(dòng)方法，有效克服了云層的干擾。

屯溪流域?yàn)橹行『恿髦攸c(diǎn)研發(fā)項(xiàng)目的典型示范流域，東方紅水庫是屯溪流域最大的一座中型水庫，其庫容、面積變化影響著該流域的應(yīng)急防控技術(shù)。由于遙感數(shù)據(jù)周期間隔較長，無法實(shí)時(shí)獲取水庫的面積信息，且研究區(qū)內(nèi)影像常年或多或少地受到云層的影響，降低了遙感數(shù)據(jù)的實(shí)用性。為了利用我國現(xiàn)有的光學(xué)遙感影像提高遙感數(shù)據(jù)的質(zhì)量，本文以東方紅水庫為研究對(duì)象，提出一種GF-1遙感影像結(jié)合DEM消除云層干擾的連續(xù)水體重建法，對(duì)于連續(xù)水體局部受云層影響的影像，結(jié)合GF-1遙感影像和高精度等高線數(shù)據(jù)，通過地理配準(zhǔn)還原遙感影像水體面積，從而實(shí)現(xiàn)云層干擾下水體的自動(dòng)提取。

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

1.1 研究區(qū)概況

東方紅水庫位于新安江水系二級(jí)支流東亭河(虞山溪)上(圖1)，是以防洪、灌溉為主，結(jié)合發(fā)電、旅游供水等綜合利用的中型水利樞紐工程。水庫集水區(qū)域處在117°E～119°E、29°N～31°N之間，集水面積為60.0 km2，主河道長11.8 km。東方紅水庫總庫容0.244 5億m3，設(shè)計(jì)灌溉面積0.3萬hm2，電站總裝機(jī)容量400 kW，多年平均發(fā)電量131萬kW·h。工程規(guī)模為安徽省重點(diǎn)中型水庫，工程等別為Ⅲ等，攔河壩、溢洪道、引水涵洞，以及由導(dǎo)流隧洞改建的放水底孔等均為3級(jí)建筑物。

圖1 研究區(qū)位置和地形Fig.1 Location and topography of the study area

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

遙感數(shù)據(jù)來自中國資源衛(wèi)星應(yīng)用中心(http：//www.cresda.com/CN/index.shtml)的高分一號(hào)寬幅傳感器(GF1_WFV)數(shù)據(jù)，影像空間分辨率為16 m。通過實(shí)地查勘的方式，收集了東方紅水庫高精度等高線數(shù)據(jù)及2015年水庫詳細(xì)運(yùn)行資料。2015年為豐水年，東方紅水庫集水區(qū)域平均降水量為2 411.6 mm。通過2020年2月陸面過程資料分布式數(shù)據(jù)存檔中心LP DAAC(Land Processes Distributed Active Archive Center)公布的1 rad/s分辨率NASADEM數(shù)據(jù)集產(chǎn)品(https：//lpdaac.usgs.gov/news/release-nasadem-data-products/)提取研究流域包含的水體區(qū)域。為體現(xiàn)水庫在不同時(shí)期的蓄水面積信息，選取2015年非雨期(4、10月)、雨期(5、8、9月)共7景衛(wèi)星影像。

2 研究方法

2.1 數(shù)據(jù)處理

2.1.1 遙感影像數(shù)據(jù)處理

選擇中國資源衛(wèi)星應(yīng)用中心公布的1A級(jí)影像數(shù)據(jù)產(chǎn)品。為降低遙感影像在成像過程中大氣傳輸、傳感器本身、地球曲率等多種因素的影響，需要對(duì)原始遙感影像進(jìn)行預(yù)處理。結(jié)合IDL(interactive data language)與ENVI(the environment for visualizing images)對(duì)數(shù)據(jù)進(jìn)行批量預(yù)處理，預(yù)處理的方式主要包括：輻射定標(biāo)、大氣校正[27]、正射校正等。預(yù)處理結(jié)束后，將研究區(qū)域?qū)?yīng)的矢量文件與遙感影像配準(zhǔn)，完成研究區(qū)域的裁剪。2015年7景水體局部受云層影響的遙感影像數(shù)據(jù)見表1。

表1 高分一號(hào)衛(wèi)星遙感影像數(shù)據(jù)集

2.1.2 高精度等高線數(shù)據(jù)處理

利用線性均分方式將東方紅水庫對(duì)應(yīng)的高精度等高線數(shù)據(jù)進(jìn)行插值，并選取多幅高保證率歷史遙感影像，通過最大類間方差膨脹迭代法提取水體，從而得到水庫不同水位值對(duì)應(yīng)的水體邊緣輪廓線，以此為參考，利用GIS軟件對(duì)插值得到的高精度等高線數(shù)據(jù)進(jìn)行人工解譯，精準(zhǔn)修正等高線數(shù)據(jù)，并確保高程數(shù)據(jù)完全覆蓋水體范圍。

2.2 水體提取與云層消除

2.2.1 最大類間方差膨脹迭代法提取水體

歸一化差分水體指數(shù)(normalized difference water index，NDWI)法[9]利用遙感影像的綠波與近紅外波兩個(gè)特定波段進(jìn)行歸一化差值處理，從而突顯影像中的水體。將經(jīng)過預(yù)處理的遙感影像利用NDWI法初步提取水體。

最大類間方差法由日本學(xué)者Ostu等[28]提出，適用于雙峰情況的閾值分割，可以實(shí)現(xiàn)自動(dòng)求取閾值的目的。該方法根據(jù)圖像的灰度值，將圖像分為背景和目標(biāo)兩部分。分割閾值從最小灰度值到最大灰度值遍歷每個(gè)值，不同的閾值分割使得背景和目標(biāo)之間的方差不同，類間方差越大說明被閾值分割的兩部分差別越大。若出現(xiàn)目標(biāo)與背景的錯(cuò)分、漏分等現(xiàn)象，均會(huì)導(dǎo)致兩部分差別變小。因此，利用類間方差最大時(shí)對(duì)應(yīng)的閾值進(jìn)行分割，被錯(cuò)分的概率應(yīng)該是最小的。

利用NDWI法提取未受云層影響的局部水體，將初始閾值T0設(shè)為-0.02，若歸一化差分水體指數(shù)INDWIT0，得到的結(jié)果為初步提取的水體。以此基于形態(tài)學(xué)膨脹算法建立水體的緩沖區(qū)。當(dāng)膨脹算法中的結(jié)構(gòu)元素大小發(fā)生變化時(shí)，可得到不同范圍的緩沖區(qū)，結(jié)構(gòu)元素過大、過小都會(huì)造成水體的誤提。根據(jù)歷史經(jīng)驗(yàn)，該研究區(qū)選取大小為15×15的結(jié)構(gòu)元素膨脹提取效果最佳，故以該結(jié)構(gòu)元素大小對(duì)初步提取的部分水體進(jìn)行膨脹，建立緩沖區(qū)，對(duì)緩沖區(qū)內(nèi)包含的區(qū)域應(yīng)用最大類間方差法確定水體分割閾值，從而實(shí)現(xiàn)水體的提取。以前一次提取的水體區(qū)域作為后一次建立緩沖區(qū)的初始范圍，重復(fù)利用最大類間方差法進(jìn)行多次迭代，自適應(yīng)地確定水體與周邊地物的分割閾值，由此得到每次迭代確定的水體面積，當(dāng)相鄰兩次水體面積差值絕對(duì)值滿足要求精度時(shí)，迭代結(jié)束。

|Sk-Sk-1|<0.01 (k=1，2,3,…)

(1)

S0=Sk

(2)

式中：Sk——第k次迭代得到的水體面積；Sk-1——第k-1次迭代得到的水體面積；S0——最大類間方差膨脹迭代法提取到的最優(yōu)水體面積。

2.2.2 基于高精度等高線數(shù)據(jù)的云層影響消除法

通過最大類間方差膨脹迭代得到東方紅水庫的局部水體范圍，與處理后的高精度等高線數(shù)據(jù)進(jìn)行地理坐標(biāo)配準(zhǔn)，可得到每個(gè)水體柵格對(duì)應(yīng)的水位高程值。某一時(shí)刻連續(xù)水體水位應(yīng)為唯一值，而由于云層的干擾以及誤差影響導(dǎo)致遙感影像提取的邊界水位并不唯一，故對(duì)云層干擾下由遙感數(shù)據(jù)得到的水體不同水位進(jìn)行數(shù)學(xué)統(tǒng)計(jì)計(jì)算。

定義受云層影響遙感影像得到的局部水體為W類水體，局部水體邊界網(wǎng)格以內(nèi)的內(nèi)部水體網(wǎng)格為W1，水庫原始外輪廓邊緣線對(duì)應(yīng)的網(wǎng)格為W2，受云層影響而產(chǎn)生的邊界線對(duì)應(yīng)的水體網(wǎng)格為W3，其他非水體網(wǎng)格為W0。

若W2類水體網(wǎng)格對(duì)應(yīng)的不同水位值共n個(gè)，對(duì)此進(jìn)行均值計(jì)算：

(3)

SM⊕SV=SV

(4)

式中：SM——受云層影響遙感影像得到的局部水體面積；SV——還原水體水位值所包圍的區(qū)域面積。兩類水體經(jīng)地理配準(zhǔn)后可實(shí)現(xiàn)云層干擾下的水體還原提取。

3 結(jié)果與分析

在IDL + ENVI編程環(huán)境中實(shí)現(xiàn)遙感影像數(shù)據(jù)的預(yù)處理。將經(jīng)過預(yù)處理的圖像在MATLAB中進(jìn)行最優(yōu)閾值的確定，再通過GIS軟件確定最優(yōu)閾值對(duì)應(yīng)的水體范圍，與水體范圍對(duì)應(yīng)的高精度等高線數(shù)據(jù)進(jìn)行地理配準(zhǔn)，從而實(shí)現(xiàn)云層干擾下水體的還原處理。選取GF-1_2015-09-08景遙感影像(云量為11%)進(jìn)行有云水體提取，各步驟處理結(jié)果如圖2所示。

圖2 水體GF-1遙感影像云消除過程Fig.2 Cloud removal process in water body extraction for the GF-1 satellite imageries

圖2(a)為東方紅水庫集水區(qū)域2015-09-08GF-1_2015-09-08景遙感影像經(jīng)過輻射定標(biāo)、大氣校正、正射校正、配準(zhǔn)及研究區(qū)裁剪等預(yù)處理后的圖像。由圖2(a)可以看出，紅色區(qū)域?yàn)樗畮焖w區(qū)域，部分水體受到云層影響。圖2(b)為云層干擾下通過最大類間方差膨脹迭代提取的東方紅水庫水體面積圖，由于受云層影響較為嚴(yán)重，水庫外輪廓線提取并不完整，且水庫內(nèi)部也缺失部分水體。圖2(c)為東方紅水庫云干擾水體提取圖與其高精度等高線數(shù)據(jù)的地理配準(zhǔn)圖，圖中水體與等高線數(shù)據(jù)產(chǎn)生輕微偏移，水位值并不唯一。圖2(d)為通過地理配準(zhǔn)還原云層干擾下的東方紅水庫水體面積圖，較為清晰地還原了水庫的基本輪廓。

根據(jù)東方紅水庫289～297 m水位區(qū)間對(duì)應(yīng)的水位-面積曲線(圖3)以及觀測資料記錄的每日水庫實(shí)測水位值(表2)，可以推導(dǎo)出對(duì)應(yīng)時(shí)間水庫的水體面積，以此作為SV的參考值。

圖3 東方紅水庫水位-面積關(guān)系曲線Fig.3 Water level-area relationship curve of Dongfanghong Reservoir

從表2可以看出，東方紅水庫GF-1_2015-09-08景遙感影像實(shí)測水位值h為291.35 m，根據(jù)水位-面積關(guān)系曲線得到的參考面積S為1.995 6 km2。然而由于云層的干擾，根據(jù)最大類間方差膨脹迭代得到的水體面積為1.060 7 km2，面積差值為0.934 9 km2，差距較大。按照本文提出的通過地理配準(zhǔn)還原云層干擾水庫水體面積，可得到水庫水位還原值為291. 95 m，還原后庫區(qū)水體面積為2.042 km2，面積差值為0.046 4 km2，可以看出，該方法簡單有效地還原了水庫水體面積。

表2 消除云層干擾后的庫區(qū)水體遙感提取面積

另選取東方紅水庫2015年6景影像并使用上述方法進(jìn)行云消除處理。圖4為最大類間方差膨脹迭代提取的水體范圍與其等高線數(shù)據(jù)的地理配準(zhǔn)圖，藍(lán)色區(qū)域?yàn)镾0。從圖4可以看出云層干擾對(duì)東方紅水庫水體提取范圍影響較大，嚴(yán)重時(shí)水庫外輪廓基本無法確定。除此之外還可以看出，庫區(qū)水體外輪廓線對(duì)應(yīng)的水位值并不唯一，因此需要通過數(shù)學(xué)統(tǒng)計(jì)計(jì)算方式確定需要還原的水位值。

圖4 地理配準(zhǔn)后的衛(wèi)星影像Fig.4 Satellite imageries after geographic referencing

圖5為經(jīng)地理配準(zhǔn)后還原的水體，與等高線數(shù)據(jù)相結(jié)合經(jīng)均值計(jì)算后，水庫對(duì)應(yīng)的水位值可唯一確定，最終結(jié)果表明該方法簡單有效地消除了云層的干擾，較準(zhǔn)確地還原了東方紅水庫在非雨期(4、10月)、雨期(5、8、9月)云層影響下的水體。從圖5(a)(b)(f)可以看出，由于水庫在4、10月來水量較少，隨著下游水量的補(bǔ)給，庫區(qū)水體范圍逐漸減小，提取的水體范圍較??；圖5(c)(d)(e)表明，5、8、9月降水量較多，水庫來水量大于出水量，蓄水量逐漸增多，故水體面積增大，水位升高，庫區(qū)之間連接的細(xì)小河道蓄水位也逐漸升高。

圖5 水體還原結(jié)果Fig.5 Results of extracted water bodies

從表2可以看出，S0與S相比差距較大，其中GF-1_20150426景遙感影像水體面積相差最大，達(dá)1.368 4 km2，難以反映真實(shí)水體。選取的7景影像經(jīng)地理配準(zhǔn)后得到的水體還原值與實(shí)測水位參考值接近，還原后的水體水位平均差值為0.438 6 m，有效水體面積平均差值為0.080 9 km2，平均相對(duì)誤差在5%以內(nèi)，較為精確地還原了庫區(qū)的水體。其中GF-1_20150512景遙感影像水體庫區(qū)面積差值僅為0.008 3 km2，較好地消除了云層干擾，還原了水體水位值及對(duì)應(yīng)的面積。

從表2還可以看出，GF-1_20150401、GF-1_20150426和GF-1_20151003景遙感影像還原的水位值與實(shí)測水位值差距略微偏高，這可能是由于所用影像分辨率較低，為16 m，在預(yù)處理過程中水體范圍整體產(chǎn)生輕微偏移，使得提取的水體存在部分誤差，且該方法對(duì)地形數(shù)據(jù)精度要求較高，本文將1 m等高線經(jīng)線性插值、人工修訂后得到的等高線與提取的遙感信息進(jìn)行地理配準(zhǔn)，不同數(shù)據(jù)源之間易出現(xiàn)誤差疊加現(xiàn)象。

4 結(jié) 語

高分一號(hào)衛(wèi)星遙感影像受云層影響限制了數(shù)據(jù)的實(shí)用性，且在云層干擾下通過最大類間方差膨脹迭代法得到的水體并不完整，面積平均差值為1.009 4 km2，無法準(zhǔn)確反應(yīng)水體的真實(shí)信息，難以直接用于陸地表面水體的提取。通過本文提出的地理配準(zhǔn)還原云層干擾下水體的方法，還原后的水體水位平均差值為0.438 6m，有效水體面積平均差值為0.080 9 km2，最小僅為0.008 3 km2，平均相對(duì)誤差在5%以內(nèi)，有效地克服了云層干擾，較為精確地還原了庫區(qū)水體。該方法簡單有效，與高精度等高線數(shù)據(jù)相結(jié)合，通過簡單的圖像計(jì)算便可克服云層的干擾，提高了衛(wèi)星遙感影像的實(shí)用性，可廣泛用于地表水資源的動(dòng)態(tài)提取，為水資源調(diào)查及地表水體動(dòng)態(tài)監(jiān)測等領(lǐng)域提供高質(zhì)量的遙感影像數(shù)據(jù)。

通過本文方法得到的水體面積雖然有效克服了云層的干擾，但該方法對(duì)地形數(shù)據(jù)精度要求較高，對(duì)于缺少水下高程數(shù)據(jù)的水體來說，水位動(dòng)態(tài)變化值常在某范圍內(nèi)上下變動(dòng)，故可根據(jù)多景無云歷史遙感影像(包含豐水期、枯水期較多高程值變化)得到的水位值對(duì)水下空缺的高程值進(jìn)行填補(bǔ)，再進(jìn)一步進(jìn)行人工修正，但精度方面問題有待提高。除此之外，對(duì)于一些復(fù)雜的水體情況，如湖心島等，該方法需要結(jié)合人工解譯才能進(jìn)行水體提取，對(duì)此仍需要進(jìn)一步的深入研究。