河流表面流场是水流内部湍动结构的外在体现，其分布特征与河流中的超大尺度流向涡^[1]、二次流^[2]等相干结构息息相关；同时，表面流场也是河流水文测量的重要参数，直接或间接控制着径流过程^[3]、污染物扩散过程^[4]、侵蚀过程^[5]和泥沙运动过程^[6]等。开展河流大范围表面流场测量有助于解译明渠水流相干结构，为研究河流水体中的生物运动或物质运输研究提供基础数据。

河流表面流场会受到插入式测流仪器本身的干扰，目前主要采用非接触式测量方法，包括电磁波测速法和光学测速法^[7]。电磁波测速法以侧扫雷达测流系统^[8]为代表，通过计算电磁波的Doppler频移测量水面短波速度以表征表面流速。侧扫雷达测流系统可以开展大范围流场测量，空间分辨率最高约10 m，最小河宽要求为30 m，设备价格昂贵且现场安装布置环境要求高。光学测速法以大尺度粒子图像测速(LSPIV)技术^[9]为代表，利用摄像机采集一段河流表面示踪体运动视频，通过图像互相关运算计算连续帧图像中示踪微团的位移来表征河流表面流速，已广泛用于室内河工模型试验^[9-10]、野外河流监测^[11]和水工建筑物过流监测^[12]等，实践表明其室内测流精度约5%^[13]，野外测流精度约10%^[14]。LSPIV常采用固定式倾斜摄影模式，适合长时段无人值守测量，但存在图像中不同像素元的空间分辨率不同和测量范围有限的不足，无法实现长河段高空间分辨率的表面流场测量。

针对天然河道大范围表面流场测量，多名研究者利用无人机在高空正视拍摄河流^[14-16]，测量河流表面流场，并提出了针对性的现场测流导则^[17]。Cao等^[14]采用无人机悬停模式采集顺直河道视频影像，通过河岸特征点提取与匹配，计算原始视频中的水面流速并转换至世界坐标系，实现悬停模式快速测流。Detert等^[15]采用无人机分段悬停模式采集弯曲河道视频影像，并基于地面控制点校正视频抖动，采用粒子图像测速(PIV)技术分段计算时均流场并拼接为全流场。Fujita等^[16]采用无人机巡航拍摄顺直河道，通过相位分析法确定兴趣区域，提取并匹配该区域的尺度不变特征校正巡航影像，并采用时空图像测速算法计算表面流场。

目前研究局限于空中悬停拍摄模式或短距离航拍模式，无法开展长距离(>300 m)巡航测量；同时，应用场景局限于开阔河道环境，假定采集影像为二维平面，采用投影矩阵进行图像校正，不适合河道周边环境存在明显高程变化时的影像正射校正。长距离复杂环境巡航测量的难点在于每帧图像对应的摄像机位置不同，需逐帧进行相机外参数标定以开展正射校正；此外，相同空间分辨率条件下长距离测量生成的全景正射图像尺寸较大，流场计算耗时长。

针对上述不足，本文基于运动恢复结构(SfM)^[18]方法和粒子图像测速技术，搭建长河段高空间分辨率表面流场测量系统，利用视觉算法解算各帧图像与参考坐标系之间的转换关系实现多帧图像的正射校正，并对连续帧局部影像进行PIV计算，叠加至统一流场网格进行数据后处理与可视化，实现长距离高空间分辨率的全河段流场重构。该系统成功用于水电站引水渠和城市河流表面流场测量。

1 系统组成

无人机图像测速(airborne image velocimetry, AIV)系统^[19]的硬件主要包括无人机及云台相机、地面控制点(GCPs)和示踪粒子，如图 1所示。无人机飞行至测量河段上方沿流向巡航并拍摄含水面示踪粒子及河岸信息的影像。基于地面控制点坐标，对影像数据进行正射校正并计算瞬时表面流场，最终拼接成为全景表面流场。

地面控制点布置于测量河段两侧，用于测量区域三维场景重构和图像正射校正，地面控制点可选用人工标记物或地面物体并精确测量其三维空间坐标。河流表面天然示踪体(表面纹理、泡漩、落叶等)稀疏且受光照及风速等环境影响较大，无法用于高空间分辨率的流场测量。本文开展河流大范围表面流场测量时，采用2~4 cm尺寸的方形薄木块或淀粉泡沫作为示踪体，此类材料对水流跟随性强，可有效表征水流的流动，且绿色环保对水体无害。

航拍过程中，无人机采用正视成像方式采集水面示踪粒子图像。相比倾斜成像模式，无人机正视成像方式采集的水面图像像素空间分辨率基本一致，可有效减小倾斜摄影带来的图像畸变和正射校正时像素插值引起的图像失真，降低流场计算的不确定性^[20]。

AIV系统的误差主要包括图像正射校正误差和流场计算误差。由于视频中提取图像不包括地理信息，图像正射校正精度仅取决于GCPs的数量、分布和实测坐标精度。受环境光和障碍物干扰，局部区域瞬时流场存在错误矢量，进而影响局部区域测流精度，需进行人工剔错。

2 技术原理

无人机巡航测量长河段表面流场技术的核心是无人机图像的批量正射校正和复杂环境条件下的流场计算，技术流程大致如图 2，具体实现步骤包括：

1) 无人机沿河道巡航时录制示踪粒子运动视频，从视频中提取时间序列的图像对，作为原始输入数据；

2) 采用基于SfM和多视角立体视觉(MVS)算法的无人机巡航图像批量正射校正方法估算各时刻的相机外参数，根据投影转换关系，将拍摄图像投影至河流水平面，批量生成正射图像；

3) 采用PIV算法计算局部河段瞬时表面流场，拼接得到全河段表面流场。

2.1 无人机巡航图像批量正射校正

无人机沿河道采集影像时，受航线偏移、高空风速和飞手操作的影响，飞行姿态沿程变化，采集的视频图像出现几何畸变。本文基于Agisoft PhotoScan Professional软件平台，采用SfM和MVS算法，利用多视角图像恢复出相机的运动参数和目标的结构信息^[21]，实现长序列无人机图像的批量正射校正。

首先，采用间隔s帧提取连续2帧图像的方法从视频图像矩阵I_i(i=1, 2, …, n)中提取测流图像矩阵序列In_j，

$ \left\{\begin{array}{l} \mathbf{I n}_{2 k-1}=\boldsymbol{I}_{(s+2)(k-1)+1}, \\ \mathbf{I n}_{2 k}=\boldsymbol{I}_{(s+2)(k-1)+2}, \end{array}\right. $

(1)

$ k=1, 2, \cdots, \text { Floor }\left(\frac{n}{s+2}\right), $

(2)

$ s=(1-\beta) \frac{L}{U_{\text {drone }}} f. $

(3)

其中：n为原始视频对应的总帧数，Floor为向下取整函数，β是图像重叠率，L为无人机单画幅拍摄区域沿航线方向的实际长度，U_drone为无人机平均飞行速度，f为视频帧率。

利用SfM计算出摄像机坐标系下每帧图像对应的相对摄像机矩阵，并结合GCPs坐标及其在MVS生成稠密点云中的坐标，将相对摄像机矩阵转换为世界坐标系下的绝对摄像机矩阵P。

根据单目视觉原理，定义空间点为X=[X, Y, Z, 1]^T，对应的像平面成像点坐标为x=[x, y, 1]^T，可建立投影方程：

$ k_r \cdot \boldsymbol{x}=\boldsymbol{P} \cdot \boldsymbol{X}. $

(4)

其中：k_r为比例系数，P为3×4维矩阵。该方程可进一步展开为

$ k_r\left[\begin{array}{l} x \\ y \\ 1 \end{array}\right]=\left[\begin{array}{llll} p_{11} & p_{12} & p_{13} & p_{14} \\ p_{21} & p_{22} & p_{23} & p_{24} \\ p_{31} & p_{32} & p_{33} & p_{34} \end{array}\right]\left[\begin{array}{c} X \\ Y \\ Z \\ 1 \end{array}\right]. $

(5)

其中：[x, y]和p为已知数，[X, Y, Z]和k_r为未知数，该方程含义为图像矩阵里的每一个像素坐标都对应于世界坐标系下的一条空间直线。

进一步，基于MVS生成的世界坐标系下的稠密点云可以生成水面与河岸的数字高程模型(DEM)，是测量河段内规则网格点的平面坐标[X, Y]与高程Z的数据集。根据式(5)确定的空间直线，遍历DEM数据集，可确定距离空间直线最近的点[X, Y, Z]，即像素坐标[x, y]对应的世界坐标[X, Y, Z]。根据所有像素点的(X, Y)坐标及颜色值，定义正射图像横纵轴方向空间分辨率λ_x和λ_y，可插值出正射投影图像矩阵Io_j。Io_j的地理信息G_j=[X_ul, Y_ul, λ_x, λ_y]，其中[X_ul, Y_ul]为图像左上角像素点的地理坐标。

2.2 表面流场计算

不同摄像机位置生成的正射影像尺寸各异，无法直接进行流场计算。大范围表面流场计算首先对连续帧图像进行尺寸补齐操作，然后计算流场并转化为参考坐标系下的流速矢量。具体过程包括图像对预处理、PIV计算和流场数据后处理。

2.2.1 图像对预处理

根据连续两帧图像矩阵Io_2k-1和Io_2k的地理信息G_2k-1和G_2k，进行像素补零，得到相同尺寸的图像矩阵Iu_2k-1和Iu_2k并保存新图像的地理信息Gn_k(X_uln, Y_uln, λ_x, λ_y)。

图像对统一尺寸后，需进行图像增强预处理，通过增强示踪粒子信号，减少背景干扰，获得更好的流场计算效果，具体的图像处理过程为：1) 使用二维Gauss滤波器对原始图像矩阵Iu_j进行卷积，生成背景图像矩阵Ib_j；2) 原始图像矩阵Iu_j减去背景图像矩阵Ib_j得到前景图像矩阵If_j；3) 前景图像If_j进行灰度拉伸，得到增强图像矩阵Ip_j。

2.2.2 PIV计算

PIV是一种可以无干扰测量流体瞬时速度场的测量技术，已广泛应用于各种尺度的流场测量。PIV算法将连续帧粒子图像划分为均匀的矩形判读窗口，选取某一位置判读窗口，与下一帧图像中对应位置附近的判读窗口逐一进行互相关运算，相关系数最大值的位置相对于窗口中心点的距离和方向即代表该判读窗口内流体微团位移的大小和方向^[22]。

假设连续2帧图像中某一位置的判读窗口尺寸为M×N，对应的窗口灰度函数为g₁(i, j)和g₂(i, j)，定义(ξ, η)为判读窗口内流体微团运动可能产生的位移量(-M≤ξ≤M，-N≤η≤N)，通过式(6)计算所有可能的位移量对应的互相关系数，最大互相关系数对应的位移量即为判读窗口中心点的位移。

$ R(\xi, \eta)=\frac{\sum\limits_{i=-M}^M \sum\limits_{j=-N}^N g_1(i, j) g_2(i+\xi, j+\eta)}{\sqrt{\sum\limits_{i=-M}^M \sum\limits_{j=-N}^N g_1{ }^2(i, j) \sum\limits_{i=-M}^M \sum\limits_{j=-N}^N g_2{ }^2(i, j)}}. $

(6)

本文采用PIVLab^[23]软件，基于连续帧增强图像矩阵Ip_2k-1和Ip_2k，采用多重网格迭代的PIV算法计算像素尺度下的第k帧位移场U_k(x, y, Δx, Δy)。

2.2.3 流场数据后处理

流场数据后处理主要包括流场数据的坐标转换、剔错和拼接。首先根据Gn_k(X_uln, Y_uln, λ_x, λ_y)，将U_k转换为世界坐标系下的瞬时流场Un_k(X, Y, u, v)：

$ \left\{\begin{array}{l} X=X_{\mathrm{uln}}+x / \lambda_x, \\ Y=Y_{\mathrm{uln}}+y / \lambda_y, \\ u=f \Delta x / \lambda_x, \\ v=f \Delta y / \lambda_y. \end{array}\right. $

(7)

进一步采用最大最小流速阈值法和归一化中值法^[24]剔除Un_k中的异常值。最后，对整个测量区域划分地理网格，将每帧流场转换或插值至统一的网格，得到全景瞬时流场Up_k。对所有时刻的Up_k进行时间平均，得到测量区域的全河段流场为Um(X, Y, U, V)，流速矢量的大小为U_bulk。

$ U_{\text {bulk }}=\sqrt{U^2+V^2}. $

(8)

3 应用实践 3.1 水电站引水渠

该系统应用于瑞士苏黎世州巴登市利马特河的Schiffmühle水电站(47.488 °N, 8.269 °E)引水渠及溢洪道，河宽约100 m，其中水电站引水渠为长约500 m的弯曲河段，如图 3。测量期间，位于水电站上游5 km处的水文站，根据水位流量关系推算断面流量约71.5 m³/s。

3.1.1 测量布置

测量区域沿岸布置有26个GCPs，并采用实时差分定位技术(RTK)测量空间坐标。试验中采用玉米淀粉泡沫作为示踪粒子，沿程桥梁、拐点和水面共布置6个示踪粒子投掷点，确保示踪粒子沿程分布均匀，以实现良好的流场覆盖度。施测现场布置如图 3。测量时，无人机沿引水渠飞行并拍摄水面示踪粒子运动视频(分辨率为4 096×2 160像素，帧率为23.98 Hz)。

为评估流场测量精度，在引水渠入口及溢洪道入口分别布设声学Doppler剖面流速仪(ADCP)，型号为River Pro 1 200 kHz，同步测量断面流速分布。

3.1.2 正射图像

基于原始视频图像序列，设置相邻图像对间的重叠率为β=80%，共提取579对典型图像对；选定空间分辨率为λ_x=λ_y=64像素/m，定义研究河段矩形区域的西南角点为坐标系原点，采用2.1节的图像批量正射校正方法逐帧输出正射图像。GCPs平均投影误差为±0.21像素，说明正射校正精度高。沿程选取4个位置的原始图像如图 4a，经过处理后输出的正射影像如图 4b，其叠加至全景图像的效果如图 4c，结果表明，正射图像与研究区域全景图吻合良好。

3.1.3 表面流场

PIV计算流场时的最终采样窗口为32×32像素(0.5 m×0.5 m)。根据2.2节的步骤处理得到579帧流场，计算全河段表面流场如图 5，其流场分布和流线走向准确刻画了引水渠弯道水流特征。

图 5表明，弯道段水流受离心力影响，在溢流道闸门与引水渠交点处形成水流分离区，主流逐渐偏向凹岸；同时，引水渠前半段过水断面沿程逐渐收缩，导致断面平均流速和水面流速沿程增加；水流出弯后，进入顺直段，表面流速分布呈抛物线型，表面主流流速大于1 m/s，局部最大流速达2 m/s；水流进入引水渠末端，主流下潜至水轮机组，表层水体受建筑物阻滞，表面流速逐渐减小。

受光照不均匀性和两岸乔木干扰，测量段的5个区域出现了不合理的流速矢量。A₁区域由于正射图像失真，导致流场散乱；A₂区域几乎无示踪粒子覆盖，导致流场缺失；B区域邻近引水渠的观测台，近岸水面下方布置有建筑物，导致此处流态复杂；C区域生长有一高耸乔木，其相对于摄像机的运动干扰了流场处理；D区域水浅流急，水面产生明显的波纹，导致图像出现强反光区域，影响图像测速计算。

3.1.4 比测结果

由于水面以下0.4 m是ADCP测量盲区，本节选取ADCP各垂直测线距离水面最近测点的流向速度作为参考值。溢洪道断面A-A′与引水渠断面B-B′的位置如图 6a(位于图 5中的红色虚线框)，其中蓝色箭头为AIV实测断面流速矢量。2个比测断面的AIV与ADCP实测结果对比如图 6b，因ADCP数据点离散程度高，此处采取10点滑动平均法平滑处理。

比测结果表明，AIV与ADCP数据吻合良好，变化趋势基本一致，76%的测点相对误差小于10%；引水渠进口断面处，AIV流速比ADCP实测流速整体偏大，与宽浅明渠水流表面流速(AIV实测值)大于水下流速(ADCP实测值)的规律一致；溢洪道出口窄深断面处，与ADCP流速相比，左侧主流区AIV流速趋势基本一致，右侧整体偏小，与右侧出现回流有关。

3.2 城市河流

该系统应用于苏黎世城区的Unterhard河段(47.391 °N, 8.525 °E)，其左岸设有水文测站，如图 7。Unterhard河段长约350 m，河宽48 m，断面呈U形，为典型顺直河型。AIV施测期间，水文站利用跨河缆道搭载旋桨流速仪和测距探头分别测量断面流速分布和断面形态，统计分析得到比测期间断面平均水深为1.68 m，断面过水流量约为96.8 m³/s，断面平均流速约为1.2 m/s。

3.2.1 测量布置

河道两侧布置有14个GCPs，并通过本地GIS系统提取空间坐标。试验选用4 cm×4 cm的薄木片作为示踪粒子，从上游人行天桥进行移动投掷，覆盖全河宽范围。试验期间，无人机在距离水面47 m高度沿流向飞行，采集示踪粒子沿程运动视频(分辨率为4 096×2 160像素，帧率为23.98 Hz)，飞行时长5 min。

3.2.2 表面流场

数据处理阶段，设置相邻图像对间隔帧数为14帧，共提取654对图像并计算得到654帧流场。系统输出的tiff格式时均表面流场如图 8，其中坐标系原点为研究区域西南角点；流场直接叠加于苏黎世州正射图像(https://maps.zh.ch/)上，表明图像正射校正质量高。全河段流场覆盖度高，流线顺直，大致与河岸平行。测量河段近桥墩区域(图中黑色虚线框)的河床上存在一根大型乔木，其作用类似斜向潜坝，导致乔木以上过水面积减小，水面流速增大。

3.2.3 比测结果

比测断面位于如图 8中的A-A′测线位置，选取旋桨流速仪垂向最高测点(水面以下0.12 m)流向速度为参考数据。受水文站跨河缆道影响，A-A′测线处的AIV流场测量存在错误矢量，无法直接用于比测，因此选取A-A′测线上下游方向分别间距1和2 m处总共4条测线流向速度均值作为AIV在比测断面处的测量值，比测结果如图 9。结果表明，比测断面处AIV与旋桨流速仪实测流向速度趋势吻合良好，呈抛物线形。根据AIV各测点的流速数据，插值出位于旋桨流速仪测点位置的AIV流速值，其相对于旋桨流速仪实测值的误差大致位于±5%范围，全部测点相对误差均方根值约4.3%。

4 结论

本文针对自然河流长河段表面流场测量，搭建了巡航模式下的无人机图像测速系统。

1) 提出了长河段高空间分辨率表面流场测量的解决方案，采用运动恢复结构算法和地面控制点坐标，解算各时刻摄像机矩阵，实现无人机巡航图像的批量正射校正；通过局部流场计算与拼接，得到全河段表面流场。

2) 应用于水电站引水渠弯曲段和城市顺直河段，得到了高覆盖度和高空间分辨率(0.5 m×0.5 m)的全河段表面流场，真实反映出弯曲河道水流分离特征和顺直河道流向速度抛物线分布特征并准确刻画局部地形对水流的影响，相比ADCP和旋桨流速仪实测值的相对误差分别为10%和5%，可以满足野外测流需求。

3) 该系统适用于水面示踪粒子均匀覆盖条件下的长河段表面流场测量，施测时应选择合适天气和时间，减少水面强反光对图像的干扰。

本文采用示踪粒子来表征河流水面流动状态，但天然河流水面存在高低速流带和泡漩等大尺度湍流结构，稀疏的示踪粒子作为漂浮物会更大概率集中于高速流带，示踪粒子运动无法涵盖全部水面湍流结构，可能引入测量误差，有待进一步研究。