鸭子荡水库地形前处理及二维水流运动模拟

黄凌霄,李春光,郑兰香,杨 程

(1.宁夏大学土木与水利工程学院,宁夏 银川 750021; 2.宁夏大学信息工程学院,宁夏 银川 750021;3.北方民族大学土木工程学院,宁夏 银川 750021; 4.宁夏大学资源环境学院,宁夏 银川 750021)



鸭子荡水库地形前处理及二维水流运动模拟

黄凌霄1,2,李春光3,郑兰香4,杨 程3

(1.宁夏大学土木与水利工程学院,宁夏 银川 750021; 2.宁夏大学信息工程学院,宁夏 银川 750021;3.北方民族大学土木工程学院,宁夏 银川 750021; 4.宁夏大学资源环境学院,宁夏 银川 750021)

为了研究鸭子荡水库的二维水流问题,用图像处理方法结合谷歌地球(GE)软件对该水库进行地形的前处理。通过图像拼接、数学形态学和实测资料相结合的方法获取该水库的初始地形,对建立的平面二维RNGk-ε紊流数学模型中的水流模块给出合理的边界条件和初始条件,并采用有限体积法离散方程,通过阵面推进法划分网格,对水库流场和入水口、取水塔附近的断面垂线平均流速进行数值模拟。结果表明,利用图像处理方法结合GE软件能够得到精确、详细的水库复杂地形。在此基础上进行了水库流场数值模拟和关键位置的断面垂线平均流速数值模拟,模拟结果和实测结果较吻合。

图像拼接;SURF算法;线性渐变融合算法;数学形态学;平面二维水流模型;鸭子荡水库地形

宁夏灵武市的鸭子荡水库是宁东供水工程的核心,更是宁东能源化工基地最重要的基础设施之一,承担着为基地提供工业用水、生态用水和生活用水的重要任务。该水库自2003年开工建设,2004年利用抽水泵站抽取黄河水引入压力管道送进水库开始蓄水。水库水域面积300 hm2,最大库容2 400万m3,设计年取水量约1.3亿m3,2005年正式供水[1]。研究该水库在各种工况下运行的水沙运移规律,对指导水库管理方采用合理的运行方式,减轻水质污染、减缓泥沙淤积、延长水库的使用寿命具有重要意义。水库地形的描绘是研究水库水沙运移规律的前提条件,但鸭子荡水库水岸边界形状极不规则,若通过传统方法实测一些库岸边界的坐标值,对这些实测值进行三次样条插值后确定的库岸曲线失真较大,导致后续建模模拟结果与实际测量结果有较大的偏差。顾立忠等[2]通过AutoCAD电子图获取河道高程信息,提取高程点数据,结合MIKE21软件生成河道地形。关见朝等[3]利用谷歌地球(Google Earth,GE)软件提取河道高程数据,生成初始地形。赵焱[4]使用ArcGIS软件对河道测量高程数据进行插值、读取、转换等操作后绘制出河道轮廓。目前大多数研究者都是通过各种软件在获取高程数据后进行相关操作,生成河道的初始地形,但详细真实的高程数据获取较困难,少量高程数据需要进行插值运算。本文通过GE软件获取水库的经纬度信息(经纬度数据比高程数据更真实有效,也比利用高程数据所绘制出的地形更加详细),并在获取水库的初始地形之后,建立平面二维紊流数学模型,对水库的水流运动进行数值模拟,找出水流运动的相关规律,旨在为后续的水质分析、泥沙运移及泥沙淤积研究奠定基础。

1 水库地形前处理

为了得到精确的水库库岸曲线,利用GE软件获取水库的局部图像,通过图像拼接技术中的加速鲁棒特性(speeded-up Robust features,SURF)算法[5]提取水库局部图像的特征点,使用k-d tree算法[6]对特征点进行匹配,并使用线性渐变融合算法将多张局部图像拼接成一张完整的水库宽视野图像。在此基础上,利用数学形态学方法[7]检测水库库岸曲线,将检测出的结果图像导入AutoCAD软件中,绘制水库库岸曲线并保存为相关文件,利用程序把该文件加载到GE软件中,可以提取水库库岸曲线的经纬度信息,结合实测资料可以精确地得到鸭子荡水库的地形。

1.1 图像拼接

在研究过程中,通过GE软件获取的全局水库图像由于拍摄高度过高,不能得到精确的水库库岸曲线，只能通过图像拼接技术将多张具有重叠部分的局部水库图像拼合成一张宽视野的完整水库图像。由于图像拼接技术中的SURF算法具有良好的尺度、光照、旋转不变等特性而被广泛应用，因此采用SURF算法进行水库局部图像特征点提取。

1.1.1 SURF算法提取特征点

SURF算法的核心是积分图像和Hessian矩阵。积分图像是某一矩形范围内的像素之和,一旦一幅图像的积分图像提前计算好,那么该图像中任意矩形范围内的像素之和都可以通过3个加(减)法运算来完成,这说明计算时间与矩形范围的大小无关,使SURF算法能够极大地缩短计算时间。Hessian矩阵是SURF算法对特征点提取的基础,对于图像I(X)上的一点X=(x,y),Hessian矩阵H(X,δ)在X处尺度为δ的定义如下:

(1)

(2)

Hessian矩阵的行列式为

(3)

SURF算法使用方框滤波器替换二阶高斯滤波器,这样可以利用积分图像加速卷积,提高计算速度。方框滤波器同图像I在3个方向的卷积分别为Dxx、Dyy和Dxy,式(3)可近似为

ΔH=DxxDyy-(0.9Dxy)2

(4)

通过计算图像中每个像素点的Hessian矩阵行列式的值,根据值的正负可以判断出该点是否为极值点。构建多尺度的图像金字塔,就可以在不同尺度上寻找特征点。对其进行非极大值抑制,可以进一步得到精确的特征点位置及所在尺度值。

为了实现特征点描述子具有旋转不变性,需要确定特征点的主方向。在以特征点为圆心、6δ为半径的圆内,计算x、y方向上的Haar小波(Haar小波的边长为4δ)响应值。以特征点为中心对这些响应值进行高斯加权,靠近特征点的权重大,远离特征点的权重小,在π/3的扇形区域内滑动窗口,计算x、y方向的响应总和,构成一个新的向量,转动π/3的扇形范围遍历整个圆,选择其中最长的向量方向作为特征点的主方向。

1.1.2k-dtree算法进行特征匹配

特征匹配就是在高维向量空间中寻找最相似的特征向量。为提高匹配点的搜索效率,采用k-d tree算法来获取匹配点。该算法对于某一特征向量,首先计算待匹配图像中所有特征向量与该向量的距离,然后寻找可能的最近邻点和次近邻点,若最近邻距离和次近邻距离的比值小于预先设定的阈值,则认为最近邻点即为该特征点的匹配点,两点组成匹配点对;若找不到或不满足条件,则认为该特征点没有匹配点。

1.1.3 线性渐变融合算法实现图像拼接

直接将水库局部图像拼接到一起,会有明显的拼接痕迹,因此需要采用合适的融合策略来解决重叠区域的拼接问题,本文采用线性渐变融合算法来实现图像的拼接，计算公式为

(5)

式中:I1(x,y)、I2(x,y)分别为图像I1和图像I2在其对应点的像素值;P1、P2为权重系数,取值范围是[0,1];xmax、xmin分别为图像I1和图像I2在重叠区域的x轴最大值和最小值。

通过以上算法可以解决重叠区域的拼接问题。

通过GE软件获取了3张水库局部图像(图1),分别为水库取水塔、水库中部和水库入水口。通过SURF算法分别对这3张水库局部图像进行特征点提取,利用k-d tree算法对图1(a)(b)以及图1(b)(c)进行特征匹配,如图2、图3所示。在此基础上,使用线性渐变融合算法对3张水库局部图像进行无缝拼接,从而得到完整水库图像,如图4所示。

图1 鸭子荡水库局部图像

图2 水库取水塔与水库中部特征匹配图像

图3 水库中部与水库入水口特征匹配图像

图4 拼接结果图像

1.2 水库库岸曲线检测

通过图像拼接技术得到完整的水库图像后,可以人工绘制出水库的库岸曲线。由于人眼视觉差别，水库库岸曲线提取有较大的偏差。因此,需要使用图像边缘检测技术对水库库岸曲线进行精确提取。近年来,图像边缘检测技术中数学形态学方法的应用越来越广泛，该方法不但能有效地去除图像中的噪声,还能较多地保留图像的边缘信息。

设原始图像A和结构元素B是整数空间Z中的集合。B对A的膨胀记为A⨁B,定义为

(6)

经过膨胀后的图像将比原始图像像素更多,该运算可以填充图像边缘的小凹陷。

腐蚀和膨胀是对偶操作。B对A的腐蚀记为AΘB,定义为

(7)

式中(B)z表示对B进行平移z(z1,z2)。经过腐蚀后的图像将比原始图像像素更少,该运算可以消除图像中小且无意义的部分。

利用数学形态学方法先将图4转换为灰度图像,选取半径为5的圆盘为结构元素,这样能保证检测到较清晰的水库库岸曲线,再对膨胀运算后的灰度图像与腐蚀运算后的灰度图像进行减法运算,结果如图5所示。将图5导入AutoCAD软件中进行描绘并保存,将保存结果转化到GE软件中拟合,结果如图6所示。

图5 水库库岸曲线检测图像

图6 水库库岸曲线拟合图像

1.3 水库地形及断面分布

提取拟合后的水库库岸曲线经纬度信息,结合现场实测的水库库岸散点经纬度信息,生成Excel数据表,导入MATLAB软件中生成鸭子荡水库的初始地形。2014年7月和11月利用船载声学多普勒流速剖面仪、RTK、GPS等仪器实测鸭子荡水库的断面资料,自入水口至取水塔方向共布设17个断面,其中7个断面由于水库复杂地形的限制被分割为左右两部分。通过提取仪器中断面散点的经纬度信息,生成Excel数据表,导入MATLAB软件中生成鸭子荡水库的断面分布,其中入水口位于断面3附近,取水塔位于断面13与断面14之间,结果如表1所示。

2 数学模型

2.1 基本方程

沿水深平均的平面二维RNGk-ε紊流数学模型[8-12]包含2个模块:水流模块和泥沙模块。本文主要研究鸭子荡水库的水流运动,而水流模块包括水流连续性方程、x方向动量方程、y方向动量方程、紊动动能k方程和紊动动能耗散率ε方程。这些方程可以统一写成:

(8)

表1 鸭子荡水库断面经纬度

式中:h为水位;Φ为通用变量;u为x方向上的流速;v为y方向上的流速;ΓΦ为扩散系数;SΦ(x,y)为源项。其中,Φ、ΓΦ和SΦ(x,y)在不同的水流模块方程中代表的变量不同,具体变量含义及其他相关公式可以参见文献[8]。

2.2 方程的离散及定解条件

使用有限体积法[13]对式(8)进行离散,其中瞬时项采用向前差分格式，对流项采用乘方格式，扩散项采用中心差分格式，源项采用局部线性化处理,具体处理过程见文献[12]。

进口边界给定流量,紊动动能k=0.01Uin2,Uin为进口边界的平均流速,可以通过进口流量与断面垂线平均流速分布公式得到,紊动动能耗散率ε=0.09k1.5/(0.05H),H为进口水深,可以通过实测值给定;出口边界给定水位和流量;近岸固壁边界按照无滑移处理。

2.3 计算区域及网格划分

通过上述数学模型对鸭子荡水库进行二维水流编程模拟计算,水库的初始水位为1 247.0 m,计算区域、入水口和取水塔等位置如图7所示,自入水口至取水塔方向共布设17个断面,其中断面5～9与断面13、14被库岸曲线分为左右两部分。采用阵面推进法[14-15]在计算区域内生成三角形网格,共有7 267个节点,11 919个单元,如图7所示。

图7 鸭子荡水库网格划分示意图

3 水流实测结果与模拟结果对比分析

根据鸭子荡水库目前实际运行方式,设计工况一如下:当水库水位降至1 247.0 m时,采取泵站抽取黄河水注入水库的方式持续进水,日进水量为60万m3,日取水量为40万m3,直至水库水位升至1 249.5 m,此时流场如图8(a)所示。随着宁东能源化工基地工业用水与生活用水的日益增加,设计工况二如下:当水库水位降至1 247.0 m时进行补水,日进水量为80万m3,日取水量为60万m3,直至水库水位升至1 249.5 m,此时流场如图8(b)所示。

图8 水流流场

分析图8可知:随着进水量和取水量的增加,2种工况运行下的水流流场都与实际情况相吻合;随着进水量的增加,进水口附近的流速增加,进水口进入库区的水流对水库上部的流场影响增大,尤其水库上部的主流区域和岔道区域流场变化明显;由于水库的中部长且浅,对进入库区的水流起到了限制和减缓作用,因此进水量增加对水库中部的流场影响较小,对取水塔的流场影响也较小;取水塔上下两侧形成2个回流区,并且下侧的回流区比上侧的大;随着取水量的增大,取水塔附近的流速增加,取水塔上下两侧2个回流区流场变化明显,水库下部的岔道区域流场变化也较明显。

2014年7月和11月利用声学多普勒流速剖面仪、RTK等仪器对水库的17个断面进行了实测,选择2014年11月实测期间的数据进行数值模拟,模拟计算的初始水位为1 247.0 m,通过上述数学模型的模拟计算,经过24 d后水位升至1 249.5 m。选择进水口附近的第3断面、第4断面的右部和取水塔附近的第14断面、第15断面进行实测值与模拟值比较分析,如图9所示。

图9 断面垂线平均流速对比

分析图9可知:模拟值与实测值总体上比较吻合,说明水库地形的前处理是合适的,数学模型的建立是有效的;从断面3、断面4右部及其他断面的实测垂线平均流速可以看出,由于断面右岸靠近入水口,断面右岸附近的流速大于断面左岸附近的流速,随着断面与入水口之间距离的增加,断面右岸的流速逐渐较小,断面左岸的流速变化很小;从断面14、断面15及其他断面的实测垂线平均流速可以看出,由于断面左岸靠近取水塔,断面左岸附近的流速大于断面右岸附近的流速,随着断面与取水塔之间距离的增加,断面左岸的流速减小得很快,断面右岸的流速变化很小。

4 结 语

以位于宁夏灵武市的鸭子荡水库为研究对象,通过图像拼接方法得到完整的鸭子荡水库图像,利用数学形态学方法和AutoCAD软件提取鸭子荡水库的库岸曲线,将提取后的经纬度信息与实测资料相结合,构建了鸭子荡水库的初始地形及断面分布。在此基础上,建立了二维水流数学模型,给出合理的边界条件和初始条件,使用有限体积法对相关方程进行离散,对计算区域进行三角形网格划分,实现了鸭子荡水库水流实测结果与模拟结果的对比分析。结果表明:图像拼接、数学形态学等方法相结合实现了水库初始地形的前处理,可以显著提高模拟过程中的地形处理效率和精度;模型计算的二维流场分布、断面垂线平均流速与实测资料比较吻合,说明所建立的数学模型是合理的,为下一步研究鸭子荡水库的泥沙运移及泥沙淤积情况奠定了基础。

[ 1 ] 齐敦哲,金满洋,赵平均,等.宁夏宁东鸭子荡水库渔业资源调查报告[J].水生态学杂志,2009,2(5):107-110.(QI Dunzhe,JIN Manyang,ZHAO Pingjun,et al.Survey on fishery resources of Yazidang Reservoir in Ningdong Ningxia[J].Journal of Hydroecology,2009,2(5):107-110.(in Chinese))

[ 2 ] 顾立忠,关琨,郑国栋,等.浅谈MIKE21模型地形前处理技巧[J].广东水利水电,2008,11(8):122-124.(GU Lizhong,GUAN Kun,ZHENG Guodong,et al.Talking About the MIKE21 model terrain pretreatment technique[J].Guangdong Water Resources and Hydropower,2008,11(8):122-124.(in Chinese))

[ 3 ] 关见朝,方春明.基于Google Earth的河流模拟地形前处理新方法[J].水利水电技术,2011,42(12):21-24.(GUAN Jianzhao,FANG Chunming.A Google Earth based new approach to pretreatment of terrain for river simulation[J].Water Resources and Hydropower Engineering,2011,42(12):21-24.(in Chinese))

[ 4 ] 赵焱.基于ArcGIS 10.0的河道地形高程数据处理方法[J].电子制作,2015(5):72-73.(ZHAO Yan.River terrain elevation data processing method based on ArcGIS 10.0[J].Practical Electronics,2015(5):72-73.(in Chinese))

[ 5 ] BAY H,ESS A,TUYTELAARS T,et al.SURF:speeded-up robust features[J].Computer Vision and Image Understanding,2008,110(3):346-359.

[ 6 ] BROWN R A.Building a balanced k-d tree in O (knlogn) time[J].Journal of Computer Graphics Techniques,2015,4(1):50-68.

[ 7 ] XU Chen,LIU Hui,CAO Wenming,et al.Multispectral image edge detection via Clifford gradient[J].Science China(Information Sciences),2012,55(2):260-269.

[ 8 ] 李春光,景何仿,吕岁菊,等.水沙水质类水利工程问题数值模拟理论与应用[M].北京:科学出版社,2015.

[ 9 ] 欧剑,马进荣,张行南,等.大通至长江口整体水动力模型[J].河海大学学报(自然科学版),2009,37(3):258-262.(OU Jian,MA Jinrong,ZHANG Xingnan,et al.Overall hydrodynamic model of Datong-Yangtze River Estuary[J].Journal of Hohai University (Natural Sciences),2009,37(3):258-262.(in Chinese))

[10] 陈学群,李福林,张瑞青,等.黄河河口平原多闸坝河道水流数学模型[J].水资源保护,2012,28(1):38-41.(CHEN Xuequn,LI Fulin,ZHANG Ruiqing,et al.Mathematical model of flow in rivers with multiple sluices and dams in plain erea in Yellow River Estuary[J].Water Resources Protection,2012,28(1):38-41.(in Chinese))

[11] 辛小康,肖洋,朱晓丹,等.基于DGA的BP神经网络及其在一维河网模拟中的应用[J].水利水电科技进展,2009,29(3):9-13.(XIN Xiaokang,XIAO Yang,ZHU Xiaodan,et al.BP neural network based on DGA and its application in one-dimensional river network simulation[J].Advances in Science and Technology of Water Resources,2009,29(3):9-13.(in Chinese))

[12] 景何仿,李春光.黄河上游连续弯道水流运动及泥沙运移数值模拟研究[M].郑州:黄河水利出版社,2012.

[13] 熊伟,翟剑峰,朱志夏,等.基于有限体积法的三维波生近岸流数值模型[J].水动力学研究与进展:A辑,2016,31(1):43-49.(XIONG Wei,ZHAI Jianfeng,ZHU Zhixia,et al.Three-dimensional nearshore currnet model based on finite volume method[J].Chinese Journal of Hydrodynamics,2016,31(1):43-49.(in Chinese))

[14] LÖHNER R.Rencent Advances in parallel advancing front grid generation[J].Archives of Computational Methods in Engineering,2014,21(2):1-14.

[15] MAZZOLARI A,TRIGO-TEIXEIRA A.Improved advancing front mesh algorithm with pseudoislands as internal fronts[J].Journal of Waterway Port Coastal & Ocean Engineering,2014,140(4):86-95.

Terrain pretreatment and two-dimensional flow simulation in Yazidang Reservoir

HUANG Lingxiao1,2, LI Chunguang3, ZHENG Lanxiang4, YANG Cheng3

(1.School of Civil and Hydraulic Engineering, Ningxia University, Yinchuan 750021, China; 2.School of Information Engineering, Ningxia University, Yinchuan 750021, China; 3.School of Civil Engineering, Beifang University of Nationalities, Yinchuan 750021, China; 4.School of Resources and Environment, Ningxia University, Yinchuan 750021, China)

In order to study the two-dimensional (2D) flow in the Yazidang Reservoir, terrain pretreatment was carried out using the image processing method and Google Earth (GE) software. The initial terrain of the reservoir was obtained through image mosaicking, mathematical morphology, and measured data. A 2D RNGk-εturbulence model was established, reasonable boundary conditions and initial conditions of the flow module of the established model were set, the equations were discretized using the finite volume method, and a mesh was divided using the advancing-front method. The flow field in the reservoir and the depth-averaged velocity at sections near the water inlet and the water power were simulated. The results show that accurate and detailed reservoir terrain can be obtained using the image processing method and the GE software. The numerical simulation results of the flow field and the depth-averaged velocity at critical positions show that the simulated and measured values are in strong agreement.

image matching; SURF algorithm; linear gradient fusion algorithm; mathematical morphology; 2D flow model; terrain of Yazidang Reservoir

国家自然科学基金地区科学基金(61362029);宁夏高等学校科学研究项目(NGY2015139);宁夏大学科学研究基金(ZR16018)

黄凌霄(1977—),男,副教授,博士研究生,主要从事计算水力学及数字图像处理研究。E-mail:hlingxiao@126.com

李春光(1964—),男,教授,博士,主要从事计算水力学研究。E-mail: cglizd@hot.mail.com

10.3880/j.issn.1006-7647.2017.04.009

TV147；TP751

A

1006-7647(2017)04-0047-06

2016-08-17 编辑：高建群)