基于无人机低空遥感和OTSU算法的黑龙江开江流凌密度研究

汪恩良，胡胜博，韩红卫，刘承前

（1.东北农业大学 水利与土木工程学院，黑龙江 哈尔滨 150030；2.黑龙江省寒区水资源与水利工程重点实验室，黑龙江 哈尔滨 150030）

1 研究背景

冬季封河和春季开河阶段，我国北纬30°以北的河流普遍存在冰情危害［1］，随着河流断面冰凌分布密度的增大，冰凌撞击水工建筑物的概率增大，且大量流凌在特定河段，如河道束窄处、浅滩、急弯、连续弯道、尚未解冻冰盖的前缘等处受阻，形成冰凌阻水体，往往严重阻塞过水断面，使上游水位显著壅高，造成凌汛灾害［2-3］。在冰情严重的年份，春季流凌的撞击会对桥墩等水工建筑物造成不同程度的破坏，每年国家和沿海各省区都要投入巨大的人力物力和财力专门用于凌汛灾害的预防和处理，故监测河流断面的冰凌变化信息，掌握冰凌分布密度数据是冰凌灾害预防的关键［4-5］。

河流断面冰凌分布密度是在指定断面内，冰凌分布面积占河流单位面积的百分比，可以反映出河流断面上不同时间冰凌流量的大小，是预测冰凌灾害的重要指标［6］。目前获得冰凌分布密度常用的方法有人工望远镜目估法和卫星遥感监测法：（1）人工望远镜目估法［7］，通过人为观测，估算河流断面冰凌分布密度。这种方式操作简单，但对观测人员要求较高，需要有比较丰富的经验，既要保持眼睛的健康，也要保持清醒的头脑，但无法对冰凌分布密度进行准确地量化，缺乏客观性；（2）卫星遥感监测技术［8］，使用卫星拍摄的冰凌图像提取冰凌分布密度，由于卫星遥感方式的时间间隔大，空间分辨率较低等因素的影响，导致冰凌分布密度监测存在较大误差。国内外学者对海冰分布密度展开了大量的研究，提出了使用船载摄像机作为观测工具，对航行途径海域的海冰状况进行观测，Mu⁃ramoto 等［9］在1993年将摄像技术首次应用到海冰观测的实践中；Hall 等［10］使用船载摄像设备对海冰进行观测，基于阈值技术提取海冰分布密度，并使用SSM/I 的数据进行验证；Worby 等［11］使用多种卫星数据结合船上观测研究南极海冰变化等；卢鹏［12］利用传统图像形态学算法，构造出一种新的海冰边缘检测算法，提高拍摄图像中单个海冰识别效率；邓霄等［13］基于摄像机和透视变换图像校正算法设计了冰凌密度监测系统，改善了传回图像存在“近大远小”的问题；摄像机冰凌监测技术能够在一定程度上满足冰凌分布密度监测的要求，且图像校正算法能够改善图像存在的畸变，但由于图像校正算法自身的限制，对摄像机布置有一定的要求，且无法彻底消除图像畸变［14］。

随着无人机（unmanned aerial vehicles，UAV）和传感器技术的不断提高，小型无人机在低空遥感测量领域得到了广泛的应用［15］。对比人工望远镜目估法和卫星遥感监测法，可搭载多传感器的无人机平台具有效率高、成本低、操作灵活、更适合复杂河道环境等特点，不仅克服了人工获取影像数据的困难，还可以从根源上避免了岸边侧视拍摄导致的图像畸变，提高原始数据精度［16-18］，为监测大范围冰凌变化信息提供了新的技术手段。

本文通过无人机低空遥感技术对目标河段进行俯视拍摄，将视频逐帧提取图片，利用matlab 仿真软件编程，通过顶帽变换算法对图片进行亮度均衡化处理，基于OTSU算法（最大类间方差法）提取最佳阈值，对图像进行阈值分割，最终通过面积滤波降噪，计算河冰分布密度，绘制河冰分布密度图，实现对河冰分布密度连续监测。

2 研究区及数据预处理

2.1 研究区域概况黑龙江位于我国最北端，是中俄界河，流域处于寒温带地区，东亚大陆季风气候特点突出［19］，且黑龙江上游流经山区，河道狭窄曲折，雨量径流比较充沛，每年流凌期，冰凌洪水以叠加组合的形式促成倒开江向下游发展，导致漠河段易发生凌汛灾害（图1）［20］，黑龙江漠河段封冻期从11月份持续到次年4月，是典型的季节性封冻河流。

图1 冰凌观测地理位置示意

2.2 遥感影像获取及预处理本文试验数据于2021年4月25日通过无人机平台获取，采用大疆DJI四旋翼小型无人机御mavic 2 pro 进行定点拍摄，无人机体积为214 mm×91 mm×84 mm，可折叠机身，方便随身携带，搭载哈苏L1D-20c（CMOS），有效像素为2000 万，可录制4k HDR 视频，每次航拍的天气条件均为晴朗，风速低于5 级，飞行高度设置为400 m，起飞点与水面高程差为5 m，镜头垂直向下，悬停拍摄。

现场通过录制冰凌运动视频，内业处理时每隔29.97 帧截取图片影像，并通过图像处理工具对影像数据进行裁剪，去除边缘异常值。

为实现对河道冰凌分布图像数据的量化分析，首先，采用无人机拍摄不同高度下A1 纸，得到无人机遥感影像拍摄高度与单位像素点所代表的实际面积的拟合公式，如图2所示。

图2 无人机高度与单位像素点代表实际面积拟合曲线

3 研究方法

利用冰凌录像数据可以掌握被监测河道的实时冰凌分布密度与流速情况。在此基础上通过matlab编程对冰凌分布图像进行量化分析，提供新的流凌分布密度监测方式，为相关部门制定防凌方案提供理论依据和决策支持。

3.1 灰度图像转化通过摄像机拍摄视频能够实现河冰分布密度数据提取，主要是因为在太阳光照射条件下，河冰与河水对可见光反射的亮度和波长不同，使得它们在摄像机拍摄出来的影像中所显示的亮度有明显区别，即河冰要比河水更亮［21］。若将河冰分布图像转化为灰度图，这种不同也体现在灰度值上的差别。采用加权平均法将拍摄的彩色图像转化为灰度图像，计算公式如下所示：

式中：I（x，y）为灰度图像；R（x，y）、G（x，y）和B（x，y）分别为彩色图像的RGB 分量值。三个加权系数是根据人的亮度感知系统调节得到的标准化参数。

3.2 图像亮度均衡化处理采用灰度形态学中的组合算法，即顶帽变换算法对河冰图像进行亮度均衡化处理［22］。灰度形态学包括四种基本算法：灰度膨胀算法、灰度腐蚀算法、开运算和闭运算。灰度膨胀运算［23］是对原灰度图像进行“加长”或“变粗”；灰度腐蚀算法［24］是在灰度图像上“收缩”或“细化”；开运算［25］将前两种算法进行结合，对灰度图像进行先腐蚀运算，再膨胀运算，可以除去图像背景中比结构元素尺寸更小的亮度明亮细节；顶帽变换［26］则是在开运算对原始图像进行处理的基础上，用原始灰度图像再减去开运算结果的运算。这种方法能够进一步增强图像的对比度，适用于处理具有暗背景（河水）、亮物体（河冰）特征的图像，顶帽变换算法推导公式为：

腐蚀运算：

膨胀运算：

开运算：

顶帽变换计算：

式中：f（x，y）为图像M（x，y）为结构元素，AM是M（x，y）的定义域，结构元素M（x，y）对图像f（x，y）的灰度膨胀记为f+M，灰度腐蚀记为f-M，开运算记为f·M，顶帽变换记为fM。

3.3 OTSU算法阈值选取是阈值分割算法的核心［27］。通过确定合理的灰度阈值，可以实现图像分割，将河冰与河水区分开，计算河冰所占比例，从而得到河冰分布密度。

OTSU算法是借助最小二乘法原理在直方图技术上推导出来，具有简单、速度快等特点，是一种常用的阈值选取方法，适合于物体目标与背景灰度差明显的情况［28］。该算法以灰度分布均匀性的度量单位为方差，方差值越大，说明构成图像的两部分差别越大［29］。当部分目标（冰凌）错分为背景（水）或部分背景（水）错分为目标（冰凌）时，就会导致类间方差变小，因此使类间方差最大的分割阈值意味着计算精度变高［30］。

将原始冰凌图像f（x，y）进行灰度化处理，得到图像灰度范围为{0， 1，…，l-1}级。灰度值为i的像素数设为ni，推导出灰度图像总像素值为从而灰度值为i像素出现的概率为：pi=ni/N，其中选择阈值t，将灰度图像划分为两类，即D1：{0，1， 2，…，t}，D2：{t+1，t+2，t+3，…，l-1}，D1和D2出现的概率分别为：

D1和D2的灰度均值分别为：

图像的总体灰度均值为：

可求出D1和D2的类间方差为：

类间方差越大，说明D1和D2之间差距越大，即冰凌区域与水分割效果越好，所以确定最佳分割阈值，就是在求两类区域的类间方差最大值，即：

求出最佳阈值，即可实现图像分割，生成二值图像。

3.4 面积去噪在野外河道环境中，河流表面波浪反光往往是影响计算结果的重要因素。由于阳光反射产生的光点灰度值远高于图像分割阈值，导致反光点被错分为冰凌，冰凌分布密度增大，但反光点面积小且较为密集，只存在于图像的局部区域，不具备如高斯噪声［31］、椒盐噪声［32］等常见图像噪声的特征，无法通过均值滤波法［33］、自适应维纳滤波法［34］、中值滤波法［35］等去除。同时，无人机能够搭载的热成像传感器分辨率较低，难以满足白天流凌监测需求，因此通过对局部采用面积去噪法剔除反光噪声［36］，即：目标面积S小于X个像素点时，将会被剔除，实现图像分割。通过反复代入X值，比较图像的去噪效果，确定X取值区间。

4 结果与分析

4.1 顶帽算法对比分析在野外环境中，由于存在阳光照射角度，光照不均匀等因素会导致原始冰凌图像背景亮度不均匀，直接采用OTSU算法进行图像分割，会导致阈值分割失败，导致计算结果存在巨大误差，如图3所示。将采用顶帽变换算法前后得到的冰凌分布密度对比可以发现，两者之间相差最大为19.00%，最小为0.71%。采用顶帽变换算法后，计算获得的流凌分布密度均变大，但由于冰凌表面存在杂质，部分碎冰过薄颜色与江面接近，导致计算结果与人工选择冰凌得到的冰凌分布密度结果相比仍偏小。标准的光源和监测环境在野外环境下不可能达成，难以与高效数字化监测流程相结合，而采用顶帽变换算法可以补偿不均匀的背景亮度，使图像背景亮度均衡化，分割结果与原始冰凌图像（图4（a））情况一致。

图3 三种方法冰凌分布密度对比

从原始灰度图像三维图（图4（b））中可以得到图中边界区域参差不齐，并且从经过开运算得到的背景灰度三维图（图4（c））中可以得到河道中心区域灰度值高于河道两侧，河道左上角区域相对于周边区域灰度值略大，这是由于光照不均匀导致背景灰度不均匀。由于原始冰凌图像背景中河道中间区域和河道左上角区域灰度值较高，在通过OTSU算法选定阈值后，由于阈值选取过高，导致河道两侧冰凌和灰度值较小的薄冰、碎冰被错误的分割为河水，即河水面积扩大，冰凌分布密度偏小，与实际情况不符。对原始冰凌图像采用顶帽算法后，补偿灰度值不均匀的背景，从经过顶帽变换处理后的灰度三维图（图4（d））中可以看出，图像亮度均衡化，冰凌与河水成功分割。通过人为选取确定冰凌分布密度为46.23%，对比顶帽变换处理前后，冰凌分布密度误差由9.00%减小至0.97%，因此顶帽算法与OTSU算法相结合，能够解决图像亮度不均衡导致计算结果误差过大的问题。

图4 冰凌图像处理过程

4.2 OTSU 和迭代算法对比分析阈值分割主要由四种方法：人工选择法、直方图技术法、OTSU算法和迭代法。人工选择法［37］是通过肉眼识别判断阈值，并反复代入阈值，比较图像分割效果，缩小阈值选取区间，直至最终确定阈值，操作简单但效率低，且缺乏客观性，受到野外环境的影响，图像总体亮度也会随时间改变而改变，不适合批量处理所有冰凌图像；直方图技术法［38］适用于处理目标与背景灰度对比大的图像，即具有典型双峰直方图特征，通过人为选择两处波峰之间波谷的灰度值为分割阈值，如图5所示原始图像灰度值不存在两个明显的波峰，所以本研究不采用直方图技术法。迭代法［39］是通过选择初始阈值，基于阈值迭代逼近思想，直至得到满足给定的准则的最佳阈值；OTSU算法［40］按照图像的灰度特性将图像分为背景和目标两部分，以目标和背景的类间方差最大为阈值选择原则，自动确定最佳阈值。

图5 原始冰凌灰度直方图

本研究从2400 张原始冰凌图像中选取流凌初期、流凌中期和流凌末期3 个阶段的冰凌图像，根据冰凌的颜色、形态，人为将冰凌从图片中分割出来，转化为二值图像，将得到的冰凌分布密度作为该情况下的标准，再分别采用迭代法和OTSU算法进行阈值分割处理，冰凌分布密度处理效果如图6所示。

图6 冰凌分布密度计算效果对比

结果表明，对比OTSU算法和迭代算法的冰凌分布密度计算结果，流凌初期，OTSU算法和迭代算法与人工提取冰凌图像计算结果相差分别为1.39%和0.03%，图像分割效果较好；流凌中期，OTSU算法与人工提取冰凌图像误差最大为1019 m2；流凌末期，冰凌分布密度较低，由于碎冰、薄冰灰度值接近河水，对算法的计算精度要求更高，OTSU算法和迭代算法流凌分布密度计算结果相差可达3.6%，迭代算法相较于人工处理的冰凌分布密度，误差高达4360 m2，而OTSU算法相较于迭代算法更适合处理冰凌分布密度低的图像，OTSU算法与人工提取冰凌图像计算结果相差仅为0.5%。且对比OTSU算法和迭代算法的运行时间，两者相差约3.5 s，在大批量实时处理原始冰凌图像情况下，OT⁃SU 算法能够节约计算时间，且不同冰凌分布密度的图像分割精度高。

4.3 面积去噪分析由于波浪与阳光等因素地存在，部分时间段的图像会存在阳光反射的情况，在本研究中选取存在阳光反射的流凌图像，如图7（a）原始局部图像所示，左下角存在高密集，小面积的反光点，影响冰凌分布密度计算结果。反光点的灰度值（188 ～ 255）远大于图像分割阈值（129），导致冰凌分布密度偏大0.15%，通过以10 个像素的速度增加X值，发现当X值小于40 时，白色像素点数量迅速减小，当X值大于40 时，白色像素点减小速度变缓，因此确定该高度下噪声面积在40 个像素点以下，剔除反光噪声效果如图7（b）（c）所示，对比去噪前后的图像，该方法对于阳光反射此类噪声去除效果较好，且只针对图像的局部区域，不存在影响冰凌图像的形态特征。该方法实现简单，去噪效率高，可根据实际情况结合OTSU算法，进一步减小误差。

图7 采用面积去噪处理反光效果

4.4 冰凌分布密度计算结果无人机在一定高度悬停，定点拍摄视频，并从视频中每隔3 s 截取图片，同时在河岸边设置摄像机，每隔30 s 自动拍摄。首先将河道以外部分裁剪舍去，采用matlab 仿真软件编写程序，对2400 张无人机拍摄图像和2235 张摄像机拍摄图像进行阈值分割，得到阈值分割图像，计算结果作为每张图像对应的冰凌分布密度。冰凌分布密度随时间变化情况如图8所示，本次冰凌观测发生在2021年4月24日至26日之间，如图8所示，4月25日8 时前由于阳光照射角度的影响，且清晨可见光线不足，通过摄像机得到的冰凌图像分割效果差，导致冰凌分布密度误差大，计算结果偏小，与实际情况不符。4月25日8 时后，光线充足，均能实现冰凌图像分割，但摄像机拍摄图像存在“近大远小”的问题，且冰凌浮出水面一定高度，遮挡了冰凌间空隙，导致计算结果相较于实际情况偏大。摄像机拍摄冰凌图像选取的范围尺度较小，冰凌分布密度容易出现急剧波动。随着时间的推移，通过无人机采集的冰凌分布密度总体呈现下降趋势，期间监测到的冰凌分布密度最大为81.05%，且在2021年4月25日15 时45 分前，冰凌分布密度稳定在60%之上，4月25日15 时左右，上游发生冰塞，到4月25日15 时45 分后观测断面冰凌分布密度迅速下降，直至4月25日16 时35 分，冰凌分布密度下降至33.2%，随后回升至60%，符合现场实际情况，在4月26日8 时17 分以后，目标断面冰凌分布密度稳定在10%以下，无大面积冰凌存在，开江结束。

图8 冰凌分布密度曲线

5 结论

为了实时监测河流断面的冰凌变化信息，预防冰凌灾害，本文探讨了使用无人机低空遥感结合图像处理技术提取冰凌分布密度，利用matlab 仿真软件编程，最终获得了河冰分布密度随时间变化的动态曲线图，得到了以下结论：

（1）阈值分割法计算量小，性能稳定，在图像分割算法中应用最为广泛。阈值选取是阈值算法的核心，通过合理确定灰度阈值，从而实现图像分割，而最大类间方差法（OTSU）因其直观性和实现的简单性，类间方差越大，说明背景和目标两部分差异越大，适合大批量提取冰凌密度图像数据。

（2）冰凌分布密度计算结果受到光照不均匀的影响。顶帽变换算法通过对流凌观测图像进行亮度均衡化处理，与OTSU算法相结合，最终将与标准的误差缩小至1.5%以内，从而解决背景亮度不均匀条件下阈值的选择。

（3）面积图像分割法能够解决由于冰凌监测影像中波浪反光的影响，提高图像识别的准确性，而且该方法实现简单，去噪效率高。

（4）通过无人机低空遥感结合图像处理技术，本次观测流凌发生在2021年4月24日至26日，期间监测到的冰凌分布密度最大为81.05%，且在4月25日15 时45 分前，冰凌分布密度稳定在60%之上，在4月25日15 时左右，上游发生冰塞，从监测结果可以发现冰凌分布密度迅速下降，符合现场实际情况，在4月25日16 时35 分左右，上游冰塞解除，冰凌分布密度回升至60%，在4月26日8 时17 分以后，观测断面流凌分布密度稳定在10%以下，开江结束。