测绘学 | 杨蕴：斑点抑制与多分辨率拓扑分析相结合的SAR图像河流水体提取

斑点抑制与多分辨率拓扑分析相结合的SAR图像河流水体提取

杨蕴,李玉,赵泉华

辽宁工程技术大学测绘与地理科学学院遥感科学与应用研究所, 辽宁 阜新 123000

基金项目：国家自然科学基金(42001286);辽宁省教育厅重点攻关项目(LJ2020ZD003)

关键词：SAR图像 河流水体提取 斑点抑制 局部阈值 凸包 多分辨率拓扑分析

引文格式：杨蕴, 李玉, 赵泉华. 斑点抑制与多分辨率拓扑分析相结合的SAR图像河流水体提取[J]. 测绘学 ，2022，51(1)：145-158. DOI: 10.11947/j.AGCS.2022.20190395

YANG Yun, LI Yu, ZHAO Quanhua. River waterbody extraction from SAR images based on speckle reduction and multi-resolution topological analysis[J]. Acta Geodaetica et Cartographica Sinica, 2022, 51(1): 145-158. DOI: 10.11947/j.AGCS.2022.20190395

阅读全文：
http://xb.sinomaps.com/article/2022/1001-1595/2022-1-145.htm

引 言

合成孔径雷达(synthetic aperture radar, SAR)以其全天时、全天候的对地观测能力，正广泛地应用于国防、国民经济建设以及防灾减灾等领域^[1]。利用SAR成像技术近实时、精确地获取河流水体信息对于水资源调查、自然灾害评估以及环境监测等具有重大意义^[2]。

由于水体表面相对平滑，其散射过程大多类似于镜向反射，因此在SAR成像过程中水体目标呈现出后向散射系数低、回波信号弱等特征^[3]，即在SAR图像中表现为暗色目标。在高分SAR图像中，河流水体在其宽度方向上通常占据若干像素，由此在几何上河流水体段表现为大小各异的多边形。综上，河流水体在高分SAR图像中表现为整体亮度较低且成弯曲分布的区域^[4-5]，这种鲜明特征使得图像分割成为对其快速识别和提取的重要途径之一。针对基于图像分割的SAR图像水体提取研究，目前已经提出如阈值、活动轮廓模型(active contour model, ACM)和边缘等方法。

基于阈值的分割法是针对SAR图像中水体具有的较低且相对一致的后向散射特性，通过选择低于给定阈值的像素来提取水体^[6-11]。文献[9]使用最大类间方差(Otsu)法设定阈值并据此提取水体，而为了降低斑点噪声造成的误提取，运用数学形态学的方法对提取结果进行后处理。文献[11]首先通过直方图统计对SAR图像进行后向散射强度分析，然后以Otsu和Kittler and Illingworth(KI)法对图像进行水体-非水体分类。试验结果表明在SAR图像水体提取中KI方法比Otsu方法更具有优势。基于阈值分割的SAR图像水体提取方法具有简单易行、性能稳定等优点，但SAR图像固有的斑点噪声严重影响其提取精度，因此该类方法需与滤波预处理或提取结果后处理的方法相结合^[12]。

ACM法包括参数ACM和几何ACM法^[13-20]。参数ACM法的主要思想是通过在水体周围定义一条初始轮廓线，在能量函数的约束下不断向水体边界演化，最终演化为提取水体的边界。文献[13]以阈值分割结果的边缘线作为初始轮廓线，由参数ACM将边缘段连接成连续的水体边界线。由于参数ACM无法处理轮廓线的自动分裂与合并，使得其无法适用于提取具有建筑遮挡的河流水体。与参数ACM法相比，几何ACM法需确定的参数较少且能够实现轮廓线的自动分裂与合并。如，文献[19]由几何ACM法在SAR图像中提取海岸线，其首先利用C-V模型进行粗分割来提供初始轮廓线，然后采用G0分布对轮廓线上每一点的邻域进行统计建模以实现海岸线精确检测，但该方法计算效率低，实现复杂。文献[20]提出了一种基于定位区域ACM(Localizing Region-based ACM, LRACM)的SAR图像水体提取方法，其用拉普拉斯核距离代替欧氏距离来计算区域能量，并结合局部和全局灰度值计算拟合中心和局部半径，从而更好地控制曲线的演化，但缺点在于演化曲线易与水陆边界附近的噪声点相连接，不利于后续处理。

边缘法是基于SAR图像水陆光谱特征差异较大这一特性利用边缘算子定位水体边界^[21-24]。如，文献[22]利用Prewitt算子结合轮廓因子进行边缘检测，以此提取河流水体。文献[23]将几种边缘检测算子(Sobel、Prewitt和Canny)应用于RadarSat-I SAR图像，并比较了它们检测水陆边界的能力，结果表明Canny算子具有最优的检测能力。对高分SAR图像，单纯采用基于边缘检测的方法难以提取水体的完整边界线，因此通常需要与其他方法相结合。如，文献[24]提出了一种结合指数加权均值比和水平集(ratio of exponentially weighted averages and level set, ROEWA-LS)的河流水体提取方法，其由ROEWA算子检测到的边缘为约束，利用ACM进行分割，结果表明该方法可较好地从高分SAR图像中提取河流水体。

虽然基于阈值、ACM、边缘及其相结合的方法在水体提取各具优点，但其均对SAR图像的斑点噪声非常敏感。故无论是在提取河流水体还是其他对象前，对SAR图像进行降噪预处理已成为提高提取精度的一个有效途径^[25-28]。SAR图像的斑点噪声与光学图像中被广泛研究的加性噪声不同，它属于乘性噪声且很大比例处于高频部分。针对SAR图像的斑点噪声抑制问题，目前已提出了如Lee、Kuan、Frost和SRAD(speckle reducing anisotropic diffusion, SRAD)等滤波器。虽然Lee、Kuan和Frost滤波器具有一定的降噪能力，但均存在对滤波窗口大小和形状敏感、边缘保持稳健性弱等缺点^[29]。针对此，一些学者提出了Lee、Kuan和Frost的增强版本^[30-32]及精致Lee和自适应精致Lee(adaptive refined Lee, AR-Lee)滤波器^[33-35]。在文献[35]中，AR-Lee的自适应表现为自适应窗口大小，它以子窗口建立中心和邻域窗口，根据边缘检测和相邻窗口相干系数区分同质和异质区域，从而选择合适的窗口大小。该方法较固定窗口大小的精致Lee，可较好抑制斑点噪声，但其仍受窗口大小范围的限制。G-MAP假定图像的概率分布函数服从Gamma分布，利用该先验信息对图像斑点噪声进行抑制^[36]。该方法克服窗口的限制，但先验的Gamma分布对不同场景不具有稳健性。相比之下，由于SRAD滤波器采用各向异性扩散的滤波方式，不受滤波窗口的影响，而且通过结合一阶和二阶梯度算子来控制边缘两侧的扩散，可在滤波过程中有效保持边缘^[37]。

1 研究方法

1.1 河流水体的SAR影像特征总结

影像目标特征是反映目标局部区域的光谱辐射特性与反映目标之间空间和几何关系的有机结合。为了更好地提取河流水体，必须了解其在高分SAR图像中的影像特征，由文献[41]的描述对其总结如下。

(1) 光谱特征：整体呈现暗色，内部光谱强度较均匀，与其相邻区域灰度反差较大。

(2) 几何形态特征：显示为弯曲的多边形，具有一定的长度，长宽比较大，宽度变化较小，方向变化较慢。

(3) 拓扑特征：具有连通性，相互间有间隙，不连续。

(4) 功能特征：通常与跨河建筑物(例如桥梁和渡槽)相伴随。

1.2 ASRAD斑点噪声滤波

设I={Ii,j, 1≤i≤M, 1≤j≤N}为原始SAR强度图像，其中，(i,j)表示像素位置；Ii,j∈{0, …, 255}为(i,j)处像素的灰度值；M和N分别为图像的行和列数。对任一像素点(i,j)，考虑其一阶邻域系统，利用前向差分来近似时间导数，可得SRAD模型的数值近似为

(1)

(2)

式中，

(3)

qi,j为瞬时扩散系数，可通过以下公式计算

(4)

q₀(ζ)为平滑尺度函数，可近似表示为

在ASRAD中引入β指标作为定量标准来控制其迭代过程。β指标利用原始图像和滤波图像的边缘相似性来评价滤波过程中边缘的保存比率，即如果边缘在滤波过程中得到很好的保留，则β指标将接近于1，定义为

(6)

式中，ΔIn和ΔIf分别表示原始图像和滤波图像的边缘图像；

式中，ε为控制迭代次数的阈值，ε越小，图像平滑度越大，信息损失越多，根据经验取ε=0.2时可得到较好保持高强度反差边缘的降噪图像。

1.3 滤波后图像河流水体段提取

1.3.1 河流水体段粗提取

基于水体在SAR图像中的光谱特征，可假设图像中只有水体和背景两类地物，同时考虑到F的局部不均匀性，利用Wolf局部阈值算法快速分割出河流水体区域。

对F内任一位置(i,j)的像素点，记以(i,j)为中心的滑动窗口大小为w×w，通过计算局部阈值Ti,j对该位置进行水体或背景判断

(8)

式中，Ei,j∈{0, 1}，当Ei,j=0时表示背景，Ei,j=1时表示水体；Ti,j是在以(i,j)为中心的w×w窗口内计算的局部阈值，可表示为

(10)

(11)

遍历所有像素，得到河流水体的粗提取结果E={Ei,j, 1≤i≤M, 1≤j≤N}。

1.3.2 伪河流水体段去除

由于在F中河流水体与沥青路面、污水池等地物具有相似的光谱特征，使得E中存在一些零星的伪河流水体，为此根据2.1节的总结，设计精细化处理对其进行去除。首先搜索E中的连通区域，对其进行标识，统计每个标识块的面积Az及其外接最小矩形的长Lz和宽Wz，结合图像分辨率和工程实践的需要设定两个阈值，若其都在阈值区间，即

图 1 伪河流水体去除算法流程 Fig. 1 Flow chart of pseudo river waterbody removal algorithm

图选项

1.4 多分辨率拓扑分析的河流水体段自动连接

由河流水体在高分SAR图像中的拓扑特征和功能特征，在提取的河流水体段之间可能有一定距离。为了得到完整的河流水体提取结果，需将E′中间断的相邻河流水体段相连接。

1.4.1 河流水体段的最小凸包计算

图 2 凸包计算流程 Fig. 2 Flow chart of convex hull calculation

图选项

依次处理E′中的各个连通区域，得到对应的凸包集P={P(z),z=1, 2，…,Z′},P(z)为E′中第z个河流水体段计算的凸包。

1.4.2 多分辨率拓扑分析的间断点定位

利用构建图像GP来获得在不同分辨率情况下河流水体拓扑结构的差异，从而自动定位河流水体段之间的间断点。

以P为输入，通过对其进行高斯平滑和子抽样来构建GP，可表示为

(13)

式中，Gk表示GP的第k层，k∈{1, 2，…,K}，

为构建GP的层数；logS(·)表示以S为底的对数；

表示下取整运算符，Hl是用于降低分辨率的二维高斯低通滤波器；?表示卷积运算；S为步长，通常取S=2；(↓S)[·]表示下采样步长为S。

通过相邻GP层之间执行加法来获得连接金字塔(connection pyramid, CP)，当k=K时

(14)

当k<K时

(15)

式中，Ck是第k层CP的表示；C′k是对Ck的扩展，以匹配k-1层上GP的Gk-1的大小；Hh是用于分辨率扩展的二维高通滤波器；(↑S)[·]表示上采样步长为S。

由上面对原始输入P以及GP的定义可知，构造的Ck中只存在0，1和2这3种像素强度值。由于原始输入P只有0和1两种像素强度值，其中等于1和0分别表示河流水体和背景，则可以做出以下分析。当Ck中任一位置(i,j)处的强度值为0时，表示相邻两层在此位置皆为背景，无须处理；当Ck中任一位置(i,j)处的强度值为2时，表示相邻两层在此位置皆为河流水体，无须处理；当Ck中任一位置(i,j)处的强度值为1时，由高斯金字塔的定义可知下层和上层在此位置的值分别为0和1，表示下层和上层分别为河流水体和背景，则此位置可能位于河流水体段之间的间断，即需处理的点。通过遍历Ck，k∈{1, 2，…,K}，可自动获得每层可能属于间断的像素点。

1.4.3 间断点生长的河流水体段连接

对所有可能的间断像素点，以其为种子点，采用区域生长策略对其真实属性进行判别分类，以此来实现河流水体段之间间断的连接。

随着k增大，Gk的分辨率逐步降低，间断也越窄，因此从CP的上层到下层搜索河道间隙的半径应逐步增大。由Ck的构造过程，可知任一种子点在其中的最大搜索半径可表示为

其中，S为构造GP的采样步长。

对Ck，k∈{1, 2，…,K}中对任一种子点，以其为中心沿8个角度向4个方向生长，如图 3所示。4个方向分别为水平(1?5)，右对角线(2?6)，垂直(3?7)和左对角线(4?8)。种子点的生长策略是从种子点开始，沿着每个生长角度搜索，直到找到强度值为2的像素，并且在搜索路径上没有强度为0的像素。以像素为单位记录每个角度的搜索路径长度rko，其中o∈{1, …, 8}是生长角度索引。

图 3 S=2，k=1的种子的搜索方向和半径 Fig. 3 Search directions and radius for seeds withS=2,k=1

图选项

统计种子点在4个方向的搜索结果，以确定种子点的属性。一旦4个方向中的任何一个满足间断条件，将此种子点与河道段连接。以水平方向为例，当水平方向的两个搜索半径的和rk¹+rk⁵≤Rk，可将此点分类为河流水体段之间的间断点，其连接过程是将间断点变为河流水体段点，即Ck(i,j)=1→Ck(i,j)=2。遍历Ck内可能的间断点，完成此层的间断连接并将连接结果传递到下一层以作为下层间断点的识别依据，由此逐层完成间断的连接。

当k=1时，完成C₀内间断的连接，记为C′₀，由此得到最终河流水体提取结果J=C′₀–P+E′={Ji,j, 1≤i≤M, 1≤j≤N}，其中Ji,j∈{0, 1}，Ji,j=1表示河流水体，Ji,j=0表示背景。考虑可能存在一些形态上与真实河流水体段相似的伪河流水体，使得在1.3.2节的伪河流水体去除中未能去除。为此，当完成河流水体段连接后，可再次使用1.3.2节的方法，利用较大的面积阈值对连接结果中可能存在的伪河流水体进行去除。

2 试验结果与分析

2.1 试验设置

通过MATLAB R2018b编写程序进行提取试验，试验运行环境为Windows10 32位专业版操作系统，处理器为Intel(R)Core(TM)CPU 32 G的个人计算机。

图 4 原始SAR强度图像 Fig. 4 Original SAR intensity images

图选项

为了验证提出方法的稳健性，在GF-3卫星上获取的SAR图像上选取多幅较大尺度的城镇河流水体场景图像，限于篇幅，选择具有最大尺度的两幅图像进行展示和分析，其分辨率为3 m、超精细条带成像模式、大小分别为2000×6000像素和4400×7700像素。

2.2 ASRAD滤波结果与分析

图 5 斑点滤波结果对比 Fig. 5 Comparison of speckle filtering results

图选项

由图 5的对比可以看出，G-MAP滤波法对斑点噪声有一定的抑制作用，但对不同场景的滤波效果不一致。AR-Lee法由于能够根据边缘检测信息和邻域相干性高低程度调整滑动窗口大小，在斑点噪声抑制和边缘信息保持方面都得到了提高，但在如图 5(b)中的非水体区域仍存在大量的杂散点，这严重影响后续处理。SRAD法对同质区域内部噪声取得较好的抑制作用，但在固定滤波迭代次数下，如图 5(g)中的一些弱强度河流水体边界被破坏，相比之下，ASRAD自适应地调节滤波次数，可在最大限度降噪的同时保持水陆边界信息，有助于河流水体的更精确提取。对于滤波结果的定量评价，选取MSE，SNR，PSNR和式(6)的β度量。MSE值越低，表示降噪图像越接近原始图像，算法具有更好的性能，SNR和PSNR反映真实与估计图像误差之间的关系，其值越大表示效果越好。各评价指标见表 1，由表中可以看出，ASRAD法具有较小的MSE值，较大的SNR和PSNR值，说明它能较好地滤除噪声，且具有较大的β值，说明在滤除噪声的同时也能有效地保持边缘。

表 1 滤波结果的定量评价 Tab. 1 Quantitative evaluation on filtering results

图像	方法	MSE	SNR	PSNR	β
图 4(a)	G-MAP	0.385	65.144	64.253	0.695
	AR-Lee	0.418	63.546	61.436	0.634
	SRAD	0.325	67.375	67.598	0.714
	ASRAD	0.316	66.168	66.040	0.746
图 4(b)	G-MAP	0.457	66.818	65.024	0.614
	AR-Lee	0.393	68.462	65.863	0.642
	SRAD	0.335	74.311	71.971	0.691
	ASRAD	0.297	77.462	76.143	0.735
图 4(c)	G-MAP	0.311	67.659	64.984	0.624
	AR-Lee	0.362	68.646	66.473	0.647
	SRAD	0.313	72.354	70.436	0.695
	ASRAD	0.284	76.456	74.789	0.722

表选项

2.3 滤波图像河流水体段提取结果与分析

利用Wolf局部阈值分割算法对滤波后的图像进行河流水体段粗提取，结果如图 6(a)-(c)所示。由于存在同谱异质现象，使得粗提取结果中存在一些零星的伪河流水体，为此使用连通区域标识法对其精细化处理，结果如图 6(d)-(f)所示。