1. 引言
在全球气候变化的大背景下,水灾的发生渐趋频繁,已是世界上对人类生存环境危胁最为严重的自然灾害之一。我国幅员辽阔,从北至南约有50度的纬度差异,受多种气候条件影响。特别是我国南部地区,常年受热带季风气候和亚热带季风气候的影响,夏季内陆地区温度上升迅速气压降低,海洋温度上升较慢相对内陆气压较高,暖湿空气受到压力作用向大陆运动,冷热空气对流易形成强降雨。强降雨和上游地区的冰雪消融造成的水灾每年对我国造成巨大的损失,因此导致的洪水灾害每年造成500亿以上的经济损失。其中1998年特大洪水灾害直接造成了1660亿的经济损失和4150人死亡,受灾人口达2.23亿。有效而快速的监测水灾的形成,确定水灾发生的区域,评估造成的损失对保障国民人身财产安全具有十分重要的意义。
由于强降雨的发生常伴随极端恶劣天气,成像条件较差,传统的光学传感器无法实时准确地获得云层下的信息,这就制约了光学遥感技术在水情监测中的应用。合成孔径雷达(SAR)具有全天时全天候对地成像的能力,不受气候条件的限制获取地物的后向散射信息,在水域范围的精确提取方面应用潜力巨大。目前在轨运行的合成孔径雷达卫星主要有中国发射的GF-3、美国宇航局(NASA)发射的Seasat-A、欧空局(ESA) ERS/Envisat/Sentinel系列卫星、意大利的COSMO-Skymed卫星、德国宇航局(DLR)的TerraSAR-X/TanDEM-X卫星、加拿大航天局(CAS)的RadarSAT卫星以及日本的JERS和ALOS卫星等。其中Sentinel-1是由欧空局发射的由两颗卫星组成的卫星星座,包括分别于2014年4月和2016年4月发射升空的Sentinel-1A和Sentinel-1B。卫星采用双极化模式对地表成像,两颗卫星协同工作可将重访周期缩短至6天,且数据免费对外开放,为世界各国及时开展水灾和各类地质灾害的应急监测提供了科学有效的数据支撑 [1] 。
目前,国内外对于极化SAR图像的分类的研究聚焦在目标分解问题,主要解决思路分“基于极化散射矩阵的分解”和“基于相干矩阵或协方差矩阵的分解”两大方向 [2] 。因研究对象具有时变复杂特征,故多采用相干矩阵和协方差矩阵实施目标分解 [3] 。如Cloude等于1997年提出的特征分解法 [4] 与Freeman提出的基于散射模型的三元分解法 [5] ,在分类的过程中均取得了较好的效果。当前针对双极化信息的分解研究尚不多见,国防科学技术大学的吴永辉等 [6] 提出的双极化修正方法在传统特征分解法的基础上对有效区域边界进行修正,在双极化信息分解的过程中具有较好的适应性。为充分利用Sentinel-1系列卫星具备双极化通道成像的优势,本文研究通过利用Sentinel-1系列卫星的双极化影像来提升水体提取和变化监测的精度。并针对之前利用极化信息的研究区域较为单一,结果可能存在特殊性的现象,进一步探究基于极化分解的分类结果在实际应用中可行性。
2. 基于Sentinel-1系列双极化数据的水域提取与变化监测方法
为准确提取双极化SAR数据中的水体目标,本文将基于Sentinel-1系列影像数据预处理得到的标准化后向散射信息,依次开展极化协方差矩阵的分解和基于H-A-Alpha值的非监督分类,提取不同时相的水域范围并获取其变化的客观信息 [7] 。上述基于Sentinel-1双极化SAR数据与变化监测的总体流程如图1所示。
2.1 影像预处理
影像预处理部分主要包括影像裁剪、辐射定标、滤波和多视。其中辐射定标是极化SAR系统相较于传统的单极化系统的预处理过程的最大区别,对单极化系统应用来说,辐射定标可能不是必要,但对于多极化SAR系统来说是必须的 [8] 。因为极化SAR系统要求不同极化通道的数据联合完成极化合成,不同极化通道的信息必须处于同一标准,才可以确保不同极化波之间的相干性。
受制于合成孔径雷达的成像机理,相干斑点噪声是其固有的原理性缺陷。SAR图像中的斑点噪声主要是由雷达波束照射到地面,经不同的散射体散射的后向散射信号之间相互干涉,合成的信号矢量相加结果有强有弱,反映在图像上为类似斑点的噪声。选择合适的滤波方法去除斑点噪声就显得尤为必要。精致极化Lee滤波(Refined Lee Filter)能够很好的解决极化协方差矩阵次对角线元素的噪声既不为乘性模型也不为加性模型的滤波问题 [9] [10] ,并且很好的保持极化信息。因此,本文选用精致极化Lee滤波对研究区域进行滤波处理。
Figure 1. Flow chart of water extraction from dual-polarization information
图1. 双极化信息提取水体变化流程图
2.2 基于双极化信息的极化分解技术
在极化数据的分析中主要针对分布式目标的后向散射特性分析,常将目标的后向散射矩阵表示为散射矢量
的形式。
(2.1)
其中S表示某一像元的后向散射系数矩阵,
下标
表示以极化方式j发射,极化方式i接收的后向散射系数。Trace表示矩阵的迹。
表示正交单位矩阵,根据常用的两组正交单位矩阵可将全极化数据散射矢量可以分别表示为:
(2.2)
(2.3)
将
分别与其共轭转置
做外积可以得到极化协方差矩阵和极化相干矩阵,我们通常使用这两种矩阵来描述散射过程。因为极化相干矩阵更容易解释其物理意义,所以在实际中更多采用的是极化相干矩阵 [11] [12] 。
Sentinel-1系统的数据只提供双极化信息,所以要对其散射矢量的表达做修正,基于VH/VV极化方式的后向散射系数矩阵可表示为
(2.4)
采用Lexicographic基和Pauli基分别对S矢量化,并基于总功率不变的前提化简得到:
(2.5)
(2.6)
由此得到极化协方差/相干矩阵为:
(2.7)
易得到:
(2.8)
上式中
分别表示协方差矩阵中的元素。
极化相干矩阵中包含充分的目标的极化散射信息,且与极化协方差矩阵也存在明显的对应关系。为进一步实现检测,分类和识别等目的,我们还需对这些数据做进一步的分析,建立极化协方差矩阵、极化相干矩阵与后向散射信号之间的联系。本文拟采用H-A-Alpha分解的结果基于wishart距离最短原则对研究区域进行分类。
分解在双极化分解中需对边界区域做修正 [6] [8] [13] [14] 。
(2.9)
将
平面有效边界区域修正为如下图2所示,由曲线1和曲线2包围,关于
对称。
区域Z1对应低熵面散射,包括平静水面,海冰面,非常光滑的陆地表面等。
区域Z2对应低熵偶极子散射,常来自于各向异性的植被。
区域Z3对应低熵二面角散射,常来自于孤立的电介质和含金属的二面角。
区域Z4对应中熵面散射。
区域Z5对应中熵偶极子散射,表现为被植被覆盖的地表散射。
区域Z6对应中熵二面角散射,在林业应用中表现为树冠对散射过程得影响。
区域Z7对应高熵多次散射。
区域Z8对应高熵二次散射,一般来源于植被。
其中第九个部分为高熵表面散射,表明目标的去极化效应很强,接近于随机噪声,基于散射机制的认知这个区域通常被认为是不可实现的区域。
2.3. 基于Wishart距离最短的非监督分类
非监督分类分类速度快,且无需任何先验知识,但分类的结果不能确定类别的属性。基于Wishart距离最短的非监督分类是计算每个聚类中心的相干矩阵:
(2.10)
其中
为类
中像素的个数。然后计算每个像素的相干矩阵
到各个类中心的Wishart距离,并重新分类:
(2.11)
将非监督分类与H-Alpha分解的结果结合可以将类别与物理散射机制相关联。当确定具体地物与当前环境下的散射的机制的联系,就可以得到非监督分类结果的属性。
Figure 2. H-Alpha effective area boundary and division
图2. H-Alpha有效区域边界及划分
2.4. 水域变化叠置分析
叠置分析是将两层或多层地图要素进行叠置产生一个新要素层的操作,其结果将原来要素分割成新的要素,新要素在综合了多层要素原有属性的基础上能够反映新的属性特征。本文在进行洪灾过后的水域变化检测时,首先通过对汛期前后两个时相的雷达影像分别进行面向对象分类,获得洪灾区域的水体部分,然后采用逻辑差运算的叠置分析方法,旨在获取洪灾后影像相对于灾前影像水域增加的位置和面积。
3. 鄱阳湖水域2018年汛期遥感监测与分析结果
3.1. 研究区域与实验数据
2018年7月3日长江防总启动防汛Ⅳ级应急响应,5日长江委水文局召开的防汛综合会发布消息称“长江第一号洪水”已经形成,江西省是本次水灾的重灾区。截止8日,灾害已造成江西省南昌,九江,抚州等9个设区市39县91.5万人受灾。本次洪水危害面积较广,十分典型。因此,本次实验将研究区域选定在江西省九江市周边范围研究区域(如图3所示)。
通过汛前6月25日和汛期7月7日两景Sentinel-1A IW成像模式雷达影像开展分析处理。提取汛期来临前和汛期的水域分布,用于开展变化监测。Sentinel-1搭载C波段SAR传感器,轨道高约700 km,工作频率5.4 GHz可采用单极化或双极化方式,IW成像模式只提供VH/VV极化方式成像。表1列举了Sentinel-1系列卫星的主要性能指标,并与同类卫星做了简要对比。
经与我国高分系列SAR卫星GF-3、加拿大新型商业卫星RadarSAT-2比较分析表明,Sentinel-1因为其较短的重访周期,高质量的成像能力在灾害发生时能够及时的获取可靠的灾区影像数据,用以研究水害灾情分析十分有利,可为科学的抗洪救灾提供可靠的依据。因此本文选用Sentinel-1的数据作为后续研究的数据源。为充分体现将双极化信息引入水体提取的效果,本文拟将最终实验结果与VV极化方式得到的监督分类结果进行比对。
Figure 3. Schematic diagram of the study area (the base image is lansdsat8 optical image)
图3. 研究区域示意图(底图为lansdsat8光学影像)
Table 1. Comparison of Sentinel-1 satellite and international similar satellite SAR load indicators
表1. Sentinel-1卫星与国际同类卫星SAR载荷指标对比
3.2. 实验结果分析
从原始的单极化强度图(图4a)来看,水体在图像中呈现暗色调与周边地物形成明显的差异,整体轮廓较为清晰,采用阈值分割的方法能够快速提取水体,其对应关系明确,物理解释也较为简单,但对表面光滑的地物和泥沙含量较高的水体不能有效判别。右图为Pauli假彩色图像,可以观察到在简单利用双极化信息后水体与非水体之间的边界更加容易识别,可进一步降低分析的难度。
为进一步利用极化反射特性,将图像根据极化协方差矩阵分解为熵值(H),平均反射角(Alpha),反熵(A)来分析水体在极化特性上的表现,如图5所示。从研究区域熵值图不难发现,水体在实验中熵值较正常水准略高,集中在中熵区域,部分出现在高熵区域。经过分析,我们认为造成水体熵值升高的主要原因是由于降雨和汛期水流加速引起的。对此我们查询了九江市6月25日和7月7日天气状况。6月25日阴天,7月7日小雨转中雨。比较了前后两景影像的熵值图(见图6)也反映了这一结论。6月25日熵值图如图左,蓝色框内水体熵值较低;而7月7日熵值图显示(如图右所示),水体熵值有明显增加。
水体的熵值偏高,为取得更好的结果,需要我们重新定义在该实验条件下,水体所对应的极化散射机制。对此,我们选择部分水体作为感兴趣区域来观察其在H-Alpha空间的分布,以此为依据来对初始的空间划分做一定程度的修正,使水体能够更好地被区分,其H-Alpha空间分布如图7所示。
结合水体在H-Alpha空间的分布,将整个区域划分为如图7(b)所示的8个区域,低熵区为0~0.5;中熵区为0.5~0.95;高熵区为0.95~1;表面散射区为0˚~40˚左右;偶极子散射区约为40˚~50˚;偶次散射区约为50˚~90˚;和反熵组成16个初始聚类中心,基于wishart距离最短原则进行分类,经过4次迭代计算和类别合并最终获取得到研究区域空间分辨率约为15米的水体提取结果如图8(b)所示。
将基于H-A-Alpha的非监督分类结果,与基于强度图的阈值分割提取水体的结果做分析比对。并采用混淆矩阵和Kappa系数做精度评估。在采用阈值分割得到的结果中Ⅰ区域水体未能被识别,但在基于
(a) VV通道SAR影像 (b) Pauli假彩色合成影像
Figure 4. Sentinel-1A image data imaging renderings
图4. Sentinel-1A影像数据成像效果图
(a) H值图 (b) Alpha值图 (c) Anisotropy
Figure 5. July 7 decomposition results grayscale
图5. 7月7日分解结果灰度图
(a) 6月25日熵值图 (b) 7月7日日熵值图
Figure 6. Entropy map of the main lake area of Poyang Lake
图6. 鄱阳湖主湖区熵值图
(a) 水体在H-Alpha空间的分布 (b) H-Alpha空间划分
Figure 7. 7 H-Alpha scatter plot
图7. H-Alpha散点图
(a) 阈值分割 (b) H-A-Alpha
Figure 8. Water extraction results in Poyang Lake area of Jiujiang City
图8. 九江市鄱阳湖区域水体提取结果
极化目标分解的结果中成功将其归类为水体。在II区域,可以明显看到,部分地物表现出较低的后向散射强度,从而使阈值分割将其识别为水体,而H-A-Alpha非监督分类的结果有效降低了错分的数量。因此引入极化信息可以有效的提升水体提取的精度。混淆矩阵和Kappa系数(见表2)也较好的证明了这一点。
显然,H-A-Alpha在水体提取精度上得到了一定程度的提升,识别的准确度达到了94.13%,较阈值提取结果有一定程度的额提升。对其他目标正确识别的准确度达到了99.67%,可以看到在防止其他地物错分方面较为突出,对于波浪较大和泥沙混杂的水域,能做到较好的的判识和保留相对于传统的阈值提取方法优势明显。但是,该结果只能与极化散射机制相对应,仍不能确定真实的地物类别。因真实的地物与极化散射机制的对应关系受外界环境以及雷达系统的影响较大,最终类别的确定仍需要一定先验知识。最后采用H-A-Alpha分类结果做为最终的水域提取结果。得到2018年6月25日水体面积约为1626平方公里,2018年7月7日水体面积为1682平方公里,可知受“长江第一号洪水”和江西本地强降雨的影响,九江市鄱阳湖区域水体面积增加迅速,12天内的新增水域面积为56平方公里见图9。
Table 2. Comparative analysis of water extraction results in the study area
表2. 研究区域水体提取结果精度对比
Figure 9. Water changes in the Poyang Lake area of Jiujiang City
图9. 九江市鄱阳湖区域水体变化
4. 结论
合成孔径雷达影像本身具有全天候观测的优势,本文为进一步挖掘极化信息在水体提取中的潜力,围绕极化信息的分解,应用以及解决实际问题的能力等问题开展研究。选取位于长江中下游的九江市鄱阳湖区域为研究区域开展实验。针对洪水爆发前6月25日和洪水高峰期7月7日两幅雷达影像,采用双极化分解结果作为遥感解译提取水体的基础研究数据开展分析。结合基于Wishart距离的非监督分类提取水体,旨在充分利用极化信息,降低阴影区域和泥沙含量较高对提取结果的干扰,并取得了良好的效果。随后,我们对汛期前后获取的雷达影像进行叠置分析,根据差异值评判洪灾期间的淹没区域和淹没面积。基于15米分辨率的水域提取结果,经变化检测得出整个九江市鄱阳湖区域在汛期的新增淹没面积达56平方公里。实验结果不仅论证了方法与流程的可行性,也充分反映出多极化雷达影像数据在洪灾应急监测应用中高分、精准的性能优势。
致谢
本文获得了国家重点研发计划“地球观测与导航”领域重点专项课题(2017YFB0502700),国家自然科学基金青年科学基金项目(41601503),四川省科技计划项目(2018JY0564),西南交通大学理工类科技创新项目(2682016CX087),西南交通大学“雏鹰学者”人才计划项目(2682016CY19)的联合资助。此外,欧洲空间局为本研究提供了Sentinel-1A卫星影像数据,在此一并表示感谢。