1、第39卷第3期2006年6月武漢大學學報(工學版)EnineerinournalofWuhanUniversitggJyVol.39No.3June2006文章編號:()1671-8844200603-016-05基于MATLAB三維數字地形圖在庫區沖淤分析中的應用孫東坡,廖小龍,王鵬濤,孫 羽(華北水利水電學院,河南鄭州 4)50011摘要:根據現狀大多數水庫庫區水下地形復雜、缺乏詳盡地形資料的特點,提出了基于MATLAB的兩種地形高程插值方法,并闡述了綜合運用這兩種插值方法制作庫區三維數字地形圖及誤差分析方法.通過應用實例分析結果表明,這種方法適合現狀水庫地形測量資料,操作簡單,程序實現容

2、易,所得三維數字地形圖可以滿足工程上對水庫沖淤分析計算的精度要求.關鍵詞:庫區;三維數字地圖;三角單元法;四邊形單元法MATLAB;中圖分類號:TV145.1 文獻標識碼:AAlicationofthree-dimensionaldiitalmaoanalsisofppgptyreservoirserosionanddeositionbasedonMATLABp,WANGP,SUNDonoLIAOXiaolonentaoSUNYugpgg(,)NorthChinaInstituteofWaterConservancndHdroelectricPowerZhenzhou450011,Chinay

3、ayg:AAbstractsmostofreservoirlandformisveromlexandwearelackofthedetailedinformationa-ycpboutthem,sothispaerreferredtwoinsertinaluemethodsbasedonMATLAB,andexlainedhowpgvptoalhesemethodstomakethediitalreliefmanerroranalticalmethodalsoisreferred.Theppytgp.Ayalicationinstanceindicatesthatthesemethodsare

4、suitableforthepresentsituationinmostofreser-ppvoirsinChina.Theresimleratedandcanbeeasiloneinproram.Thethree-dimensionaldi-yapyopydggitalreliefmaadebhesemethodsissatisfiedtheprecisionforanalzinhereservoir'serosionandpmytygtdeosition.p:;Keordsreservoirthree-dimensionaldiitalreliefmaMATLAB;trianula

5、relementmethoduad-gp;gqywranleelementmethodg水庫清淤是延長水 在水資源日益短缺的今天,庫壽命、提高水庫攔洪蓄水能力的有效措施.了解水庫淤積數量及分布規律,摸清清淤區域和深度,對確定清淤方案是至關重要的.水庫淤積分析的傳統處理方法是通過不同時期的實測大斷面對比分此方法簡單易析和計算來確定水庫淤積分布狀況.行,但分析計算成果的精度受庫區分段布局合理性且只能分析庫區沖淤的平面分及疏密程度的影響,布規律,難以確定庫區沖淤的空間分布規律.如果能方便取得庫區不同時期準確的水下數字地形資料,則不僅可以準確計算水庫不同時期內的沖淤量,而且可以了解庫區沖淤量空間分布規

6、人們過去多是直接在庫區普通地形圖上進行人律.工網格剖分提取地形數據,而今可采用比較合理的收稿日期:2005-09-28作者簡介:孫東坡(),男,河南開封市人,教授,主要從事環境水利及河流動力學教學與研究.1952-基金項目:國家自然科學基金()資助項目.50279009第3期孫東坡等:基于MATLAB三維數字地形圖在庫區沖淤分析中的應用17插值方法,通過計算機自動剖分提取數據來制作數字地形圖,并進行庫區沖淤分析.MATLAB具有很強的數據處理功能、靈活的插值方法、強大的圖是較理想的庫區三維數字地形像顯示和處理能力,制作軟件工具.本文針對庫區水下地形特點、立足工程實際情況,闡述如何使用MATLA

7、B實現水庫淤積數量及淤積分布規律的計算方法.資料,首先用測量學閉合差法刪除不合理數據,然后將數據整理為計算機程序易讀取的向量形式:X,Y,Z.1.2.2 圖紙數據整理七、八十年代的大多數地形資料用圖紙存儲.即使是現在,工程上常常得到的地形資料也是圖紙.對這類數據的處理方法:首先掃描轉為電子圖像,然后用矢量軟件轉為C最后按電子地AD圖,1 利用MATLAB處理水下地形的方法.1 MATLAB在地形處理方面的相關函數MathWorks公司開發的高性能工程計算軟件MATLAB,主要包括核心函數和工具箱.在地形處理方面,常用函數有:interp1,interp2,gridda-a,mesh和surf1

8、.(1)yi=interp1(x,y,xi,method)是用于確定邊界點位置的一維插值函數,x,y為已知坐標向量,xi為需插值坐標向量,method為插值方法.MATLAB提供了4種插值方法:nearest,linear,pline和cubic.(2)zi=interp2(x,y,z,xi,yi,method)是根據有結構四邊形網格數據點插值的二維插值函數,y,z為已知坐標矩陣,且x和y為單調向量,xi,i為需插值坐標矩陣,zi為插值結果.若已知矩陣形式的地形數據,則此函數是比較理想的插值方法.(3)zi=griddata(x,y,z,xi,yi,method)是根據離散數據點插值的二維插值

9、函數,其變量含義與nterp2類似,但x和y可以為非單調向量,因而可進行離散數據點插值.此外griddata提供了“V4”插值方法.若已知較密集且分布較均勻的數據點,用“V4”插值可得到準確、光滑地形.(4)mesh(x,y,z)和surf(x,y,z)為常用的地形重現函數,mesh函數用于生成反映網格布設的網眼圖,surf函數用于生成表現地形形態的三維色帶圖.2 數據前處理庫區地形測繪技術從人工測繪發展到GPS及遙感測繪,前后數據形式的差別較大.從目前數據形式看,大致可以分為:零散人工觀測數據、圖紙數據和電子地圖.對這三類形式數據處理方法如下:.2.1 零散人工觀測數據整理人工觀測數據是分析

10、庫區早期沖淤的第一手圖方法處理.因圖紙常常收縮變形,而矢量化過程又是按像素進行的,所以需要校正CAD圖坐標.1.2.3 電子地圖數據整理常遇到的電子地圖是CAD圖,地形高程常用高程點和等高線表示,而這兩類數據容易用CAD的二次開發工具提取.這里給出提取以TEXT對象表示高程點的VBA程序核心片段:textset.SelectacSelectionSetAll,fType,fDa-ta創建數據集textsetForEachadTextIntextsetIfIsNumeric(adText.TextString)Then判斷是否是數值型數據IfVal(adText.TextString)>=

11、ZminAndVal(adText.TextString)<=ZmaxThen判斷是否滿足高程要求Print#1,adText.InsertionPoint(0),adText.InsertionPoint(1),adText.Text-String1.3 插值方法數據整理過程沒改變數據點空間分布的無結構性,而現在所有數字圖都基于有結構網格,所以需要數據插值.雖然MATLAB提供了兩個空間插值函數,但它們的插值方法都是對整個區域插值,且在已知數據點較稀區域常無法插值和高程點突變,因此需要根據地形特點尋求合適的插值方法.庫區水下地形及其沖淤變化與水流條件直接相關,一般在水庫近岸區和上游進口

12、區的水位及流速變化較大,沖淤幅度大,水下地形相對陡峻;而在庫區中心深水區,沖淤主要受異重流與細沙沉積影響,地形平緩.根據庫區地形以上特點,可以采用三角單元法和四邊形單元法分別處理,操作流程如圖1所示.1.3.1 三角單元法三角單元法是針對淺水沖淤地形提出的,其基本假定:待插值區域地形是由已知數據點系列確定1tsxyi1118武漢大學學報(工學版)2006一系列已知節點誤差平方和最小的四邊形平面組成.由最小二乘法理論可知,節點誤差平方和最小的四邊形平面就是4個節點按最小二乘法確定的2平面按最小二乘法原理確定四邊形平面和插.值,步驟如下:圖1 插值流程首先,生成四邊形單元.四邊形單元可以通過合并三

13、角單元獲得.為了使4個節點對其構成區域內各插值點具有最強的約束力,所以合并原則是:的一組三角面組成.根據基本假定,落在每一個三角形內待插值點高程的確定步驟如下:首先,生成三角單元網格系.根據整理得到已知數據點坐標向量X,Y,Z和MATLAB提供的三角單元生成函數delaunay(X,Y),就可以方便地生成三角單元網格系.其次,確定三角單元形函數.對任意一個三角單元的3個節點Pi(xi,yi,zi),Pj(xj,yj,zj),Pk(xk,yk,zk),根據前面假定,它們唯一確定一個平面S:z=ax+by+c,寫成矩陣形式為:z=(x,y,1)·N,其中N=(a,b,c)T為形函數矩陣,

14、是節點的坐標函數.將節點坐標代入平面矩陣方程,整理得:AN=C(2)xiyi1i其中,A=xjyj1,C=jxkyk1k因為平面S唯一存在,所以方程組(2)有唯一解,所以有N=A-1C(3)MATLAB提供了線性方程組的求解功能,方程組(2)的求解形式為:N=A/C.最后,三角單元內插值.設空間三角形平面Sm在xoy面的投影域為Dm,所以對xoy平面上任意一點(xn,yn).Dm,都有zn=(xn,yn,1)·Nm,空間點Pn(xn,yn,zn)唯一確定.從本質上說,三角單元法仍屬于線性插值方法,但它是對三角單元區域而不是整個插值區域插值,所以操作靈活,插值速度較快.3.2 四邊形單

15、元法從理論上說,在水庫深水區仍可用三角單元法插值.但有兩個缺點:首先若進行水流模擬,對比較穩定的深水區水流,三角形網格插值所得地形顯得有些粗糙;其次由圖紙獲得數據不可避免存在誤差,而三角單元法不能合理調整誤差,使插值所得地形總誤差最小.為解決上述兩個問題,在深水區宜采用四邊形單元法,該方法的基本假定是:插值區域地形是由合并所得四邊形具有最小外接圓.其次,確定節點誤差平方和最小的四邊形平面.對任意一個四邊形平面,其矩陣方程為:z=(x,y,1)·N,其中N=(a,b,c)T為形函數矩陣,是節點的坐標函數.將四邊形節點Pi(xi,yi,zi),Pj(xj,yj,zj),Pk(xk,yk,

16、zk)和Pl(xl,yl,zl)代入矩陣方程,整理得:AN=C(4)xiyi1zi其中,A=xjyj1xkyk1,C=jkxlyl1zl顯然線性方程組(4)為不相容方程組,其最小二乘法解為N=(ATA)-1ATC(5)MATLAB提供了求解線性不相容方程組最小二乘法解功能,其求解形式:N=A/C.最后,四邊形單元內插值.設空間四邊形平面Sm在xoy面的投影域為Dm,所以對xoy平面上任意一點(xn,yn).Dm,都有zn=(xnyn1)·Nm,空間點Pn(xn,yn,zn)唯一確定.從上面的求解過程可以看出:(1)由于四邊形單元法采用最小二乘法求解,能比較有效地限制誤差數據點影響;(

17、2)四邊形單元法只適合地勢變化較平緩的區域.4 邊界處理如果庫區邊界不是十分復雜,能夠由已知數據點完全確定,采用三角單元法和四邊形單元法插值時,則邊界可以不用另作處理.但在進行庫區沖淤分析時,要求相對比兩個時段的沖淤地形有共同參照控制邊界.這時就必須另加邊界線,邊界處理的方法是添加邊界分段函數y=f(x).若采用MAT-LAB程序,邊界離散處理方法是:根據采集的邊界特征位置數據點,應用interp1函數,采用分段插值,標定各邊界點.11第3期孫東坡等:基于MATLAB三維數字地形圖在庫區沖淤分析中的應用191.5 插值結果誤差分析三角單元法和四邊形單元法都能有效地將插值點高程控制在已知數據點高

18、程范圍內,但還需了解插值精度.下面從統計學角度出發,通過樣本誤差分析來評定總體誤差.從已知數據集中隨機抽取少量數據點組成樣本集S而剩余數據組成插值集S以樣本集Sa,i.a中的數據點為待插值點,對插值集Si作用插值方法f,則可生成與樣本集S這a相對應的結果集Sr.析時,采用了傳統的庫容分析法、大斷面分析計算法進行不同時期庫區淤積狀況分析,同時還采用了基于MATLAB的庫區數學模型計算法進行不同時期庫區淤積總量及淤積分布狀況計算.為了檢驗庫區數學模型的可信度,對采用的三種分析方法的4成果進行合理性對比分析.庫區數學模型計算法的核心是:建立不同時期的庫區三維數字地圖,利用MA確定庫區TLAB的功能,

19、沖淤在空間的分布,計算各個時期總體沖淤量和各高樣,對任意一個數據點Pi(xi,yi,zi).Sa,對集合r都有唯一的確定點P'i(xi,yi,z'i)與Pi對應,則插值誤差為:Xi=zi-z'i,則樣本集Sa插值誤差的平均值/X和標準差S為n/X=nXi(6)i=1S=nn-1(Xi-/X)2(7)i=1由統計規律可知:當樣本較大時,樣本誤差將近似服從零均值正態分布3;由此根據/X和S就可評定樣本插值結果的合理性和誤差離散程度,進而估計總體插值結果誤差情況.應用實例對某中型水庫的淤積特點和發展趨勢進行分圖2 某水庫1979年庫區地形圖 結合三維數字地形和沖淤分布函數,還

20、可以計算出庫區某高程以下的沖淤量,做出淤積量與高程關系曲線,圖4為不同年份庫區淤積量與高程關系曲線.分析淤積量與高程關系曲線變化特點,可以了解不同時期庫區淤積的空間分布特點:可以看到自1987年以后高程在160m以下的庫區幾乎沒有新的淤積,主要淤積都集中分布在165m以上高程的庫區.基于MATLAB的庫區數學模型可以提供詳盡的淤積分布態勢,對預估庫區淤積發展趨程段沖淤量,從而分析水庫淤積規律和發展趨勢.采用庫區數學模型計算法,首先要建立庫區地形差值模型;然后按照前述數字地形圖制作方法做出該水庫3個典型年份(1979、1987、2001年)的庫區數字地形.這里僅給出1979年庫區三維地形圖,如圖

21、2所示.有了3個典型年份的庫區數字地圖,由公式(8)就可以得到相鄰兩年份的沖淤分布函數Zij(x,y),根據沖淤分布函數利用編程MAT-AB的功能,就可以進一步得到庫區沖淤量分布,同時根據公式(9),可以計算出相鄰兩年份間庫區沖淤總量Vij.圖3為19791987年期間某水庫庫區的沖淤分布圖.Zij(x,y)=Zi(x,y)-Zj(x,y)(8)Vij=0Zij(x,y)ds(9)D圖3 某水庫1979-1987年淤積分布圖勢、采取水庫減淤應對措施都是十分重要的.為了驗證數字地形圖的可靠性,可以采用兩種方法:一是前面介紹的誤差分析方法,二是不同淤積分析方法成果對比.采用誤差分析方法時,將已知數

22、據點隨機分為20組,然后對每組數據進行誤差分析,圖為1979年數字地形的誤差分析成果,表明總體誤差較小,插值的數字地形是可靠的.采用3種淤積分析方法計算19792001年期間某水庫庫區總淤積量,具體成果統計見表1.通SL2520武漢大學學報(工學版)2006圖4 某水庫不同時期庫區淤積與高程關系曲線 圖5 某水庫1979年數字地圖誤差分析成果過對比分析3種計算結果可知:庫區數學模型計算度;法成果恰在庫容分析法計算成果與大斷面分析計算法之間;相對誤差很小.一般講庫容分析法計算法與大斷面分析計算法是基于實測資料的常規方法,相對比較可靠,因而可以認為利用數字地形計算的淤積成果總體上也是準確可信的.表1 兩個不同時期庫區總淤積量統計104m3計算方法庫容斷面數學分析法分析法模型法備 注計1979-1987年107713001210相對誤差算相對誤差/%9