李 穎， 陳懷亮， 葉昊天， 周兆基

(1.中國(guó)氣象局·河南省農(nóng)業(yè)氣象保障與應(yīng)用技術(shù)重點(diǎn)開(kāi)放實(shí)驗(yàn)室，鄭州 450003；2.河南省氣象科學(xué)研究所，鄭州 450003； 3.哈爾濱市氣象局，哈爾濱 150028； 4.河南省氣象局，鄭州 450003)

引 言

近年來(lái)，全球氣候變化加劇[1]，極端天氣現(xiàn)象頻繁出現(xiàn)。2021年7月17日至22日,河南省出現(xiàn)了歷史罕見(jiàn)的極端暴雨過(guò)程[2]。由于強(qiáng)降水引發(fā)的洪澇災(zāi)害發(fā)生頻率較高、影響范圍較大，因此在災(zāi)后快速獲取大面積洪水分布信息對(duì)于災(zāi)害評(píng)估、救援等具有重要意義。遙感影像具有宏觀、客觀、及時(shí)獲取地表信息的優(yōu)勢(shì)，可以克服災(zāi)情現(xiàn)場(chǎng)調(diào)查時(shí)交通不便和淹沒(méi)區(qū)域難以抵達(dá)的困難，節(jié)約人力成本和時(shí)間成本，在洪水面積提取方面具有重要的應(yīng)用價(jià)值。

晴空條件下中高分辨率光學(xué)遙感影像在水體信息提取方面具有較高精度，已得到廣泛應(yīng)用[3-4]。由于光學(xué)影像難以穿透云層，而洪水發(fā)生前后通常為云雨霧天氣，因此使用光學(xué)遙感影像難以及時(shí)獲取地表信息。合成孔徑雷達(dá)(SAR)具有穿透云層的能力，并可通過(guò)地物在不同極化方式下的后向反射特性識(shí)別地表類型，在災(zāi)害性天氣發(fā)生過(guò)程中可以全天時(shí)、全天候持續(xù)監(jiān)測(cè)地表變化情況[5]?，F(xiàn)有研究中，將SAR影像用于水體信息提取的技術(shù)已較為成熟，最廣泛應(yīng)用的方法主要有基于極化灰度圖、干涉相干圖、紋理圖等的閾值法[6-8]和自動(dòng)聚類算法[9-10]。相關(guān)研究包括Marti-Cardona等[11]通過(guò)入射角分析識(shí)別濕地洪水區(qū)和洪水造成的開(kāi)放水體；孫亞勇等[12]通過(guò)山體陰影剔除，提高了山區(qū)水體提取的精度；谷鑫志等[13]在閾值分割的基礎(chǔ)上，綜合利用了GF-3影像的多極化信息與空間上下文信息，利用最大后驗(yàn)概率準(zhǔn)則輸出水體分布圖等；眭海剛等[14]利用深度特征和語(yǔ)義信息實(shí)現(xiàn)了災(zāi)前光學(xué)影像和災(zāi)后SAR影像的高精度配準(zhǔn),進(jìn)而實(shí)現(xiàn)洪水變化監(jiān)測(cè)和災(zāi)損信息提?。粐?yán)夏青等[15]基于Sentinel-1數(shù)據(jù)提取長(zhǎng)江干流水域面積并開(kāi)展動(dòng)態(tài)變化監(jiān)測(cè)。

總體來(lái)看，在基于SAR影像的水體識(shí)別和洪水面積提取方面，國(guó)內(nèi)研究成果相對(duì)較少。一方面，對(duì)包含水體在內(nèi)的典型地物進(jìn)行后向散射特征分析時(shí)，多數(shù)研究通常僅針對(duì)不同極化方式進(jìn)行分析，對(duì)局地入射角對(duì)典型地物后向散射特性影響的研究較少；另一方面，在洪水發(fā)生后往往僅對(duì)新增開(kāi)放水體進(jìn)行識(shí)別，沒(méi)有開(kāi)展新增濕地的識(shí)別，故提取洪水區(qū)面積偏小。針對(duì)上述問(wèn)題，本研究統(tǒng)計(jì)分析了不同典型地物的后向散射系數(shù)隨局地入射角的變化情況，并在洪水發(fā)生后提取新增水體的基礎(chǔ)上，進(jìn)行了新增濕地與易混淆地物的極化分解分析，進(jìn)而提取了新增濕地，綜合新增開(kāi)放水體與新增濕地得到更為完整的洪水區(qū)面積。

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

1.1 研究區(qū)

美國(guó)2018年的Florence颶風(fēng)為30年來(lái)威脅南北卡羅來(lái)納州的最強(qiáng)級(jí)別颶風(fēng)。該颶風(fēng)于2018年9月13日登陸美國(guó)東海岸，隨后給南北卡羅來(lái)納州帶來(lái)了強(qiáng)降水和災(zāi)難性洪水，其洪水至9月26日仍未完全消退。本文研究區(qū)為覆蓋南北卡羅來(lái)納州部分地區(qū)的一景Sentinel-1影像，包含Catawba河和Congaree河兩條由Florence颶風(fēng)造成的洪水泛濫的河流的部分區(qū)域。

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

1.2.1 Sentinel-1數(shù)據(jù)

Sentinel-1載荷為C波段SAR傳感器，工作頻率為5.4 GHz，由Sentinel-1A和Sentinel-1B兩顆衛(wèi)星組網(wǎng)，單顆星重訪周期為12 d，雙星組網(wǎng)的重訪周期為6 d。Sentinel-1的成像模式包括條帶、干涉寬幅、超寬幅和波浪。本文共使用了6幅Sentinel-1A和Sentinel-1B干涉寬幅模式的level-1級(jí)地距影像(Ground Range Detected，GRD)產(chǎn)品，該產(chǎn)品采用WGS84橢球?qū)?cè)視成像的斜距圖像校正為地距影像，使影像的方位向和距離向的分辨率一致，均為10 m，具有VV和VH兩種極化方式。根據(jù)數(shù)據(jù)質(zhì)量和軌道情況(選擇同軌道數(shù)據(jù))，影像成像時(shí)間為2018年7月26日、8月7日、8月19日、8月31日、9月12日和9月18日，包含洪水前的5景影像和泛濫期的1景影像，以支持研究區(qū)Florence颶風(fēng)造成的洪水淹沒(méi)遙感監(jiān)測(cè)。

1.2.2 雷達(dá)數(shù)據(jù)預(yù)處理

使用SNAP軟件和ENVI軟件對(duì)Sentinel-1 level-1級(jí)地距影像進(jìn)行預(yù)處理，主要步驟包括熱噪聲消除、多視處理、輻射定標(biāo)、斑噪去除、地形校正和直方圖拉伸。其中，多視處理的等效視數(shù)設(shè)為4.9，輻射定標(biāo)得到Sigma后向散射系數(shù)，斑噪去除選擇Gamma-Map濾波器(濾波窗口設(shè)為5×5)，地形校正采用SRTM DEM數(shù)據(jù)(空間分辨率為90 m)。詳細(xì)的處理方法見(jiàn)參考文獻(xiàn)[16-17]。預(yù)處理后斑點(diǎn)噪聲得到較好抑制，水域邊界清晰光滑，9月18日的影像較洪水前的2景影像可見(jiàn)明顯的泛濫區(qū)(圖1)。

圖1 預(yù)處理后研究區(qū)局部地區(qū)VV極化、VH極化、VV極化RGB合成圖

1.2.3 NLCD2016數(shù)據(jù)及預(yù)處理

本文使用美國(guó)地質(zhì)調(diào)查局(USGS)發(fā)布的2016年國(guó)家土地覆被數(shù)據(jù)庫(kù)(NLCD)分類數(shù)據(jù)(NLCD2016)作為洪水前的地表真實(shí)覆被數(shù)據(jù)，用于SAR圖像中典型地物特征分析與解譯。裁剪出覆蓋研究區(qū)范圍的NLCD2016數(shù)據(jù)，根據(jù)研究區(qū)實(shí)際情況對(duì)地物類別進(jìn)行合并，歸納為水體、裸土、建成區(qū)、農(nóng)作物、灌木、森林和濕地等7種典型地物(圖2)。

圖2 研究區(qū)土地利用類型分類圖

2 典型地物特征分析

2.1 后向散射系數(shù)均值分析

利用NLCD2016數(shù)據(jù)對(duì)洪水前的5景影像進(jìn)行分類，計(jì)算不同典型地物在兩種極化方式下后向散射系數(shù)的均值(表1)。由表1可看出，各種地物類型VV極化的后向散射系數(shù)均高于VH極化的后向散射系數(shù)，同時(shí)VV極化后向散射系數(shù)較高的地物VH極化后向散射系數(shù)也較高，反之亦然。兩種極化方式下，水體后向散射系數(shù)的均值均最低，其次較低的是裸土的，最高的則是濕地的。在VV極化下，建成區(qū)、農(nóng)作物、灌木、森林等均值的可分性較差；在VH極化下，森林后向散射系數(shù)的均值則明顯高于另外三種地物的均值。

表1 典型地物的后向散射系數(shù)均值 dB

2.2 局地入射角分析

在不考慮角度效應(yīng)的情況下，使用VV和VH兩種極化方式，水體、裸土、森林、濕地的可分性相對(duì)較好。但考慮地形起伏后，雷達(dá)波的局地入射角(Local Incidence Angle, LIA)亦會(huì)在一定范圍內(nèi)變化。LIA是影響地物散射特性的一個(gè)重要因素。本文統(tǒng)計(jì)不同典型地物的后向散射系數(shù)隨LIA的變化情況，結(jié)果顯示，LIA越小，地物的回波強(qiáng)度越大；當(dāng)LIA在0°～90°均勻變化時(shí)，所有地物類型在VV和VH極化方式下均無(wú)法僅使用后向散射系數(shù)進(jìn)行有效區(qū)分(圖3)。

圖3 典型地物的后向散射系數(shù)隨局地入射角的變化情況

實(shí)際上，大部分星載雷達(dá)的入射角設(shè)計(jì)在40°左右，因此雷達(dá)波的LIA均勻變化在0°～90°的可能性較小。本文統(tǒng)計(jì)了研究區(qū)6景雷達(dá)影像中LIA的真實(shí)分布情況，統(tǒng)計(jì)結(jié)果顯示，各景影像中90%左右像元的雷達(dá)波LIA分布在30°～45°(圖4)。分析LIA在30°～45°時(shí)典型地物的后向散射系數(shù)特征可見(jiàn)，水體在兩種極化方式下較其他地物均有較好的可分性，特別是在VV極化方式下，使用后向散射系數(shù)可有效識(shí)別開(kāi)放水體。

圖4 局地入射角在30°～45°時(shí)典型地物

2.3 紋理特征分析

現(xiàn)有基于雷達(dá)影像的地物分類、識(shí)別研究中，紋理特征使用較為廣泛。本文亦使用灰度共生矩陣(GLCM)分析研究區(qū)各類典型地物的多種紋理特征。GLCM是一種通過(guò)圖像灰度空間的相關(guān)性描述紋理的方法。其定義為：在圖像內(nèi)，統(tǒng)計(jì)所有從灰度級(jí)為i的像元點(diǎn)，沿移動(dòng)方向移動(dòng)一定距離后，到達(dá)灰度級(jí)為j的像元的概率為p(i,j)時(shí)所形成的方陣。不同地物類型的7種GLCM紋理特征，包括方差(Variance)、均質(zhì)性(Homogeneity)、對(duì)比度(Contrast)、熵(Entropy)、角二階矩(ASM)和相關(guān)性(Correlation)，計(jì)算公式如下：

(1)

(2)

(3)

(4)

(5)

(6)

其中，N為圖像的量化等級(jí)，μ為GLCM的均值，σ為GLCM的標(biāo)準(zhǔn)差，其計(jì)算公式如下：

(7)

(8)

本文計(jì)算GLCM紋理特征時(shí)選擇3×3的窗口，移動(dòng)方向設(shè)為0°、45°、90°和135°，圖像量化等級(jí)設(shè)為64級(jí)，步長(zhǎng)設(shè)為1。分別統(tǒng)計(jì)7種典型地物在VV極化方式下的6種紋理特征(表2)，結(jié)果顯示：建成區(qū)在方差和對(duì)比度上明顯大于其他地物類型的；水體則在6種紋理特征上均區(qū)別于其他地物類型的，特別是方差與對(duì)比度明顯小于其他地物類型的；在均質(zhì)性、熵和角二階矩上，裸土和建成區(qū)的均值接近，而農(nóng)作物、灌木、森林和濕地的均值接近。由以上分析可知，紋理特征可有效應(yīng)用于水體識(shí)別，并可幫助區(qū)分植被與非植被。

表2 典型地物的6種紋理特征

3 洪水區(qū)域提取方法

3.1 ISODATA算法

ISODATA算法屬于非監(jiān)督分類方法，是在K-均值算法的基礎(chǔ)上改進(jìn)的迭代自組織聚類方法，該方法在開(kāi)放水體的提取中已得到廣泛應(yīng)用[18]。通過(guò)上文分析，使用VV和VH兩種極化方式下的后向散射系數(shù)與方差和對(duì)比度兩種紋理特征，可有效識(shí)別開(kāi)放水體和其他地物類型，故選擇上述4種特征值進(jìn)行ISODATA迭代分類。分類后進(jìn)行類別合并，將水體作為一類，其他地物類型合并為一類，之后運(yùn)用數(shù)學(xué)形態(tài)學(xué)算子，通過(guò)聚類處理消除開(kāi)放水體中的斑點(diǎn)。

3.2 極化分解

發(fā)生洪水后，一方面包括裸土、道路、農(nóng)作物、牧草等在內(nèi)的低矮地物可能被淹沒(méi)形成開(kāi)放水體，另一方面森林和灌木的底部可能被淹沒(méi)，從而形成濕地。利用Freeman極化目標(biāo)分解理論分析洪水發(fā)生后由森林和灌木變化形成的濕地區(qū)域的后向散射特征變化。Freeman分解由Freeman和Durden于1998年提出[19-20]，該方法在假定散射體滿足反射對(duì)稱性的條件下，對(duì)極化協(xié)方差矩陣建立包括適度粗糙表面的Bragg散射項(xiàng)、隨機(jī)方向偶極子組成的體散射項(xiàng)和由正交平面構(gòu)成的偶次散射項(xiàng)等3個(gè)分量的基本散射機(jī)制模型。Bragg散射的典型地物包括水體、裸土、硬化地面等，體散射的典型地物包括森林、灌木等，偶次散射的典型地物則包括復(fù)雜建筑物、建筑與地面、冠層與地面、冠層與水面等能夠形成角散射的地物組合。其中，冠層與水面形成角散射的情況主要發(fā)生在濕地。

進(jìn)行Freeman分解時(shí)，森林和灌木的主要散射分量為體散射，另外兩種分量貢獻(xiàn)較小，在非定量研究中較之體散射可忽略。濕地則存在來(lái)自森林、灌木冠層的體散射和冠層與水面之間的偶次散射兩種主要散射分量，且其偶次散射項(xiàng)的回波信號(hào)相比體散射的后向散射系數(shù)較高。因此，針對(duì)洪水發(fā)生前后的森林和灌木，消除回波信號(hào)中的體散射項(xiàng)后，若仍存在較高的回波信號(hào)，則說(shuō)明發(fā)生了偶次散射，即被洪水淹沒(méi)。Freeman分解中體散射的極化協(xié)方差矩陣形式為：

(9)

式中，fv表示體散射分量的權(quán)重。

在Freeman分解理論中，體散射的同極化后向散射系數(shù)是交叉極化的3倍。因偶次散射項(xiàng)理論上不存在交叉極化，因此應(yīng)用Freeman極化目標(biāo)分解理論，計(jì)算極化分解系數(shù)

σ3v=SVV-3SVH[21-22]

(10)

其中，SVV為VV同極化后向散射能量，SVH為VH交叉極化后向散射能量。通過(guò)計(jì)算σ3v，體散射項(xiàng)可以得到較大程度的消除，而偶次散射項(xiàng)不受影響，即森林與灌木的極化分解系數(shù)σ3v很低，而濕地的σ3v較森林與灌木明顯較高。表3為洪水發(fā)生前后森林、灌木和濕地的極化分解系數(shù)。由表3可看出，洪水發(fā)生前5幅影像中森林的平均σ3v最低，其次是灌木的，濕地的平均σ3v明顯高于前兩者的。洪水發(fā)生后，原森林區(qū)域和灌木區(qū)域的平均σ3v明顯升高，顯示研究區(qū)域中確有大面積的灌木和森林因洪水淹沒(méi)變?yōu)闈竦?。在上述分析的基礎(chǔ)上，通過(guò)判讀洪水前森林、灌木、濕地的極化分解系數(shù)的變化特征，結(jié)合VV后向散射系數(shù)設(shè)置變化檢測(cè)閾值，可提取洪水后形成的濕地。

表3 洪水發(fā)生前后典型地物的極化分解系數(shù)

4 結(jié)果與分析

4.1 開(kāi)放水體提取

利用ISODATA算法和VV、VH兩種極化方式下的后向散射系數(shù)與方差、對(duì)比度兩種紋理特征，分別提取洪水發(fā)生前雷達(dá)影像中的開(kāi)放水體和洪水發(fā)生后雷達(dá)影像中的開(kāi)放水體，計(jì)算其差值并做去噪聲點(diǎn)等后處理，即得到洪水發(fā)生后新增開(kāi)放水體面積(共計(jì)449938個(gè)像元，新增水體面積45 km2)，局部提取效果見(jiàn)圖5，圖例與圖2相同。洪水淹沒(méi)區(qū)域主要分布在Catawba河沿岸，Congaree河沿岸也有部分區(qū)域被淹沒(méi)。對(duì)照NLCD2016土地覆被數(shù)據(jù)，淹沒(méi)區(qū)域的土地利用類型主要是農(nóng)作物、灌木和濕地。

圖5 洪水后新增開(kāi)放水體提取結(jié)果

4.2 新增濕地提取

在分析洪水發(fā)生前后森林、灌木、濕地的極化分解系數(shù)變化特征的基礎(chǔ)上，排除洪水發(fā)生前的濕地面積，利用閾值判別法和VV極化方式下的后向散射系數(shù)提取新增濕地。因新增濕地區(qū)域出現(xiàn)在洪水前的森林和灌木土地覆蓋類型上，通過(guò)典型地物雷達(dá)回波特征分析，這兩類土地覆蓋類型在變化為濕地后，VV極化后向散射系數(shù)將明顯提升。結(jié)合圖3與圖4中典型地物的后向散射系數(shù)隨局地入射角的變化和分布情況，在極化分解原理的基礎(chǔ)上，經(jīng)人機(jī)交互實(shí)驗(yàn)，提取洪水發(fā)生后森林和灌木VV極化后向散射系數(shù)值增大3 dB以上的區(qū)域?yàn)樾略鰸竦孛娣e，做去噪聲點(diǎn)等后處理，共計(jì)得到2124106個(gè)像元，面積約212.4 km2，局部提取效果見(jiàn)圖6，圖例與圖2相同。新增濕地區(qū)域主要分布在現(xiàn)有濕地的周邊。

圖6 洪水后新增濕地提取結(jié)果

5 結(jié)論與討論

本研究將洪水發(fā)生后淹沒(méi)區(qū)分為兩種情況：一種是發(fā)生洪水后，裸土、道路、農(nóng)作物、灌木、濕地等在內(nèi)的低矮地物被淹沒(méi)形成的開(kāi)放水體;另一種是發(fā)生洪水后，森林和灌木下部被淹沒(méi)形成新的濕地。通過(guò)典型地物特征分析，開(kāi)放水體在VV和VH兩種極化方式下較其他地物具有良好的可分性，其方差與對(duì)比度兩種紋理特征也明顯小于其他地物類型的。利用ISODATA算法進(jìn)行非監(jiān)督分類提取水體，結(jié)果顯示，洪水淹沒(méi)區(qū)域面積約45 km2，主要分布在Catawba河沿岸和Congaree河沿岸，即新增開(kāi)放水體出現(xiàn)在原有開(kāi)放水體的鄰近位置，分析原因可能是河道水位升高，進(jìn)而洪水漫灌或因河谷地帶地勢(shì)較低容易積水，洪水匯聚形成開(kāi)放水體。

通過(guò)典型地物特征分析，濕地在VV極化方式下后向散射系數(shù)的均值明顯高于其他地物類型的均值，但其紋理特征與其他地物類型，特別是與森林和灌木的可分性較差。采用極化分解分析顯示，洪水發(fā)生后原森林區(qū)域和灌木區(qū)域的平均極化分解系數(shù)明顯升高，表明研究區(qū)域中確有大面積的灌木和森林因洪水淹沒(méi)變?yōu)闈竦?。?jīng)人機(jī)交互實(shí)驗(yàn)，確定洪水發(fā)生后森林和灌木區(qū)域VV極化后向散射系數(shù)值增大3 dB以上的區(qū)域?yàn)樾略鰸竦貐^(qū)域，提取新增濕地面積約212.4 km2。

本研究較之國(guó)內(nèi)現(xiàn)有研究往往僅進(jìn)行新增開(kāi)放水體的提取，且較少考慮入射角對(duì)不同典型地物的后向散射系數(shù)影響的問(wèn)題，在利用ISODATA算法提取開(kāi)放水體的基礎(chǔ)上，分析了不同典型地物的后向散射系數(shù)隨局地入射角的變化情況和不同地物類型的紋理特征，進(jìn)行了洪水發(fā)生后新增濕地與灌木、森林的極化分解分析，利用閾值法提取了新增濕地面積，綜合新增開(kāi)放水體與新增濕地得到更為完整的洪水區(qū)面積。因本研究采用歷史衛(wèi)星遙感影像研究了國(guó)外的一次洪水過(guò)程，難以開(kāi)展實(shí)地災(zāi)情調(diào)查，故利用目視解譯分析評(píng)價(jià)本方法提取洪水區(qū)位置的合理性。下一步的研究中，擬利用本方法對(duì)國(guó)內(nèi)洪水過(guò)程開(kāi)展洪水區(qū)面積的實(shí)時(shí)提取，并開(kāi)展實(shí)地調(diào)查，進(jìn)一步評(píng)價(jià)本方法的精度與適用性。