强降雨下地表汇水的时空分布信息提取方法研究

2022-03-29葛艳琴李彦荣

无线电工程 2022年3期

葛艳琴,李彦荣

(太原理工大学矿业工程学院，山西太原 030024)

0 引言

中国黄土高原气候相对干旱、降雨集中、雨热同期，晋陕甘宁是中国黄土分布的中心区域，约占全国黄土总面积的72%。黄土地区因具有土结构疏松、沟壑纵横、地形破碎以及气候干燥等特点，其地质环境极为脆弱，在降雨集中时易造成水土流失和生态环境恶化[1-2]。流域是气象降雨的地表水流汇集地，强降雨条件下流域周边衍生的落水洞加剧了黄土地区的地貌侵蚀问题，易发生坍塌、滑坡等黄土地质灾害，因而需探明降雨导致的地表汇水对流域周边地貌侵蚀的影响。研究该影响的首要任务是高效提取大范围地表汇水的时空分布信息，鉴于遥感对地观测技术在监测地表水环境的传统优势，亟需利用水体提取方法开展地表汇水的遥感监测研究。

在区域发生强降雨(24 h降雨量超50 mm)现象时，地表径流极易造成短时的局部汇水，汇水形成的地表水体分布较为离散，低分辨率遥感影像难以准确地识别其边界。水体信息遥感提取方法中应用最为广泛的是阈值分析法。单通道的阈值分析依据水体在近红外和中红外光谱通道上显著的吸收特性[3]，但无法有效区分水体与山体阴影[4]；基于多谱段关系分析和水体指数的多通道阈值分析是利用水体在不同光谱通道上的组合关系或计算模型，得到更高解译质量的特征影像。多谱段关系分析能够应用于地形起伏地小的地区[5-6]，与单通道法结合时可应用于山体阴影的去除[7]。水体指数法的代表性工作是McFeeters等[8]提出的归一化差异水体指数，既可消除建筑物和裸土影响，也部分抑制了植被、阴影等的地表弱反射特征。徐涵秋[9]提出了适用于城市水体提取的改进归一化差异水体指数，闫霈等[10]提出了修订归一化差异水体指数可抑制山体阴影和混合像元问题，沈占锋等[11]提出基于高斯转换的归一化差异水体指数来保存水体提取形态完整性。增强水体指数[12]、新型水体指数[13]、更具普适性的自动水体提取指数[14]和阴影水体指数[15]等方法也被发展用于特定场景的水体提取问题或削弱阴影、建筑物等噪声的干扰。

基于分类思想的水体提取方法中，决策树、支持向量机(Support Vector Machines,SVM)、面向对象是基于学习的主要方法。监督分类器[16]能够去除水体掩膜中的高反射特征地物，可利用相邻景影像重叠区域[17]或进行自组织训练[18]拓展样本。决策树分类是由上至下逐级细分得到水体像元，可与IHS(Intensity Hue Saturation)变换[19]、密度聚类、水体指数等进行组合分析，解决样本噪声与属性缺失问题[20]。SVM法在水体类型标记样本量很少的情况下能保持好的分类质量[21]，在水体分布差异复杂度高时，比指数分析法和面向对象法的提取精度更高[22]。面向对象方法可解决提取结果破碎化与混合像元问题，如结合水体指数与多尺度分割的方法[23]、结合光谱特征与拓扑关系的方法[24]等。此外，深度学习方法充分利用影像中水体的空谱信息，如深度学习与SVM分类[25]、多尺度融合[26]等方法可联合进行水体提取，并可改进从水体视觉特征到语义特征的自下向上的识别结构[27]。

现有研究表明，SAR传感器数据在水体信息遥感提取方面具有显著优势，但应用瓶颈在于传统水体提取方法难以自适应确定水体与非水体像元的光谱分割阈值。有鉴于此，结合分贝化与种子点扩散思想，提出了一种基于SAR影像的自适应水体提取方法，该方法通过SAR影像直方图分析能够自动化确定水体分割阈值，可用于高效监测黄土地区流域周边地表汇水的时空分布信息。

1 研究区与数据源

实验选取的研究区位于汾河下游新绛段，属暖温带大陆性气候，地理坐标为东经111°11′32″～111°14′37″、北纬35°35′11″～35°37′48″，平均海拔约392 m，分布有黄土丘陵植被、平原植被和汾河沼泽植被等多类型植被区。区内河道为汾河高低阶地构成的冲/湖积平原区，河道整体呈东西向展布，长约6 000 m、宽约200 m，年平均流量48.7 m3/s，年径流量约15.4亿m3，南北两侧为农田和新绛县城区。

为监测强降雨条件下区内地表汇水的变化过程，选取了高分一号(GF-1)多光谱数据、高分三号(GF-3)和哨兵一号(Sentinel-1)C波段SAR数据以及无人机航摄多光谱数据来进行水体信息的遥感提取与验证，其中SAR数据用于地表汇水信息的自动化提取，卫星与航空多光谱数据作为参考数据源用于对部分提取结果进行精度评价。所采用的实验数据源及其元数据信息如表1所示。

表1 研究区实验数据源

在利用Sentinel-1和GF-3 SAR卫星数据进行水体提取前，需要进行初步的数据预处理，主要包括多视处理、滤波处理、辐射定标以及地理编码。首先将斜距复数数据转换为SAR强度数据，弱化相干斑的影响，增强2种SAR卫星影像的解译特征；然后利用Refined Lee滤波器对SAR强度数据进行滤波，消除其中的斑点噪声，提升图像信噪比；Sentinel-1和GF-3影像中通常存在显著的辐射偏差，需借助辐射定标处理将SAR影像强度值转换为地距后向散射系数值；最后，为消除因SAR卫星侧视成像造成其对地形起伏较为敏感且易产生几何畸变的问题，应对SAR影像进行几何校正，以消除透视收缩、叠掩和阴影等问题，即地理编码处理。

采用GF-1卫星和无人机航空摄影2种方式对2021年10月12日(对应无人机航摄数据)和17日(对应GF-1卫星数据)的Sentinel-1 SAR卫星影像水体提取结果进行验证。GF-1卫星全色-多光谱(PAN-Multispectral，PMS)数据在完成几何与辐射处理后，采用基于Gram-Schmidt的Pansharpening方法将全色波段(2 m)和多光谱波段(4 m)融合成2 m分辨率的多光谱影像。无人机的航摄载荷为多光谱相机，采集流程包括新绛县河道周边的测区踏勘、航线设计、像控点布设以及现场航空摄影测量等步骤，像控点信息需利用RTK测量方法采集，采集到的航空相片空间分辨率优于10 cm；然后，对采集到的原始无人机影像依次进行几何畸变纠正、特征提取与匹配、空三加密、影像镶嵌等处理，生成研究区的多光谱数字正射影像图(Digital Orthophoto Map，DOM)。此外，本研究中涉及的所有多源多模态实验数据均需进行几何纠正与影像配准，并下采样到统一的像元尺寸下，便于后续水体提取结果的对比分析与验证。

2 研究方法

监测强降雨条件下汾河下游新绛段地表汇水时空分布信息的整体思路为：先将经预处理后的多时相SAR影像进行分贝化处理，然后根据直方图分析法确定水体与非水体的阈值分割点，再利用阈值分割结果作为水体种子点构建区域生长模型，最后基于区域生长算法实现对多时相SAR数据的水体提取，得到研究区地表汇水的时空分布信息。具体提取方案如图1所示。

图1 本文发展的流域洪灾汇水信息提取方案Fig.1 Proposed scheme of water body extraction

2.1 基于智能种子点扩散的地表汇水信息提取

在利用种子点扩散方法进行水体提取前，对SAR影像进行分贝化处理(即对数变换)，使其能够对SAR影像中的水体类别进行特征增强，提高其解译标志的显著性与可判别性。分贝化处理可表达为：

dB(σ0)=10×lgσ0，

(1)

式中，σ0为SAR卫星数据的同极化后向散射系数；dB(·)为分贝化后的分贝值。SAR影像经分贝化处理后，其数据分布呈现出明显的高斯分布特征。

水体种子点确定是区域生长算法提取水体的起始步骤，它对于水体区域的生长至关重要。传统的种子点确定方法是基于人机交互解译的手动选取方式，该方法执行效率低且生长结果受主观因素影响，不适用于区域级的水体提取研究。本文通过分析SAR影像dB值的统计分布特征，提出了一种结合直方图分析与区域种子点扩散的自动化水体提取方法，该方法利用直方图分析方法来确定水体种子像元的分割阈值，在筛选出水体种子区域的基础上利用种子点扩散方法得到最终的水体分布掩膜图。

首先对经分贝化的dB值SAR影像进行直方图分析，研究区时序SAR影像dB值的直方图分布既存在单峰分布也存在双峰结构。研究区部分时相下SAR影像单峰和双峰直方图结构的分析图。如图2和图3所示。

图2 基于单峰直方图分析的SAR影像水体种子像元提取Fig.2 Seed extraction of water body in SAR images based on unimodal histogram analysis

图3 基于双峰直方图分析的SAR影像水体种子像元提取Fig.3 Seed extraction of water body in SAR images based on bimodal histogram analysis

对于单峰直方图分布情形(图2)，在dB值符合高斯正态分布的前提下，其两侧拐点指示了显著的地物类别区分特征，根据水体的dB值分布特点，选取其左侧拐点作为水体目标与非水体目标的阈值分割点，该位置所对应的dB值通常为影像标准差(STD)的负值，可表达为-STD，则水体种子像元的判别依据可表达为：

dB(σ0)≤-STD。

(2)

对于SAR影像dB值呈双峰结构的直方图分布情形(图3)，其dB值分布范围定义为[σ1,σn]，采用双峰法分析其直方图分布特征，即通过计算2个波峰(对应图3中σi和σj)间峰谷σt的位置来确定阈值分割的基本依据。具体实现途径是：① 搜索直方图中最大峰值σj；② 直方图平滑/滤波处理；③ 统计最大峰值左侧的其他峰值；④ 基于其他各峰值及其与最大峰值距离找出第2个波峰σi；⑤ 根据最大和第2个峰值计算两峰之间的峰谷值σt。则水体种子像元可根据式(3)判别：

dB(σ0)≤σt。

(3)

在利用直方图分析方法确定出水体种子点像元的分割阈值后，可得到初步的水体“种子区域”掩膜图，并以筛选出的水体种子像元作为区域生长的起始生长点。区域生长是从水体“种子区域”对应的灰度矩阵中多个种子点开始，当水体种子点像元Pi与其邻近像元Pj间满足如下条件时，则将像元Pj加入到水体种子点集合中：

|dB(Pi)-dB(Pj)|≤K，

(4)

式中，K为区域生长预设的生长阈值。

在首次扩散得到扩展后的水体种子点集合后，按照更新后的种子点集合重新进行区域生长，重复上述生长过程直到遍历所有SAR影像像元，最终筛选出所有符合条件的水体类别像元。

2.2 地表汇水信息提取精度评价

为精确检验2.1中所述方法的水体分布信息提取质量，主要采用基于人工目视解译的方法从光学遥感影像(GF-1 PMS和无人机多光谱数据)中精确提取水体分布边界，作为提取质量评价的标准参考掩膜。在对水体目标进行目视提取前，可利用归一化水体指数(Normalized Difference Water Index，NDWI)将卫星或无人机多光谱遥感数据转换为水体专题图，并通过阈值分析等方法作为多光谱数据的解译参考。某多光谱像元的水体指数NDWI的生成算法表示为：

(5)

式中，ρgreen和ρNIR分别为经辐射定标与大气校正处理后绿光波段、近红外波段的像元反射率，由于无人机成像时距离地面高度较小，因此在航摄数据采集时只需利用定标板进行基础的辐射定标处理即可。

在对本文算法的水体提取结果与标准参考掩膜进行质量评价时，采用总体一致性(精度)po、漏检率、误检率以及Kappa系数等质量评价指标进行综合评估。其中，总体一致性po是待处理图像中识别出的与参考掩膜一致的正确水体像元占图像总像元数量的比值，漏检率是待处理图像中未识别出的参考掩膜中的水体像元占图像总像元数量的比值，误检率是待处理图像中识别错误的水体像元占图像总像元数量的比值。Kappa系数的计算式为：

(6)

(7)

式中，Wa和Wr分别为标准参考掩膜中的水体类别像元个数以及其中算法识别出的水体像元个数；NWa和NWr分别为标准参考掩膜中的非水体类别像元个数以及其中算法识别出的非水体像元个数；n为水体掩膜图像中的像元总数。

3 实验结果与分析

实验采用所发展的基于智能种子点扩散的水体信息提取方法，将依次采集于2021年9月30日、10月5日、10月12日、10月17日的Sentinel-1时间序列SAR卫星影像以及10月10日下的GF-3卫星SAR影像全部作为研究区地表汇水信息提取实验的基础研究数据，利用采集于10月12日的无人机航摄多光谱影像和10月17日下的GF-1 PMS影像进行地表汇水分布信息的人工目视解译，目视解译的结果用于对10月12日和17日的Sentinel-1卫星SAR影像水体提取结果进行精度检验。此外，上述所有实验数据均已经过相应的数据预处理，包括几何配准与重采样操作。

图4为2021年9月30日、10月5日、10月10日、10月12日以及10月17日共5个时相下Sentinel-1或GF-3卫星的分贝化影像、水体种子点(区域)以及水体提取结果，其中分贝化影像(图4(a)，(d)，(g)，(j)，(m))上叠加的红色实线为水体提取结果掩膜的矢量边界;图5展示了研究区10月12日的UAV航摄DOM产品、10月17日GF-1 PMS融合影像及其目视解译结果，该结果将作为参考的精细水体提取结果来评价图4中的SAR卫星影像水体提取掩膜(图4(c)，(f)，(i)，(l)，(o))；图6为10月12日和17日的SAR影像水体提取结果与参考水体掩膜的对比分析图，包括本文水体提取结果在2个时相下的正确检测掩膜(图6(a)，(d))、漏检测掩膜(图6(b)，(e))、误检测掩膜(图6(c)，(f))以及这3类掩膜在UAV多光谱影像(图6(g))和高分一号多光谱影像(图6(h))上的叠加展示图；其中，图6(g)，(h)中标注的红色实线为对应时相下SAR影像正确提取的水体掩膜边界、绿色实线为漏提取的水体掩膜边界、蓝色实线为误提取的水体掩膜边界。

表2为研究区在2021年10月12日和17日下本文算法提取结果与参考水体掩膜的质量评价指标，包括总体检测精度、误检率、漏检率以及Kappa系数等；图7为图4中所示的5个时相(2021年9月30日、10月5日、10月10日、10月12日以及10月17日)下基于SAR卫星影像提取的水体分布面积变化曲线，它直接展示了汾河下游新绛县段在强降雨条件下地表汇水信息的时序变化及影响范围。其中，10月5日相对于9月30日的地表汇水面积增加了0.32 km2，10月10日相对于10月5日的地表汇水面积增加了2.3 km2，10月12日相对于10月5日的地表汇水面积减少了0.21 km2，10月17日相对于10月12日的水体提取面积进一步减少了1.57 km2。图8为研究区气象地基雷达的降雨强度扫描图，指示了区内的降雨时点和降雨强度，其中降雨强度以组合反射率(单位dBz)形式呈现，红色星标注区域为研究区地理位置。

表2 研究区水体提取结果质量评价指标

(a) 9月30日Sentinel-1影像

(b)9月30日水体种子点

(e)10月5日水体种子点

(f)10月5日水体提取结果

(g)10月10日GF-3影像

(h)10月10日水体种子点

(i)10月10日水体提取结果

(j)10月12日Sentinel-1影像

(k)10月12日水体种子点

(l)10月12日水体提取结果

(m)10月17日Sentinel-1影像

(n)10月17日水体种子点

(a)10月12日UAV多光谱影像

(b)10月12日水体目视解译结果

(d)10月17日水体目视解译结果

(a)10月12日正确检测掩膜

(b)10月12日漏检测掩膜

(d)10月17日正确检测掩膜

(e)10月17日漏检测掩膜

(f)10月17日误检测掩膜

(g)10月12日多掩膜叠加

(h)10月17日多掩膜叠加

图7 基于SAR影像提取的地表汇水信息变化曲线Fig.7 Variation curve of surface catchment distribution derived from SAR images

(a) 9-30 21:48

(b) 10-3 21:30

(d) 10-5 17:36

(e) 10-6 15:00

(f) 10-7 21:12

图4中时序SAR影像的地表汇水提取结果表明：本文发展的基于智能种子点扩散的水体提取方法在种子点生成过程中既能够较为全面地捕捉研究区内明显的水体种子像元，这些水体种子像元在接下来的区域生长构建中连通了成团簇的水体分布区域，实现了黄土地区流域周边地表汇水分布信息的精确监测，同时避免了显著的区域漏检与误检问题。另一方面，由于Sentinel-1卫星影像在成像几何上的时序辐射偏差及其与GF-3卫星影像在成像几何与空间分辨率(Sentinel-1分辨率为15 m，GF-3分辨率为1 m)上的差异性，5个时相下的SAR卫星影像在水体提取结果上也存在少量的不一致性，但实验结果也表明该问题对水体提取精度无显著影响。

对比图4中10月12日、17日下Sentinel-1 SAR影像水体提取结果与10月12日UAV航摄影像及10月17日GF-1多光谱影像的水体目视解译参考结果可发现，在检测精度上，SAR影像的总体水体提取精度优于87%，与GF-1光学卫星数据解译结果的总体一致性达90%以上，这表明光学卫星(分辨率为2 m)得到的水体目视解译结果更多地被SAR卫星影像捕捉到；但SAR影像水体提取结果与无人机航摄影像得到的水体掩膜的Kappa系数显著高于其与GF-1光学影像，这说明前二者在水体和非水体识别上的整体一致性更好。此外，SAR影像水体解译结果与无人机航摄影像解译结果间的漏检率更高(>10%)，则主要是因SAR卫星数据(分辨率为15 m)与参考数据源(UAV分辨率为0.1 m、GF-1 PMS分辨率为1 m)的空间尺度差异程度不同所引起的。在水体类别像元的形态分布特征上，基于人工目视解译的水体掩膜在聚集性和完整性上具有天然优势，但难以精细刻画水体在复杂地形地物条件下的分布轮廓，这主要受解译过程中的专家经验与解译时效性(效率)的影响；基于本文算法的水体提取结果能够检测出参考水体掩膜中的绝大部分水体像元，水体分布区域的边界特征不够平滑，这一方面是由于区域生长过程的算法结构造成，但另一方面更为重要的是由于作为水体提取实验的SAR卫星数据源空间分辨率大大低于作为参考水体掩膜的无人机和光学卫星高分辨率遥感影像，而利用人工目视解译方法从SAR影像中提取水体的复杂性也高于人眼观测的光学影像。因此，不考虑研究区数据采集条件的情况，综合这些多源、多模态、多尺度的遥感数据源进行地表汇水分布信息监测，其本身就存在系统观测误差。

在上述质量评价及不确定性分析的基础上，通过对研究区9月30日—10月17日间5个观测时相下SAR卫星影像的水体智能提取结果分析表明，研究区地表汇水面积从9月30日(约0.78 km2)起整体呈现出先增加后减少的分布特征。灾害前期的地表汇水面积呈加速递增的发展趋势，在10月10日的SAR影像提取得到的汇水面积达到峰值(约3.41 km2)，此后汇水面积急速减少至1.63 km2(10月17日)，这表明强降雨下黄土地区流域周边地表汇水具有短时、高汇流特点。对比图8的气象地基雷达扫描图可发现，研究区主要降雨时点均发生在10月8日之前，期间的持续强降雨造成了10月10日的地表汇水峰值现象。