鴨子蕩水庫地形前處理及二維水流運動模擬

黃凌霄,李春光,鄭蘭香,楊 程

(1.寧夏大學土木與水利工程學院,寧夏 銀川 750021; 2.寧夏大學信息工程學院,寧夏 銀川 750021;3.北方民族大學土木工程學院,寧夏 銀川 750021; 4.寧夏大學資源環(huán)境學院,寧夏 銀川 750021)

?

鴨子蕩水庫地形前處理及二維水流運動模擬

黃凌霄1,2,李春光3,鄭蘭香4,楊 程3

(1.寧夏大學土木與水利工程學院,寧夏 銀川 750021; 2.寧夏大學信息工程學院,寧夏 銀川 750021;3.北方民族大學土木工程學院,寧夏 銀川 750021; 4.寧夏大學資源環(huán)境學院,寧夏 銀川 750021)

為了研究鴨子蕩水庫的二維水流問題,用圖像處理方法結(jié)合谷歌地球(GE)軟件對該水庫進行地形的前處理。通過圖像拼接、數(shù)學形態(tài)學和實測資料相結(jié)合的方法獲取該水庫的初始地形,對建立的平面二維RNGk-ε紊流數(shù)學模型中的水流模塊給出合理的邊界條件和初始條件,并采用有限體積法離散方程,通過陣面推進法劃分網(wǎng)格,對水庫流場和入水口、取水塔附近的斷面垂線平均流速進行數(shù)值模擬。結(jié)果表明,利用圖像處理方法結(jié)合GE軟件能夠得到精確、詳細的水庫復雜地形。在此基礎(chǔ)上進行了水庫流場數(shù)值模擬和關(guān)鍵位置的斷面垂線平均流速數(shù)值模擬,模擬結(jié)果和實測結(jié)果較吻合。

圖像拼接;SURF算法;線性漸變?nèi)诤纤惴?數(shù)學形態(tài)學;平面二維水流模型;鴨子蕩水庫地形

寧夏靈武市的鴨子蕩水庫是寧東供水工程的核心,更是寧東能源化工基地最重要的基礎(chǔ)設施之一,承擔著為基地提供工業(yè)用水、生態(tài)用水和生活用水的重要任務。該水庫自2003年開工建設,2004年利用抽水泵站抽取黃河水引入壓力管道送進水庫開始蓄水。水庫水域面積300 hm2,最大庫容2 400萬m3,設計年取水量約1.3億m3,2005年正式供水[1]。研究該水庫在各種工況下運行的水沙運移規(guī)律,對指導水庫管理方采用合理的運行方式,減輕水質(zhì)污染、減緩泥沙淤積、延長水庫的使用壽命具有重要意義。水庫地形的描繪是研究水庫水沙運移規(guī)律的前提條件,但鴨子蕩水庫水岸邊界形狀極不規(guī)則,若通過傳統(tǒng)方法實測一些庫岸邊界的坐標值,對這些實測值進行三次樣條插值后確定的庫岸曲線失真較大,導致后續(xù)建模模擬結(jié)果與實際測量結(jié)果有較大的偏差。顧立忠等[2]通過AutoCAD電子圖獲取河道高程信息,提取高程點數(shù)據(jù),結(jié)合MIKE21軟件生成河道地形。關(guān)見朝等[3]利用谷歌地球(Google Earth,GE)軟件提取河道高程數(shù)據(jù),生成初始地形。趙焱[4]使用ArcGIS軟件對河道測量高程數(shù)據(jù)進行插值、讀取、轉(zhuǎn)換等操作后繪制出河道輪廓。目前大多數(shù)研究者都是通過各種軟件在獲取高程數(shù)據(jù)后進行相關(guān)操作,生成河道的初始地形,但詳細真實的高程數(shù)據(jù)獲取較困難,少量高程數(shù)據(jù)需要進行插值運算。本文通過GE軟件獲取水庫的經(jīng)緯度信息(經(jīng)緯度數(shù)據(jù)比高程數(shù)據(jù)更真實有效,也比利用高程數(shù)據(jù)所繪制出的地形更加詳細),并在獲取水庫的初始地形之后,建立平面二維紊流數(shù)學模型,對水庫的水流運動進行數(shù)值模擬,找出水流運動的相關(guān)規(guī)律,旨在為后續(xù)的水質(zhì)分析、泥沙運移及泥沙淤積研究奠定基礎(chǔ)。

1 水庫地形前處理

為了得到精確的水庫庫岸曲線,利用GE軟件獲取水庫的局部圖像,通過圖像拼接技術(shù)中的加速魯棒特性(speeded-up Robust features,SURF)算法[5]提取水庫局部圖像的特征點,使用k-d tree算法[6]對特征點進行匹配,并使用線性漸變?nèi)诤纤惴▽⒍鄰埦植繄D像拼接成一張完整的水庫寬視野圖像。在此基礎(chǔ)上,利用數(shù)學形態(tài)學方法[7]檢測水庫庫岸曲線,將檢測出的結(jié)果圖像導入AutoCAD軟件中,繪制水庫庫岸曲線并保存為相關(guān)文件,利用程序把該文件加載到GE軟件中,可以提取水庫庫岸曲線的經(jīng)緯度信息,結(jié)合實測資料可以精確地得到鴨子蕩水庫的地形。

1.1 圖像拼接

在研究過程中,通過GE軟件獲取的全局水庫圖像由于拍攝高度過高,不能得到精確的水庫庫岸曲線，只能通過圖像拼接技術(shù)將多張具有重疊部分的局部水庫圖像拼合成一張寬視野的完整水庫圖像。由于圖像拼接技術(shù)中的SURF算法具有良好的尺度、光照、旋轉(zhuǎn)不變等特性而被廣泛應用，因此采用SURF算法進行水庫局部圖像特征點提取。

1.1.1 SURF算法提取特征點

SURF算法的核心是積分圖像和Hessian矩陣。積分圖像是某一矩形范圍內(nèi)的像素之和,一旦一幅圖像的積分圖像提前計算好,那么該圖像中任意矩形范圍內(nèi)的像素之和都可以通過3個加(減)法運算來完成,這說明計算時間與矩形范圍的大小無關(guān),使SURF算法能夠極大地縮短計算時間。Hessian矩陣是SURF算法對特征點提取的基礎(chǔ),對于圖像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矩陣行列式的值,根據(jù)值的正負可以判斷出該點是否為極值點。構(gòu)建多尺度的圖像金字塔,就可以在不同尺度上尋找特征點。對其進行非極大值抑制,可以進一步得到精確的特征點位置及所在尺度值。

為了實現(xiàn)特征點描述子具有旋轉(zhuǎn)不變性,需要確定特征點的主方向。在以特征點為圓心、6δ為半徑的圓內(nèi),計算x、y方向上的Haar小波(Haar小波的邊長為4δ)響應值。以特征點為中心對這些響應值進行高斯加權(quán),靠近特征點的權(quán)重大,遠離特征點的權(quán)重小,在π/3的扇形區(qū)域內(nèi)滑動窗口,計算x、y方向的響應總和,構(gòu)成一個新的向量,轉(zhuǎn)動π/3的扇形范圍遍歷整個圓,選擇其中最長的向量方向作為特征點的主方向。

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

特征匹配就是在高維向量空間中尋找最相似的特征向量。為提高匹配點的搜索效率,采用k-d tree算法來獲取匹配點。該算法對于某一特征向量,首先計算待匹配圖像中所有特征向量與該向量的距離,然后尋找可能的最近鄰點和次近鄰點,若最近鄰距離和次近鄰距離的比值小于預先設定的閾值,則認為最近鄰點即為該特征點的匹配點,兩點組成匹配點對;若找不到或不滿足條件,則認為該特征點沒有匹配點。

1.1.3 線性漸變?nèi)诤纤惴▽崿F(xiàn)圖像拼接

直接將水庫局部圖像拼接到一起,會有明顯的拼接痕跡,因此需要采用合適的融合策略來解決重疊區(qū)域的拼接問題,本文采用線性漸變?nèi)诤纤惴▉韺崿F(xiàn)圖像的拼接，計算公式為

(5)

式中:I1(x,y)、I2(x,y)分別為圖像I1和圖像I2在其對應點的像素值;P1、P2為權(quán)重系數(shù),取值范圍是[0,1];xmax、xmin分別為圖像I1和圖像I2在重疊區(qū)域的x軸最大值和最小值。

通過以上算法可以解決重疊區(qū)域的拼接問題。

通過GE軟件獲取了3張水庫局部圖像(圖1),分別為水庫取水塔、水庫中部和水庫入水口。通過SURF算法分別對這3張水庫局部圖像進行特征點提取,利用k-d tree算法對圖1(a)(b)以及圖1(b)(c)進行特征匹配,如圖2、圖3所示。在此基礎(chǔ)上,使用線性漸變?nèi)诤纤惴▽?張水庫局部圖像進行無縫拼接,從而得到完整水庫圖像,如圖4所示。

圖1 鴨子蕩水庫局部圖像

圖2 水庫取水塔與水庫中部特征匹配圖像

圖3 水庫中部與水庫入水口特征匹配圖像

圖4 拼接結(jié)果圖像

1.2 水庫庫岸曲線檢測

通過圖像拼接技術(shù)得到完整的水庫圖像后,可以人工繪制出水庫的庫岸曲線。由于人眼視覺差別，水庫庫岸曲線提取有較大的偏差。因此,需要使用圖像邊緣檢測技術(shù)對水庫庫岸曲線進行精確提取。近年來,圖像邊緣檢測技術(shù)中數(shù)學形態(tài)學方法的應用越來越廣泛，該方法不但能有效地去除圖像中的噪聲,還能較多地保留圖像的邊緣信息。

設原始圖像A和結(jié)構(gòu)元素B是整數(shù)空間Z中的集合。B對A的膨脹記為A?B,定義為

(6)

經(jīng)過膨脹后的圖像將比原始圖像像素更多,該運算可以填充圖像邊緣的小凹陷。

腐蝕和膨脹是對偶操作。B對A的腐蝕記為AΘB,定義為

(7)

式中(B)z表示對B進行平移z(z1,z2)。經(jīng)過腐蝕后的圖像將比原始圖像像素更少,該運算可以消除圖像中小且無意義的部分。

利用數(shù)學形態(tài)學方法先將圖4轉(zhuǎn)換為灰度圖像,選取半徑為5的圓盤為結(jié)構(gòu)元素,這樣能保證檢測到較清晰的水庫庫岸曲線,再對膨脹運算后的灰度圖像與腐蝕運算后的灰度圖像進行減法運算,結(jié)果如圖5所示。將圖5導入AutoCAD軟件中進行描繪并保存,將保存結(jié)果轉(zhuǎn)化到GE軟件中擬合,結(jié)果如圖6所示。

圖5 水庫庫岸曲線檢測圖像

圖6 水庫庫岸曲線擬合圖像

1.3 水庫地形及斷面分布

提取擬合后的水庫庫岸曲線經(jīng)緯度信息,結(jié)合現(xiàn)場實測的水庫庫岸散點經(jīng)緯度信息,生成Excel數(shù)據(jù)表,導入MATLAB軟件中生成鴨子蕩水庫的初始地形。2014年7月和11月利用船載聲學多普勒流速剖面儀、RTK、GPS等儀器實測鴨子蕩水庫的斷面資料,自入水口至取水塔方向共布設17個斷面,其中7個斷面由于水庫復雜地形的限制被分割為左右兩部分。通過提取儀器中斷面散點的經(jīng)緯度信息,生成Excel數(shù)據(jù)表,導入MATLAB軟件中生成鴨子蕩水庫的斷面分布,其中入水口位于斷面3附近,取水塔位于斷面13與斷面14之間,結(jié)果如表1所示。

2 數(shù)學模型

2.1 基本方程

沿水深平均的平面二維RNGk-ε紊流數(shù)學模型[8-12]包含2個模塊:水流模塊和泥沙模塊。本文主要研究鴨子蕩水庫的水流運動,而水流模塊包括水流連續(xù)性方程、x方向動量方程、y方向動量方程、紊動動能k方程和紊動動能耗散率ε方程。這些方程可以統(tǒng)一寫成:

(8)

表1 鴨子蕩水庫斷面經(jīng)緯度

式中:h為水位;Φ為通用變量;u為x方向上的流速;v為y方向上的流速;ΓΦ為擴散系數(shù);SΦ(x,y)為源項。其中,Φ、ΓΦ和SΦ(x,y)在不同的水流模塊方程中代表的變量不同,具體變量含義及其他相關(guān)公式可以參見文獻[8]。

2.2 方程的離散及定解條件

使用有限體積法[13]對式(8)進行離散,其中瞬時項采用向前差分格式，對流項采用乘方格式，擴散項采用中心差分格式，源項采用局部線性化處理,具體處理過程見文獻[12]。

進口邊界給定流量,紊動動能k=0.01Uin2,Uin為進口邊界的平均流速,可以通過進口流量與斷面垂線平均流速分布公式得到,紊動動能耗散率ε=0.09k1.5/(0.05H),H為進口水深,可以通過實測值給定;出口邊界給定水位和流量;近岸固壁邊界按照無滑移處理。

2.3 計算區(qū)域及網(wǎng)格劃分

通過上述數(shù)學模型對鴨子蕩水庫進行二維水流編程模擬計算,水庫的初始水位為1 247.0 m,計算區(qū)域、入水口和取水塔等位置如圖7所示,自入水口至取水塔方向共布設17個斷面,其中斷面5～9與斷面13、14被庫岸曲線分為左右兩部分。采用陣面推進法[14-15]在計算區(qū)域內(nèi)生成三角形網(wǎng)格,共有7 267個節(jié)點,11 919個單元,如圖7所示。

圖7 鴨子蕩水庫網(wǎng)格劃分示意圖

3 水流實測結(jié)果與模擬結(jié)果對比分析

根據(jù)鴨子蕩水庫目前實際運行方式,設計工況一如下:當水庫水位降至1 247.0 m時,采取泵站抽取黃河水注入水庫的方式持續(xù)進水,日進水量為60萬m3,日取水量為40萬m3,直至水庫水位升至1 249.5 m,此時流場如圖8(a)所示。隨著寧東能源化工基地工業(yè)用水與生活用水的日益增加,設計工況二如下:當水庫水位降至1 247.0 m時進行補水,日進水量為80萬m3,日取水量為60萬m3,直至水庫水位升至1 249.5 m,此時流場如圖8(b)所示。

圖8 水流流場

分析圖8可知:隨著進水量和取水量的增加,2種工況運行下的水流流場都與實際情況相吻合;隨著進水量的增加,進水口附近的流速增加,進水口進入庫區(qū)的水流對水庫上部的流場影響增大,尤其水庫上部的主流區(qū)域和岔道區(qū)域流場變化明顯;由于水庫的中部長且淺,對進入庫區(qū)的水流起到了限制和減緩作用,因此進水量增加對水庫中部的流場影響較小,對取水塔的流場影響也較小;取水塔上下兩側(cè)形成2個回流區(qū),并且下側(cè)的回流區(qū)比上側(cè)的大;隨著取水量的增大,取水塔附近的流速增加,取水塔上下兩側(cè)2個回流區(qū)流場變化明顯,水庫下部的岔道區(qū)域流場變化也較明顯。

2014年7月和11月利用聲學多普勒流速剖面儀、RTK等儀器對水庫的17個斷面進行了實測,選擇2014年11月實測期間的數(shù)據(jù)進行數(shù)值模擬,模擬計算的初始水位為1 247.0 m,通過上述數(shù)學模型的模擬計算,經(jīng)過24 d后水位升至1 249.5 m。選擇進水口附近的第3斷面、第4斷面的右部和取水塔附近的第14斷面、第15斷面進行實測值與模擬值比較分析,如圖9所示。

圖9 斷面垂線平均流速對比

分析圖9可知:模擬值與實測值總體上比較吻合,說明水庫地形的前處理是合適的,數(shù)學模型的建立是有效的;從斷面3、斷面4右部及其他斷面的實測垂線平均流速可以看出,由于斷面右岸靠近入水口,斷面右岸附近的流速大于斷面左岸附近的流速,隨著斷面與入水口之間距離的增加,斷面右岸的流速逐漸較小,斷面左岸的流速變化很小;從斷面14、斷面15及其他斷面的實測垂線平均流速可以看出,由于斷面左岸靠近取水塔,斷面左岸附近的流速大于斷面右岸附近的流速,隨著斷面與取水塔之間距離的增加,斷面左岸的流速減小得很快,斷面右岸的流速變化很小。

4 結(jié) 語

以位于寧夏靈武市的鴨子蕩水庫為研究對象,通過圖像拼接方法得到完整的鴨子蕩水庫圖像,利用數(shù)學形態(tài)學方法和AutoCAD軟件提取鴨子蕩水庫的庫岸曲線,將提取后的經(jīng)緯度信息與實測資料相結(jié)合,構(gòu)建了鴨子蕩水庫的初始地形及斷面分布。在此基礎(chǔ)上,建立了二維水流數(shù)學模型,給出合理的邊界條件和初始條件,使用有限體積法對相關(guān)方程進行離散,對計算區(qū)域進行三角形網(wǎng)格劃分,實現(xiàn)了鴨子蕩水庫水流實測結(jié)果與模擬結(jié)果的對比分析。結(jié)果表明:圖像拼接、數(shù)學形態(tài)學等方法相結(jié)合實現(xiàn)了水庫初始地形的前處理,可以顯著提高模擬過程中的地形處理效率和精度;模型計算的二維流場分布、斷面垂線平均流速與實測資料比較吻合,說明所建立的數(shù)學模型是合理的,為下一步研究鴨子蕩水庫的泥沙運移及泥沙淤積情況奠定了基礎(chǔ)。

[ 1 ] 齊敦哲,金滿洋,趙平均,等.寧夏寧東鴨子蕩水庫漁業(yè)資源調(diào)查報告[J].水生態(tài)學雜志,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 ] 顧立忠,關(guān)琨,鄭國棟,等.淺談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 ] 關(guān)見朝,方春明.基于Google Earth的河流模擬地形前處理新方法[J].水利水電技術(shù),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的河道地形高程數(shù)據(jù)處理方法[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 ] 李春光,景何仿,呂歲菊,等.水沙水質(zhì)類水利工程問題數(shù)值模擬理論與應用[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] 陳學群,李福林,張瑞青,等.黃河河口平原多閘壩河道水流數(shù)學模型[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īng)網(wǎng)絡及其在一維河網(wǎng)模擬中的應用[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] 景何仿,李春光.黃河上游連續(xù)彎道水流運動及泥沙運移數(shù)值模擬研究[M].鄭州:黃河水利出版社,2012.

[13] 熊偉,翟劍峰,朱志夏,等.基于有限體積法的三維波生近岸流數(shù)值模型[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

國家自然科學基金地區(qū)科學基金(61362029);寧夏高等學?？茖W研究項目(NGY2015139);寧夏大學科學研究基金(ZR16018)

黃凌霄(1977—),男,副教授,博士研究生,主要從事計算水力學及數(shù)字圖像處理研究。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 編輯：高建群)