黄河内蒙古河段封冻期冰晶结构特征分析

2021-06-16李志军张邀丹

人民黄河 2021年6期

王娟，黄樾，李志军，邓宇，张邀丹

（1.郑州大学水利科学与工程学院，河南郑州450001；2.大连理工大学海岸和近海工程国家重点实验室，辽宁大连116042；3.黄河水利科学研究院，河南郑州450003）

每年冬季，我国北方的河流、湖泊、水库都会发生冰凌现象。冰凌会阻碍正常的航运、发电、供水，破坏大坝、护坡、桥梁等结构，还可能产生冰塞、冰坝等凌汛灾害，对人民的生命财产造成威胁。其中，针对黄河冰凌理论的研究，以河冰的生消演变和冰塞冰坝的形成机理为主［1-3］。河冰作为一种天然复合材料，不同地区、不同时间形成的河冰性质存在差异，这主要是由时空的差异造成了河冰细观结构的差别，而冰晶体结构很大程度上决定了河冰的物理力学性质［4-5］。研究表明，粒状冰的脆性强度与其晶体直径的平方根成反比［6］，多晶冰的极限抗压强度与粒径平方根的倒数成正比例关系［7］。因此，对晶体结构进行研究，有利于了解冰的物理性质，揭示冰的内在破坏机理［8］。

相关学者尝试过许多方法对冰晶体结构进行观测。19世纪中期，Seligman使用了在纸上制作铅笔拓片的方法［9］，使风化的晶界清晰地显示出来。近年来，基于偏振光的原理观测晶体结构的相关设备不断发展，如自动冰组构分析仪［10］、费氏台［11］等。2006年，Kipfstuhl等［12］使用电子显微镜和CCD摄像机对晶体结构进行了观察，通过显微结构映射法得到了晶体的高分辨率图片。其中，费氏台因操作方便、便于携带、试验条件要求低等优点，成为观测冰晶结构常采用的仪器。

冰晶观测技术一直在提高，但对于晶体图片的后续处理技术仍停留在人工操作的阶段。目前常用的方法是用Adobe Photo Shop CS软件配合手工勾选出晶体图片中晶体的边界，然后使用图像处理软件计算每个冰晶粒的尺寸［13］。该方法需要花费大量的时间，同时手工勾选边界存在误差。基于以上背景，本文基于MATLAB的数字图像分析能力，提出使用边界提取功能分割冰晶图像并计算晶粒尺寸的方法，分析黄河冰晶的结构分布特征。

1 冰晶图像的获取

2019年2月，在黄河内蒙古典型河段什四份弯道处（见图1（a））采集了4组冰胚。什四份弯道位于黄河内蒙古段纬度较高的位置，海拔990 m，每年冰期持续100 d左右。冰胚的取样点选取在黄河北岸距岸边60～150 m处。该处为黄河冬季主河道，冰层表面平整无积雪，冰下水深6～7 m，冰盖泥沙含量较低，排除冰凌堆积成冰。采集冰坯的具体位置见表1。采冰时首先利用GPS确定正北方向，按照约30 cm宽、30 cm长在冰面上切割冰块，使用电锯和人工板锯将冰胚切断取出。

图1 黄河冰样提取及冰晶图像获取

取样完成后，将冰样沿垂直方向分层，用锯骨机切成厚度相同的小块，并使用刨刀将冰块的一面修平，贴在温热的玻璃片上，使冰块与玻璃片充分黏结。随后，再使用刨刀将冰块削成厚度约1 mm的薄冰片，在费氏台偏光镜（见图1（b））下观测冰薄片，拍摄冰晶图片，并记录每张图片对应的实际尺寸。

表1 取样点坐标

2 冰晶图像处理

根据黄河冰晶结构观测结果，截取黄河冰晶图片的有效部分并转化为灰度图像，对图像进行滤波除噪和图像分割。基于边缘检测的方法提取冰晶体的边界，对部分边界特征不明显的区域进行校正补充。提取出完整边界之后，通过图像分析获取黄河冰晶的结构特征。

2.1 图像前处理

黄河冰的晶体结构原始图见图2（a），在原图中截取出最大内接矩形（见图2（b）），同时对照刻度尺记录矩形图片的实际尺寸，用于后期的单位换算。将图片转化为灰度图，提高图片处理的效率。使用偏光镜拍摄的冰晶图片具有较大的噪声，常见的空间域去噪方法有邻域平均法、选择平均法、中值滤波法、空间低通滤波法［14］。通过效果对比，选用中值滤波法对图像进行降噪处理。中值滤波的原理是在每一个像素点的周围取一个邻域，将该点的像素值取为邻域的中值，从而消除与周围像素相差较大的噪声。基于中值滤波的原理，分别使用边长为4像素、8像素、16像素、24像素、30像素的邻域对冰晶灰度图进行降噪处理。通过比较，最终选用边长为24像素的邻域进行滤波，见图2（c）。

图2 冰晶图像预处理

传统的图像分割方法有边缘检测法［15］、阈值分割法［16］、区域法［17］等。传统的阈值法适用于将图像分割为有限的两个或多个区域，不适用于晶粒较多的冰晶图像的分割。同时，冰晶图像噪声大，部分冰晶体的特征不明显，存在像素值接近、边界不明显的相邻晶粒。这些现象给使用区域法和聚类算法带来较大难度。基于上述情况，最终选择边缘检测方法对冰晶图像进行分割。常见的边缘检测算子有Robert算子、Sobel算子、Prewitt算子、Laplacian算子、Canny算子［18］等。前4种算子都属于局域窗口梯度算子，对噪声的抵抗性差，而Canny算子对噪声的抵抗性较强，故采用Canny算子对冰晶图像进行边界提取。从图3（a）可看出，大部分冰晶边界已经被提取出来，但仍存在少量边界的缺失和少量噪音。进一步处理后，边界清晰连续而且没有明显错位的现象（见图3（b））。

图3 冰晶图像分割

2.2 冰晶图像分析

采用MATLAB，通过编程将图3（b）中的连通区域进行编号并赋予不同的灰度值，每一个连通区域代表一个冰晶粒。通过灰度值的不同，计算每个连通区域像素点的数量以及其边界的像素点数量，根据记录的图像宽度，将像素面积和像素周长换算为实际面积和周长。由于冰晶粒形状不规则，因此用与晶粒面积相等的等效圆直径来描述冰晶粒的直径，面积和直径的计算公式如下:

式中:S为实际面积；D为等效直径；S n为区域像素点数量；k为每个像素点对应的实际长度。

由于直线边界和斜线边界每个像素点对应的实际长度不同，因此在计算实际周长L时需要将斜边界上的像素点额外乘以系数见式（3），其中L直和L斜分别为直线边界和斜线边界的像素点数量。

3 结果分析

3.1 冰晶粒总体分布规律

对共计50张黄河冰晶图片进行处理，获取7 000多个冰晶粒尺寸数据。等效直径小于3 mm的冰晶粒占12.99%，等效直径为3～6 mm的冰晶粒最多（见图4），占计算总量的31.94%。随着等效直径的增大，冰晶粒数量逐渐减少，等效直径大于21 mm的冰晶粒仅占计算总量的5.67%。

图4 冰晶粒尺寸总体分布

冰晶体类型可分为粒状冰和柱状冰两种，对两种不同的冰晶结构进行分析，见图5。粒状冰的晶粒分布与图4中总体分布的规律相似，其中:等效直径0～3 mm的晶粒数量小于3～6 mm的晶粒数量，等效直径大于3 mm的晶粒占比随着晶粒尺寸的增大逐渐减小。粒状冰的大尺寸晶粒较少，等效直径大于12 mm的晶粒仅占计算总量的6.57%。柱状冰尺寸偏大，小于3 mm的冰晶数量较少（仅占计算总量的6.57%），等效直径3～12 mm晶粒数量随晶粒等效直径的增大而减少，但大于12 mm的晶粒占比比粒状冰的多，晶粒数量下降的幅度相对平缓。这说明柱状冰的大尺寸冰晶粒数量更多，且尺寸分布较为均匀。

图5 粒状冰与柱状冰晶粒尺寸分布

3.2 冰晶粒水平分布规律

图6 展示了两组代表性冰胚的冰晶水平图像处理结果，左边冰样从上到下分别是1号冰胚8、16、24、32、40、48、56、63 cm深度的水平切片，右边冰样从上到下分别是2号冰胚5、15、25、35、45、55、68 cm深度的水平切片。从图6可看出，晶粒分割图与原图基本一致，边界连续完整，误差很小。图7（a）为计算得到的1号冰胚每层水平切片的等效直径分布，可以看出，8 cm和24 cm深度处的晶粒等效直径分布集中在0～3 mm和3～6 mm区间，大于6 mm的晶粒占比不到10%，原因是1号冰胚0～24 cm深度是粒状冰，晶粒较小。40～63 cm之间的冰晶体是柱状冰，随着深度的增加，大尺寸冰晶粒呈现逐渐增加的趋势。63 cm深度处的晶粒异常的原因是晶粒等效直径大，图片拍摄的晶粒数量少，统计规律不明显。图7（b）为计算得到的2号冰胚每层水平切片的等效直径分布，2号冰胚0～3 cm和20～35 cm深度是粒状冰，粒状冰层等效直径小于6 mm的晶粒占比达到了70%。随着深度的增加，柱状冰层晶粒占比增大。

图6 水平方向冰晶体图像识别

通过对两组冰晶体水平方向晶粒等效直径分布的对比可知，所有粒状冰层的晶粒占比均在3～6 mm区间达到最大，随后大幅度减小；所有柱状冰层的晶粒占比分布较为均匀，差异并不显著。

3.3 冰晶粒竖直分布规律

图7 不同深度水平方向冰晶粒尺寸分布

两组代表性的冰晶体垂直分布结果见图8。图8（a）为1号冰胚垂直切片的等效直径分布。深度0～24 cm是粒状冰层，晶粒等效直径维持在3～4 mm之间；深度24～63 cm是柱状冰层，晶粒等效直径随着深度的增加而增大，并在56 cm深度达到最大值11.5 mm；在深度63 cm处晶粒等效直径出现了小幅的减小，推测是冰花凝结的作用导致。图8（b）为2号冰胚垂直切片的等效直径分布，晶粒等效直径在5 cm深度达到极大值6.2 mm，随后又在25 cm深度减小至3.38 mm，再往下晶粒等效直径逐渐增大，并在68 cm深度达到最大值12.2 mm。冰晶粒等效直径的变化规律整体反映了冰晶结构粒状—柱状—粒状—柱状的交替变化趋势。