基于SAR图像桥梁多次散射的河流水位变化监测方法

丁泽刚 李莫凡 李 根 卫扬铠 胡子浛

（北京理工大学信息与电子学院雷达技术研究所，北京 100081）

1 引言

洪涝灾害是当今世界上最具破坏性的自然灾害之一，对国民经济和人民的生命财产安全产生了严重威胁。为了对洪涝灾害进行有效的管理，降低洪涝灾害的影响范围，有必要对河流进行持续的水位监测。

传统的河流水位监测方式是在河流周围设置水位测量仪器，包括超声波水位计［1-2］、雷达水位计［3-4］、激光水位计［5］等。这些水位计将传感器悬挂在水面上并发射信号脉冲，通过记录脉冲返回时间来计算水体表面的高度变化。然而，安装这些设施需要额外的人力和财力。在许多欠发达国家或交通不便的偏远地区，现场测量设备难以安装和维护［6］。

遥感技术是远程获取水文数据的一种有效手段，具有范围广、周期短、时效强且不受地面监测条件限制的特点，在洪涝灾害监测中发挥着极为重要的作用。近几十年来，人们一直在研究如何利用遥感技术监测河流［7-9］，随着卫星对地遥感技术的进步，出现了一些基于遥感图像的水位监测技术，主要分为两类：图像变化检测方法［10］和InSAR 方法［11-13］。第一种方法从光学图像或SAR图像中提取代表水体的像素，通过拟合水域面积与水位变化的数学关系，实现水位监测。这类方法的精度约为0.7 m［14］，其缺点是在拟合过程中水面高度变化的定量信息必须通过遥感图像以外的手段获得，如激光测高［15］或雷达测高［16］。第二种方法可以通过In-SAR 技术对特定区域的水位进行定量监测，该方法利用重轨卫星获取的湿地中非淹没水生植被区图像进行干涉处理，通过相位差计算卫星两次经过同一地点时的水位变化，其处理精度约为7 cm［17］。由于InSAR 对地物相干性的要求较高［18-19］，该方法只能应用于水面位置处分布有较强散射物体的场景，如含植被的湿地区域［20-22］或包含堤坝的水库区域［23-24］，目前其可应用场景有限。

基于桥梁的多次散射机理，本文提出了一种利用时间序列SAR 图像获取河流水位的方法，实现了桥下河流水位的定量监测，将遥感水位监测应用于架设桥梁的一般河流。同时，为了避免船只强散射信号的干扰，本文提出了一种桥梁能量累积算法，基于形态学的级联拟合处理，可提取SAR 图像中的桥梁方向并校正其在方位向的徙动，接着对桥梁散射信号沿方位向进行能量积累，大大提高了桥梁散射信号的信杂噪比。通过准确提取时序SAR 图像中的桥梁散射信号位置变化，可以反演得到河流水位变化。分析结果表明，该方法可以在目前的高分辨率星载时序干涉SAR 图像上实现分米级测量精度。最后，通过Sentinel-1A 数据和COSMO-SkyMed数据验证了所提方法的有效性，实现了分米级的估计精度。未来，随着星载SAR 分辨率的提高和重访时间的减少，该方法有望实现厘米级河流水位测量精度与数天级监测频率。

本文的后续章节安排如下：第2 节讨论了桥梁多次散射的机理和几何模型；第3 节介绍了基于能量累积算法的河流水位变化提取方法的处理流程；第4节中分析了本方法在星载SAR图像应用中的精度；第5 节进行了水上桥梁星载SAR 图像的处理与验证；第6 节讨论了本方法目前的处理能力和未来潜力；最后，第7节为结论。

2 基于桥梁能量累积算法的水位变化监测机理

SAR 信号脉冲的回波延迟决定了目标在SAR图像中的距离向坐标。来自同一目标的多次散射信号会由于时间延迟不同而出现在图像中的不同位置。桥梁在SAR 图像中通常表现出较高的散射强度和明显的线型特征，在河流区域能观察到明显的多径效应。

基于多次散射机理，来自桥梁的总后向散射信号（经算法聚焦后）可表示为：

其中：τ和η分别为距离向和方位向采样时刻；pr(t)和pa(t)分别为距离向和方位向的冲激函数；n=1，2，…，N是多次散射的次数（最高为N次）；A0为后向散射系数；R(n)为目标的n次散射在SAR 图像中的对应斜距；c为光速；f0为信号载频；fηc为非零多普勒中心。

多次散射几何如图1所示。其中θ为入射角，H为散射体的高度。根据图1 中的几何关系，来自同一目标的一次散射信号和二次散射信号之间的斜距差为Hcosθ。雷达发射脉冲在桥梁和水面之间还可以多次往返弹射，例如雷达脉冲沿着水面-桥-水面路径弹射形成三次散射，因此除一次散射和二次散射外，还可能会在图像中观察到相对低强度的高次散射信号。第n次散射的斜距可以表示为：

图1 多次散射路径示意图：一次散射和二次散射Fig.1 Sketches of signal paths：single-bounce（direct reflection）and double-bounce

其中Rb(η)为方位向η处桥梁的一次散射信号斜距，称为桥梁方位向徙动，由桥梁和SAR 卫星之间的空间位置关系决定。

假设桥梁的后向散射信号能够被准确识别，则在某一指定方位时刻η0时，来自桥梁的全部后向散射信号可以表示为：

因此，方位时刻η0时桥梁一次散射信号和二次散射信号的距离向采样时刻τ0和应分别表示为：

上述结果表明，多次散射的回波延时与桥梁距水面的高度相关。

对于一组配准后时间序列SAR 图像，不同SAR图像中特定桥梁的绝对高度（一次散射回波延时）可以视为常数，而对应的二次散射回波延时会由于河流水面高度变化而改变。因此，可以用以下方法来估计河流水位的变化：

在高分辨率的SAR 图像中可以观察到桥梁的细致结构以及其多次散射信号。为了详细讨论桥梁的多次散射现象，图2 展示了雷达脉冲在SAR 图像中不同散射路径的位置关系。其中凹形结构为桥梁在SAR 距离向与地面法线构成的平面中的截面，位于该截面上的4 个典型散射点Pi(i=1，2，3，4)用于说明多次散射机理：P1和P2为靠近雷达一侧的上下顶点，P3和P4为远离雷达一侧的上下顶点。Wi和分别为过Pi的信号脉冲在水面上的散射位置和Pi在水中的镜像位置。

图2 三种多次散射的斜距关系.红色为一次散射，橙色为二次散射，黄色为三次散射Fig.2 Relative location of three kinds of scattering in slant range geometry.Red：single-bounce；orange：double-bounce；yellow：triple-bounce

黑色箭头表示雷达脉冲的传播路径，θ为入射角，散射信号与其对应的路径起点用S标记，其中Si为Pi处的一次散射（直接后向散射返回传感器），S1、S2和S3在SAR 图像中的斜距值依次增大。是Pi处的二次散射，沿着桥-水面（或水面-桥）的散射路径返回到传感器。以为例，由于P′2是P2在水中的镜像位置，因此有：

其中O是P1（以及P2）在水面上的垂直投影，根据中位线定理可以进一步得到：

考虑到雷达的延迟是双程的，则Sdb和S′2两路径的差值为：

由公式（9）和公式（10）可知δτ=0，即两条路径Sdb和的时延相同。同样，可以证明Sdb和两条路径的时延是相同的，因此上述任何一种二次散射路径都可以等效为一次散射路径Sdb。可以看出，P2在SAR 图像中一次散射和二次散射的斜距差由P2到水面的高度Hb决定，并满足关系：

可采用公式（11）反演桥梁距离水面的高度。为Pi处的三次散射，根据水面-桥-水面的路径返回到传感器，由于：

因此上述两条传播路径可以分别等价于和处的一次散射。

水上桥梁的SAR幅度图像如图3所示。图像像素对应的斜距从左到右增大，彩色箭头展示了上述分析中三种桥梁散射信号在SAR 图像中的对应关系。红色标记的线性特征表示一次散射；橙色标记的二次散射在视觉上位于桥正下方的水面上，与先前分析吻合；用黄色标记的三次散射由于桥底部密集分布了强散射的材料，呈现出带状特征。

图3 SAR幅度图像Fig.3 SAR amplitude image of the bridge

3 基于桥梁能量累积算法的水位变化提取方法

本文基于自定义的桥梁多次散射识别算法，即能量累积算法（EAA），提出了一种水位变化提取方法。

图4介绍了利用时间序列SAR图像估计河流水位变化的核心思路。为简单起见，假设桥梁与SAR图像的方位轴平行（即桥梁方位向徙动Rb(η)为常数），第n幅SAR 图像的获取日期为tn，τRef和分别为tn处的一次散射和二次散射。通过计算不同时刻τRef与之间的斜距差值，可以利用公式（11）求出各时刻的桥梁距水面高度。目前星载SAR 的观测能够长期覆盖全球大部分地区，且重访周期较短，如Sentinel-1A 可小于12 天，因此该方法可以实现对河流水位变化的持续监测，也可以解译在水位监测任务开始前的历史水位变化。

图4 水位反演示意图Fig.4 The sketch of river level inversion

上述方法实施的前提是准确提取SAR 图像中桥梁散射信号的斜距信息。理想情况下，任意提取一行SAR 图像的数据都能观测到能量峰值，如公式（3）和图2 所示。但在实际应用中，由于散斑噪声、船只强散射干扰或水面波动等原因，仅从一行数据很难检测出各种桥梁散射回波的峰值能量，自动识别桥梁散射位置变得十分困难，如图5所示。

图5 特定方位时刻单行数据的散射能量分布示例Fig.5 Examples of energy distribution in azimuth line

图5（a）为桥梁区域SAR 图像，图5（b）为沿Row1提取的特定方位时刻单行幅度数据。在Row1中，桥梁的散射能量比船只的散射能量弱得多，难以提取一次散射和二次散射信号。注意到-τRef与方位采样时间η无关，可以将不同方位时刻的桥梁散射信号叠加以突出来自桥梁的能量，其计算方法为：

该表达式称为能量曲线，其中η∈[ηstart，ηend]为覆盖桥梁信号的方位采样时间，(η)为Rb(η)的估计。由于船只散射的能量并非沿着桥梁方向分布，因此在累积能量曲线的过程中，桥梁散射能量得到集中，而船只的能量被分散到不同位置，提高了桥梁散射信号的信杂噪比。本文基于形态学方法，实现了Rb(η)的精确估计，并对SAR 图像中的桥梁散射信号进行了累积。具体流程如图6所示。

图6 水位监测流程图Fig.6 The block scheme of river level monitoring

算法流程的输入是时间序列SAR 图像数据。在预处理步骤中，选择时间序列SAR 图像中与其他图像之间基线最短的图像作为主图像，并将其他图像配准到该单一主图像上。配准后所有图像都在距离向上进行升采样。然后取出幅度图像数据进行中值滤波，消除孤立噪声点的影响。幅度图像和中值滤波图像如图7（a）和（b）所示。

图7 算法流程的中间结果Fig.7 Illustrations of some steps

为了抑制船只散射和水面波动的干扰，本文采用EAA 对桥梁散射信号进行提取。EAA 的关键是准确提取桥梁方向，即实现Rb(η)的估计，此估计的准确度也受到船只强散射的影响。考虑到桥梁大多是直线型的，本文提出了一种基于形态学粗差消除的拟合方法，称为级联拟合（CF）：首先从中值滤波后的图像中搜索每个方位时刻中（即每一行）幅度最大的像素，将这些像素的值设置为1，将其他像素的值设置为0，新的二值图像称为最大值矩阵。然后对最大值矩阵进行形态学闭运算，使得桥梁区域的最大值像素相连，此时提取图像最大连通域中的最大值点，将此类点的方位距离坐标作为样本进行线性拟合，得到初始拟合结果。最后在拟合线上进行膨胀运算（膨胀尺度由样本方差决定），将位于膨胀区域的最大值点重置为拟合样本，得到更精确的拟合直线。不断按上述步骤迭代，直到没有新的最大值点加入，此时桥梁方向的估计结果即为最终结果。

根据CF 的处理结果，对SAR 图像中桥梁信号的方位向徙动进行矫正，使桥梁的散射能量垂直距离向分布。矫正结果如图7（c）所示。将校正后图像能量沿方位向累加，得到累积的能量曲线，如图7（d）所示。标有红色星号的峰值对应图3中的一次散射τ2，而标有橙色星号的峰值对应图3 中的二次散射。最后，计算两峰之间的距离，利用公式（11）得到各图像中水面与主图像的高度差，实现水位监测。

如图8（a）所示，应用本文提出的CF 算法后，对桥梁方向的估计相比采用CF算法前更准确，能够更稳健地提取散射信号的位置。将桥梁信号沿图8（a）中的拟合方向矫正，对应的能量曲线如图8（b）所示，随着迭代次数增加，能量曲线的峰值不断升高，说明桥梁方向估计得到不断优化。

图8 级联拟合的处理效果.红：首次拟合；蓝：一级级联拟合；绿：二级级联拟合Fig.8 Performance of cascade fitting.Red：first fitting；blue：CF epoch-1；green：CF epoch-2

4 精度分析

本文所提方法主要误差来源分为两类，一类是成像算法中由于多次散射现象造成的图像散焦，另一类是图像中桥梁散射信号的提取误差。主要误差构成如下图9所示。

图9 方法误差构成Fig.9 Error composition of the proposed method

本方法的图像散焦问题是由于多次散射现象导致的二次相位误差引起的。

为讨论该问题，建立图10 所示SAR 成像场景，该场景由SAR、桥梁侧沿和水面三部分组成。为简单起见，假设SAR 平台正侧视，飞行高度Hs并沿Y轴方向以速度Vr做匀速直线运动，桥梁侧沿和水面分别位于YZ平面和XY平面。如图10（a）所示，设方位时刻η=0 时，距水面高度Hb的桥梁散射目标A和SAR 平台达到最短斜距Ra=Ra0，此时A点在水面的垂直投影与其二次散射出射点分别为O和B，桥梁的二次散射等效单程斜距可以表示为：

图10 桥梁区域SAR成像场景Fig.10 SAR imaging for bridge scene

其中Rmb=Hb/cosθ为桥梁入射点A和水面出射点B之间的距离，Rb=Rb0为B点和SAR平台之间的距离。

由于SAR 成像处理是针对静止目标计算匹配滤波器调频率来实现方位向聚焦的，即认为A点和B点位置不随时间发生变化，则SAR 平台在方位时刻η收到A点的回波数据时，对应的斜距历程应为：

其中Rmb不随方位时刻变化，且考虑到Ra0和Rb0均远大于Vrη，因此对R(η)采取了抛物线近似。

但在SAR 平台运动过程中，SAR 平台位置与平面AOB不满足共面关系，A点的出射方向存在Y轴方向的分量，当偏转角度α≈Vrη/Ra0较小时，二次散射的水面出射点B′将近似沿Y轴方向匀速直线运动，如图10（b）所示。根据图中几何关系可知，B′沿SAR平台相反方向运动，且满足关系：

因此SAR 在采集A点回波数据时，真实斜距历程应为：

可见，在成像过程中基于静止目标假设设计的匹配滤波器将与实际的二次散射信号失配，调频率误差可以表示为：

由于调频率失配造成的二次相位误差将引起方位向旁瓣展宽，进而导致散焦信号在EAA 过程中积累，造成提取误差。根据二次相位误差和目标照射时间T的定义［25］可以计算对应的二次相位误差：

其中ρa和λ分别为雷达的方位向分辨率和工作波长。根据公式（20），在星载SAR 观测几何下，由二次相位误差引起的旁瓣展宽［25］不足1%，故可以被忽略。

本方法的提取误差是由于图像配准误差和取整误差引起的。

图像配准主要依赖高相干性的桥梁一次散射作为配准同名点，根据图4中的水位反演工作原理，当配准出现一定误差时不同图像中的τRef会发生偏移，进而影响水位的提取精度。图像配准精度通常要求达到1/8 像素的精度，因此处理误差与图像的距离分辨率有关，可表达为：

其中ρr为图像的距离分辨率。配准误差与ρr呈线性关系，斜率由1 cosθ决定。处理误差与距离方向分辨率的关系如图11（a）所示。例如，X 波段和L波段星载SAR 分辨率分别约为1 m 和6 m［26］，对应的配准误差分别约为0.12 m和0.70 m。

图11 提取误差分析Fig.11 Extraction error analysis

本方法通过CF 算法成功在距离升采样后的图像中精确提取了桥梁方向，并将图像的每行数据沿拟合桥梁方向进行搬移，完成徙动矫正。由于搬移时对徙动量进行了取整操作，导致每行搬移存在误差，进而导致多次散射信号的提取误差。该误差可表示为：

其中Nup为图像的距离向升采样倍数。为了保证在不同分辨率的图像中数据处理量接近，可以设定每个像素单元的斜距间距δp=ρr/Nup为定值，则取整误差是与δp和θ有关的一组曲线，如图11（b）所示。例如，将图像升采样至δp=0.1 m，对应的取整误差通常在0.15 m以下。

5 实验结果

为评估所提出方法的有效性，本节将详细介绍两个实测案例研究。案例1 和案例2 分别展示了使用COSMO-SkyMed 数据集和使用Sentinel-1A 数据集在中国洞庭湖的处理结果。

5.1 案例1：COSMO-SkyMed数据集

实验数据由COSMO-SkyMed 卫星在中国湖南洞庭湖获取，包括杭瑞高速桥和蒙华铁路桥两座水上桥梁。对应的光学图像和SAR 图像分别如图12和图13所示。数据的相关参数如表1所示。

表1 测试数据相关参数Tab.1 Characteristics of the test data

图12 关注区域的光学图像Fig.12 Optical image of the mentioned area

图13 COSMO-SkyMed获取的关注区域SAR图像Fig.13 SAR image of the mentioned area（from COSMO-SkyMed datasets）

数据集由2019 年1 月至2020 年12 月拍摄的N=24 张SAR 图像组成，每次观测间隔时间约为30 天。采用所提出的方法对两座桥下的水位进行估计，处理流程如图6所示。

数据集中的幅度图像与中值滤波结果如图7（a）和图7（b）所示，根据CF 算法的处理结果，对SAR 图像中桥的方位向徙动进行了校正，如图7（c）所示。图7（d）为图7（c）的能量曲线。将估计得到的河流水位与实地测量的水位真值进行比较，用于评估方法的性能。实地水位数据来自于岳阳当地的城陵矶水文站。此外，为便于比较，河流水位估计结果的参考高程应与实地水位数据的参考高程一致。

图14 展示了两座桥梁下的水位变化。红线表示使用本文方法估计的河流水位，记为，蓝线表示河流水位真值，记为，其中i=1，2…N表示第i次观测。从图中看出估计的水位变化曲线与实际水位变化曲线比较吻合。此外，2020 年7 月，在图像获取地点发生了由强降雨引起的洪涝灾害，洪水导致的水位上涨情况也被所提的方法捕捉到。

图14 案例1中的水位图，其中红色代表水位估计值，蓝色代表实际水位Fig.14 River levels of case 1.Red lines are predicted values whereas the blue lines represent ground truths

图15 展示了案例1 中估计结果的散点图，表示水位真值和估计水位之间的误差。蓝线表示估计值与真实值之间的线性回归结果，该结果与对角线（虚线）基本重叠，表示估计值非常接近真实值。

图15 估计结果的误差散点图Fig.15 Scatterplot of ground truth and estimated water level

为了量化分析实验结果，使用RMSE（均方根误差）和R（皮尔逊相关系数）对所提方法的性能进行评价，结果如表2所示。RMSE和R可表示为：

表2 处理精度指标Tab.2 Performance indicators

根据第4节的讨论，在上述情况下，本文所提方法的理论估计误差为：

表2 中的RMSE 结果接近此分析值，说明了本方法的有效性。剩余误差可能是由以下两个原因引起的：（1）水文站与桥梁之间有几公里的距离，影响了验证结果的准确性；（2）在SAR 图像获取过程中，由于风浪的影响，桥梁多次散射信号在SAR 图像中的位置不稳定。

5.2 案例2：Sentinel-1A数据集

实验数据由Sentinel-1A 卫星在中国湖南洞庭湖获取，包括杭瑞高速桥和蒙华铁路桥两座水上桥梁。对应的光学图像和SAR图像分别如图12和图16所示。数据的相关参数如表3所示。

表3 测试数据相关参数Tab.3 Characteristics of the test data

图16 Sentinel-1A获取的关注区域SAR图像Fig.16 SAR image of the mentioned area（from Sentinel-1A datasets）

数据集由2019 年1 月至2020 年12 月拍摄的N=42 张SAR 图像组成，每次观测间隔时间约为18 天。采用所提出的方法对两座桥下的水位进行估计，处理流程如图6所示。

数据集中的幅度图像与中值滤波结果如图17（a）和图17（b）所示，根据CF 算法的处理结果，对SAR图像中桥的方位向徙动进行了校正，如图17（c）所示。图17（d）为图17（c）的能量曲线。将估计得到的河流水位与实地测量的水位数据进行比较，实地水位数据来自于岳阳当地的城陵矶水文站。

图17 案例2的中间结果Fig.17 Results of some steps in case 2

图18 展示了两座桥梁下的水位变化。红线表示使用本文方法估计的河流水位，记为，蓝线表示河流水位真值，记为，其中i=1，2…N表示第i次观测。从图中看出估计的水位变化曲线与实际水位变化曲线比较吻合。

图18 案例2中的水位图，其中红色代表水位估计值，蓝色代表实际水位Fig.18 River levels of case 2.Red lines are predicted values whereas the blue lines represent ground truths

图19 展示了案例2 中估计结果的散点图，表示水位真值和估计水位之间的误差。蓝线表示估计值与真实值之间的线性回归结果，该结果与对角线（虚线）基本重叠，表示估计值非常接近真实值。

图19 估计结果的误差散点图Fig.19 Scatterplot of ground truth and estimated water level

使用RMSE（均方根误差）和R（皮尔逊相关系数）对所提方法的性能进行评价，结果如表4 所示。根据第4节的讨论，在上述情况下，本文所提方法的理论估计误差为：

表4 处理精度指标Tab.4 Performance indicators

表4 中的RMSE 结果接近此分析值，说明了本方法的有效性。

6 讨论

水位监测能力最重要的指标是测量精度和监测频率。根据第4 节的分析，本文方法测量精度的约束因素是星载SAR 的分辨率，而重访周期决定了监测频率。此外，低分辨率数据导致桥的不同散射信号相互重叠，这将从机理上降低了多次散射信号识别的精度。

一些在轨卫星的参数如表5所示。在轨卫星轨道高度均为500 km 以上，可以忽略二次散射信号的成像散焦影响；COSMO-SkyMed 可获得1 m 高分辨率InSAR图像，此时算法误差小于24 cm。基于以上分析，使用当前的星载时序干涉SAR 图像可以实现分米级监测。目前SAR卫星的重访周期主要在10～20天左右，这说明现有设备的监测频率约数周。

表5 部分在轨卫星参数Tab.5 Parameters of satellites in orbit

在未来，星载SAR 图像的距离分辨率可达到0.5 m～1 m，本文方法的精度可达厘米级。Capella SAR 星座［27］目前已经实现2.5 天的数据刷新率，这将进一步提高方法监测频率。因此，本方法有望利用先进的星载SAR 数据实现厘米级的水位测量精度和数天级监测频率，其精度与常见的现场水位监测设备相当。

7 结论

基于水上桥梁的多次散射机理，本文提出了一种利用时间序列SAR 图像监测河流水位变化的方法。船只强散射信号干扰是SAR 图像中桥梁多次散射信号识别的主要问题，本文提出了一种桥梁能量累积算法，基于形态学的级联拟合处理，可提取SAR 图像中的桥梁方向并校正其在方位向的徙动，接着对桥梁散射信号沿方位向进行能量积累，大大提高了桥梁散射信号的信杂噪比。通过准确提取时序SAR 图像中的桥梁散射信号位置变化，可以反演得到河流水位变化。精度分析表明，该方法可以在高分辨率星载时序干涉SAR 图像上可实现分米级的测量精度。COSMO-SkyMed与Sentinel-1A星载SAR 实测数据处理结果验证了所提方法的有效性。由于星载干涉SAR 图像易于获取，该方法有效拓展了遥感水位监测的应用场景，使全球更多河流能够得到及时监测。在未来，随着星载SAR 技术的发展，该方法有望实现厘米级水位测量精度与数天级监测频率。