基于Frangi滤波的长江中下游河流信息提取

2022-01-08李化良

城市勘测 2021年6期

李化良

(镇江市勘察测绘研究院，江苏镇江 212004)

1 引言

河流在全球水循环、工业热源和生物生长及人类生活需求中发挥了重要的作用[1，2]。随着全球气候的变化，河流水体的生态特征和生物、物理性质都出现了改变，影响了河水流动、泥沙传输和生物组成[3]。准确有效地提取河流信息对我们的生活和发展至关重要。随着遥感技术的不断发展，越来越多的中高分辨率遥感数据能够免费获取。这些数据在时间分辨率和空间分辨率上基本满足河流实时监测的需求，但是若要准确提取地表水体，需要高精度、鲁棒性强的河流算法。水体指数算法是最常见的河流提取方法，它根据水体光谱特征差异对目标水体特征进行增强，方法简单但精度相对较低[4～6]；光谱特征分类在河流水体中提取也很常见，例如监督分类和非监督分类算法，该算法容易忽略较少的水体[7，8]；深度学习算法在目前应用十分广泛，识别精度高但需要较多样本[9]。在利用遥感技术提取河流信息时，河流不同大小的尺度特征以及复杂的背景特征是影响河流提取精度的重要因素。基于以上两点，我们提出了一种基于Frangi滤波的河流信息提取算法[10]，考虑到近些年我国政府加大对长江流域的监管与保护，我们将研究区选取在长江中下游流域。

2 数据源及研究区

长江流域作为产业和城镇布局的重要载体长期以来都是我国重要经济区，准确、高效提取长江流域水体对于监测长江岸线和保护长江生态意义重大[11，12]。随着遥感技术的发展，近些年国内外提供了一系列分辨率较高、重访周期较短的卫星数据，如美国的Landsat系列卫星(Landsat 5、Landsat OLI)、欧空局的哨兵系列卫星(Sentinle-1，Sentinel-2A/B等)以及中国环境卫星(HJ)、资源卫星数据(ZY-1，ZY-3等)和高分系列数据(GF)。多源卫星数据为监测长江水体提供了丰富的基础数据源，通过高精度河流提取算法能够对遥感影像中的河流信息进行精确提取。为了监测全球气候变化对人类和环境影响，欧洲航空局实施了“伽利略”计划，目前共发射Sentinel-1，Sentinel-2A/B，Sentinle-3和Sentinel-5等一系列卫星。作为美国Landsat系列卫星的后续卫星，Sentinel-2A/B在空间定位和影像信噪比上相较于Landsat OLI都有较大改进[13]。Sentinel-2A/B共包括13个波段，分别为4个10m分辨率、6个 20 m分辨率和3个60 m分辨率波段。鉴于长江东西流向很长(约 6 300 km)，本文使用 10 m分辨率(波段2：490 nm，波段4：665 nm和波段8：842 nm)的Sentinel-2A/B影像来提取长江中下游区域。

3 基于Frangi滤波的河流信息提取算法

遥感影像识别地面河流水体时，河流形状和复杂的背景是影响算法提取的两个重要因素。因此，首先通过水体指数来初步增强河流水体特征；再使用Frangi滤波来进一步突出线状河流水体特征；最后利用水平集分割对河流进行提取，算法流程如图1所示。

图1 基于Frangi滤波的河流提取算法流程图

3.1 水体指数算法

水体指数是河流识别中常用方法，其原理为近红外或中红外波段水体的强吸收而导致的吸收谷，其简单有效且对面积较大的水体识别效果较好。目前，最常用的水体指数算法包括归一化差异水体指数(Normalized Difference Water Index)、改进的归一化差异水体指数(Modified NDWI，MNDWI)和增强水体指数(Enhanced Water Index，NWI)等[14]。本文利用Sentinel-2A/B MSI中3和8波进行水体指数运算：

(1)

其中band3和band8分别对应Sentinel-2A/B MSI影像的第3和第8波段。本文以江苏境内的长江流域为例来阐述河流信息提取的流程(图2(a))。经过波段运算，图2中的长江流域水体变得高亮显著。但从图2(b)中河流边界处水陆界线十分模糊，简单的阈值分割方法无法准确提取干流附近细小支流。

图2 (a)长江流域真彩色图(b)NDWI图

3.2 Frangi滤波算法

Frgangi滤波在医学影像中常常被用于视网膜、血管等线性几何结构的识别。视网膜和血管的形态与河流形状本质上是一致的，本文创新性地将Frangi滤波推广到陆地河流提取中。对原始影像做高斯运算，能够得到特定尺度下高斯滤波图像，公式如下：

Iσ(p)=Iσ(x，y)=I0(x，y)⊗Gσ(x，y)

(2)

其中，Iσ表示尺度为σ滤波处理后的图像，p=(x，y)为图像中像元位置，为卷积运算符，Gσ(x，y)代表标准偏差为σ的二维高斯核函数。其中，σ∈[σmin，σmax]，σmin和σmax分别表示最小和最大像元尺度值。

(3)

差分运算是识别数字图像中亮度变化差异明显邻域像元最常用算法，差异较大像元表现为物理性质变化、灰度特征不连续。不同尺度差分运算的基本形式由公式(3)表示：

(4)

利用二阶高斯滤波运算对水体指数图像进行处理，计算每一像元在特定尺度下高斯矩阵，为式(4)：

(5)

由于Hessian矩阵是对称矩阵，故存在两个实数特征值λ1和λ2及对应的特征向量e1和e2，特征值λ1、λ2表示像元灰度值变化大小和趋势。水体像元在图像中通常表现为高亮特征且背景低暗，此时λ1≈0且λ2<0[15]。

(6)

其中参数β，γ用来调节滤波敏感性的Rb和S。

根据尺度空间理论只有像元在最佳尺度因子时Lσ才能够达到最大值。因此，在多尺度线状特征提取过程中将不断进行多尺度运算来得到最大尺度特征值。

(7)

图3是不同尺度σ下线性河流的增强情况，可以发现σ=1时只能增强细小的水体(图3(a))；随着尺度的增加河流越来越明显清晰，在σ=6时已经能识别出来较宽的河流。由于长江太宽，使用σ=9无法对其进行完全识别，需要更大的尺度因子(图3(d))。为了节省计算时间，我们对水体指数图像(NDWI)取阈值选取较大阈值后直接提取明显河流，然后将得到的显著河流与Frangi滤波得到的结果相叠加，得到最终的图像。

图3 不同尺度σ下线状河流的增强情况

3.3 水平集分割

(8)

xN表示像元x的邻域，K(·)代表灵活窗口方程来保证临近像元非负性，b用来评估影像密度与真实情况反射率偏差。

为简化计算，该方程能够通过水平集转换进行区域离散化：

(9)

Mi(·)是聚类中心Ci类别方程，φ(·)为水平集分割算法。

对图3得到的结果做水平集分割，得到线状水体的分布，然后将所得结果与NDWI影像中提取的高亮河流合并，然后进行噪声去除得到最终结果。

4 河流提取结果及评价

图4为最终提取结果，图中不同尺度的河流都被提取出来，细小的河流也被有效地识别。宽度较宽的区域河流呈现高亮特征，细小的区域出现比较暗丝状结构，对最终得到的结果进行精度分析。进行精度分析时，对图2(a)真彩色图像中的水体进行目视解译做矢量化结果图，作为精度分析结果的“真值”，然后在ENVI软件中使用混淆矩阵对河流提取结果进行精度分析总体精度为93%，Kappa系数为0.86。与监督分类算法中支持向量机(Support vector machine，SVM)提取长江下游区域河流(图4(a))，本算法能够提取更为细小的河流和沟渠，河流连通性也更好。最终结果基本满足政府部分对农业灌溉、城镇居民用水和长江流域水体生态保护等方面需要。