一種基于斷面高程數(shù)據(jù)的河道沖淤計算方法—曲線正交網(wǎng)格地形法*

孫中強，谷 碩，王厚杰，3，楊作升，3，畢乃雙，3**

(1.中國海洋大學(xué)海底科學(xué)與探測技術(shù)教育部重點實驗室，山東 青島 266100； 2.黃河水利委員會，水文局，河南 鄭州 450004；3.青島海洋科學(xué)與技術(shù)國家實驗室，海洋地質(zhì)過程與環(huán)境功能實驗室，山東 青島 266061)

河道沖淤演化是河流徑流與泥沙相互作用的結(jié)果，是河流穩(wěn)定與安全的關(guān)鍵指標，也是河流治理與規(guī)劃的重要參考依據(jù)。同時，作為陸源物質(zhì)與能量向海洋傳輸?shù)年P(guān)鍵通道，河道的沖淤演化也關(guān)系到關(guān)鍵元素的地球化學(xué)循環(huán)和物質(zhì)的源匯過程。目前，傳統(tǒng)的河道沖淤研究手段包括基于水文站水文泥沙觀測數(shù)據(jù)和河道地形高程測量的定量計算，即輸沙量法和斷面法[1-4]。

輸沙量法也稱輸沙量平衡法、輸沙率法，是利用水文站實測水沙資料以及區(qū)間增沙、引沙資料進行計算的[2-3]。輸沙量法不需要考慮河道地形演變，且水文站測驗時間間隔較小、資料連續(xù)，對于事件性水沙過程對河道沖淤變化的影響[5,6]具有較好的指示意義。早在1950年，Einstein就建立了均勻沙推移質(zhì)輸沙率統(tǒng)計公式[7]，1955年Colby等人提出了根據(jù)部份實測資料計算全沙輸沙率的方法[8]。然而1956年錢寧等人基于懸沙理論提出了懸移質(zhì)積點測量決定輸沙率時會引起誤差的問題[9]，此后國內(nèi)外眾多學(xué)者就河流輸沙率的修正計算進行了大量研究[2,10-14]。但是，這些計算方法多具有所需參數(shù)多、過程復(fù)雜、普適性差等缺點，在一些河道變遷頻繁或缺少歷史資料的情況下不能較好的廣泛應(yīng)用。另外，輸沙量法計算存在缺少推移質(zhì)輸沙數(shù)據(jù)、未考慮水量不平衡和區(qū)間引水引沙資料不全面等缺陷，長期計算結(jié)果存在普遍失真現(xiàn)象[15]。

斷面法又稱斷面地形法，是利用斷面地形實測資料，按照固定灘槽和間距的梯形體或錐體進行計算的[2-4]。前人研究表明，斷面法沖淤計算值僅與始末狀態(tài)有關(guān)，不存在輸沙量法的累積性誤差問題，且相對誤差隨河道累積沖淤量的增加誤差還會逐漸減小[2]；對于小范圍順直河段的計算，斷面法基本上可以反映河道實際沖淤情況[16]；而對于大范圍且河彎曲折的情況，則需在彎道、汊道、河道急劇放寬和束窄的局部河段布置相對較密的斷面[4]。由于斷面法計算僅用到斷面和斷面間距數(shù)據(jù)，其計算精度取決于實測數(shù)據(jù)的精度和斷面的密度，結(jié)果相對較粗略[17]。

隨著河道治理對于河道沖淤演化要求的提高，傳統(tǒng)方法已漸漸不能滿足實際工作中對計算精度和空間分辨能力的需求，網(wǎng)格地形法應(yīng)運而生。網(wǎng)格地形法是近年來隨著計算機技術(shù)以及遙感技術(shù)發(fā)展出現(xiàn)的一種新的河道定量沖淤計算方法，是基于衛(wèi)星遙感數(shù)據(jù)和野外實測河道地形三維數(shù)據(jù)，建立河道數(shù)字高程模型(Digital Elevation Model, DEM)進行河道沖淤計算的[4,16]。該方法不僅能夠進一步提高計算結(jié)果精度，而且擁有較強的空間分析能力，計算結(jié)果的可視化效果更加直觀，具有較好的應(yīng)用前景。

得益于近年來衛(wèi)星遙感技術(shù)的迅速發(fā)展，國內(nèi)外一些機構(gòu)(美國宇航局(NASA)、德國宇航中心(DLR)等)逐漸開放了全球陸地高程衛(wèi)星遙感數(shù)據(jù)，分辨率也相當可觀。然而，受水面波動、懸沙濃度、溫鹽變化等因素的影響，水下高程數(shù)據(jù)仍主要依靠野外實地測量獲得。河道固定斷面測量工作開展較早，發(fā)展成熟，數(shù)據(jù)相對完善，是研究地貌演化常用數(shù)據(jù)。部分國內(nèi)外學(xué)者就河道地形模型建立中的網(wǎng)格選取與生成[18]、斷面數(shù)據(jù)插值算法[19-21]、模型體積計算[17]、與遙感數(shù)據(jù)耦合[22]以及實際計算模擬應(yīng)用[1,23]等方面做了大量研究，并取得了顯著的成果。前人應(yīng)用網(wǎng)格地形法進行河道沖淤演化研究中為了保證計算結(jié)果的準確性，往往需要大量的河道地形數(shù)據(jù)作為支撐，對斷面數(shù)據(jù)的要求也較高。但是，由于野外測量實際條件的限制，進行長期、大范圍、高密度的野外測量工作并不可行。特別是早期河道斷面地形測量工作中，斷面間距、方位變化往往較大，導(dǎo)致斷面數(shù)據(jù)非常有限而且橫縱向分布極不均勻。如何利用有限的實測斷面數(shù)據(jù)結(jié)合衛(wèi)星遙感數(shù)據(jù)進行河道沖淤計算，成為現(xiàn)階段網(wǎng)格地形法急需解決的難題。

本文基于河道斷面高程數(shù)據(jù)和遙感數(shù)據(jù)，提出了一種基于DEM進行河道沖淤定量計算的普適方法，并以黃河下游部分河道為例，比較了本文方法和一般網(wǎng)格地形法、斷面法的計算結(jié)果差異，并對導(dǎo)致這一差異的主控因素進行了討論分析。

1 研究區(qū)域概況

黃河以水少沙多、水沙異源著稱，歷史上其年均向渤海輸送泥沙量高達10.8億t[24]，位于世界大河第二位。黃河下游地勢低平開闊，河道寬淺，大量泥沙淤積，使得河床抬升，形成著名的地上懸河，決口與改道頻發(fā)[25]。近年來由于自然因素以及人類活動的影響，黃河水沙量急劇減少[26]，下游河道淤積加重。為緩解下游河道淤積，提高河道行洪能力，黃河水利委員會自2002年利用小浪底等干流大型水庫在汛期開展調(diào)水調(diào)沙工作[27]。此后，下游河道持續(xù)淤積狀況改善，河床沖刷明顯[28-30]。

研究區(qū)(見圖1)位于黃河下游，本文選取了孫口-艾山河段，孫口-艾山河段位于孫口、艾山水文站之間，長約63.9 km，落差約7.5 m，河道曲折。

圖1 研究區(qū)地形及斷面分布Fig.1 Topography and settings of bathymetrical cross sections in study area

2 數(shù)據(jù)與方法

2.1 數(shù)據(jù)來源

本文計算所用數(shù)據(jù)包括河道斷面高程數(shù)據(jù)以及衛(wèi)星遙感數(shù)據(jù)。斷面高程數(shù)據(jù)為1991年10月—2012年10月黃河水利委員會統(tǒng)測黃河下游斷面高程數(shù)據(jù)(不包括1998年以后加密斷面)，孫口-艾山段共有固定斷面17個，斷面間距2～9 km。衛(wèi)星遙感數(shù)據(jù)為Landsat5的TM遙感影像數(shù)據(jù)，由美國地球資源與科學(xué)中心免費下載(http://glovis.usgs.gov/)，分辨率為30 m，其被廣泛應(yīng)用在地形地貌研究。本文選取1991—2012年每年10月前后與斷面數(shù)據(jù)測量時間最為接近的無云或少云遮擋的研究區(qū)遙感影像數(shù)據(jù)。

2.2 計算方法

2.2.1 斷面法 斷面法計算具體公式如下：

(1)

(2)

式中：為i～i+1橫斷面間河槽體積；為第i個橫斷面面積；為i～i+1橫斷面之間間距，兩測次河槽體積之差即為斷面間沖淤體積。對于任意相鄰兩斷面i、j，當且時，應(yīng)用(2)式[17]進行計算。

2.2.2 基于矩形網(wǎng)格建模的網(wǎng)格地形法 本文矩形網(wǎng)格計算基于Surfer11.0軟件完成，該軟件是美國Golden公司開發(fā)的地學(xué)常用繪圖軟件之一，其內(nèi)置多種常見插值方法，同時可以對生成的網(wǎng)格模型進行體積計算，具有方便、快捷、安裝簡單等特點，常用于數(shù)據(jù)相對充足且空間分布均勻的海洋或湖泊DEM的建立。建模具體步驟如下：(1)提取斷面測點空間三維坐標。將測點坐標用Global Mapper軟件統(tǒng)一轉(zhuǎn)換為基于高斯克呂格投影的WGS84坐標系。(2)將測點高程數(shù)據(jù)直接插值并網(wǎng)格化，得到數(shù)字高程模型。Surfer軟件直接建模生成網(wǎng)格為規(guī)則矩形網(wǎng)格，考慮斷面數(shù)據(jù)實測間隔、提取的河道邊界精度，為了使模型網(wǎng)格盡可能貼合河道邊界，計算用網(wǎng)格間距為30 m；由于斷面數(shù)據(jù)空間分布差異較大且外插數(shù)據(jù)不可靠，插值方法采用結(jié)果較為準確的三角網(wǎng)線性插值。(3)計算網(wǎng)格模型體積變化，即為河道沖淤體積。

2.2.3 基于曲線正交網(wǎng)格的網(wǎng)格地形法 本文提出了一種應(yīng)用斷面數(shù)據(jù)和衛(wèi)星遙感數(shù)據(jù)，建立基于曲線政教網(wǎng)格的DEM計算河道沖淤的方法，具體步驟如下：(1)提取研究區(qū)河道邊界。本文選擇對水陸邊界最為敏感的波段4、5、1(RGB)合成Landsat假彩色圖片，應(yīng)用Global Mapper軟件獲取歷年與斷面數(shù)據(jù)測量時間對應(yīng)的河道研究區(qū)邊界，河道固定斷面處以斷面法固定灘槽界位置控制，并勾畫研究年限內(nèi)總體河勢方向線。(2)根據(jù)河道邊界建立DEM網(wǎng)格?？紤]河道蜿蜒、曲折多變，規(guī)則矩形網(wǎng)格不能準確表示河道邊界的變化(見圖2a)，本文選用與河道邊界更為貼合的曲線正交網(wǎng)格(見圖2b)。網(wǎng)格間距一般不超過30 m，最大為100 m。(3)斷面數(shù)據(jù)插值建模。應(yīng)用自編Matlab程序首先按照反距離插值法將實測斷面高程數(shù)據(jù)插值到垂直于主河道方向的網(wǎng)格控制點，再沿河道河勢方向按河長逐行進行線性內(nèi)插，得到河道DEM(見圖3)。(4)模型體積計算。使用自編Matlab程序計算模型空間體積，模型體積變化即為河道沖淤變化體積。如圖4，DEM中x-y平面任意三個相鄰點(A，B，C)對應(yīng)的五面體(ABC-DE′F′)的體積等于三棱柱(ABD-DEF)的體積減去五面體(D-EFF′E′)的體積，即

(3)

其中：AD為最長棱，可通過排序賦值得到；DD′為D到EF的高。

((a)規(guī)則矩形網(wǎng)格Rectangular grid; (b)曲線正交網(wǎng)格Orthogonal curvilinear grid.)
圖2 不同網(wǎng)格貼合河道邊界示意圖
Fig.2 Two different DEM grids

圖3 基于曲線正交網(wǎng)格的斷面數(shù)據(jù)插值示意圖Fig.3 Schematic diagram of interpolationin calculation domain based on orthogonal curvilinear grid

圖4 DEM體積計算示意圖Fig.4 Schematic diagram for volume calculation of DEM

3 不同模型沖淤計算結(jié)果

孫口-艾山河段主河道的沖淤計算結(jié)果表明三種方法計算結(jié)果沖淤性質(zhì)基本一致，但沖淤數(shù)量存在明顯差異，尤其是矩形網(wǎng)格建模與本文方法和斷面法差異顯著(見圖5)。研究區(qū)主河道整體呈先淤后沖趨勢：除1993、1996和1998年出現(xiàn)沖刷外，至2000年10月研究區(qū)呈持續(xù)淤積，單次最大淤積出現(xiàn)在1997年；此后研究區(qū)呈持續(xù)沖刷，2004年10月與1992年10月沖淤水平基本相當，2012年10月達最大沖刷，單次最大沖刷出現(xiàn)在2003年。

結(jié)合相關(guān)性分析(見圖6)和誤差分析(見表1)可以看出，矩形網(wǎng)格建模和斷面法的累計誤差隨累計沖刷量增大而增大，矩形網(wǎng)格建模計算誤差更大。矩形網(wǎng)格建模較本文方法計算結(jié)果明顯偏小，兩者擬合直線斜率僅為0.62，線性正相關(guān)性較好；兩者絕對誤差范圍為(-5.91～4.76)×106m3，平均絕對誤差約為0.68×106m3，相對誤差范圍可達-157.3%～56.7%，平均相對誤差高達-40.7%。斷面法計算結(jié)果總體較本文方法略有偏大，兩者線性擬合斜率為1.09，呈高度線性正相關(guān)，但絕對誤差范圍為(-1.05～1.35)×106m3，平均誤差為-0.15×106m3，相對誤差范圍高達-71%～45%，平均相對誤差約為9.0%。

(正值表示淤積，負值表示沖刷。Positive values indicate accumulation, and negative values indicate erosion.)
圖5 孫口-艾山主河道累計沖淤曲線
Fig.5 Cumulative erosion-accumulation sediment volume in main channel between bathymetrical cross sectionsSunkou and Aishan

圖6 計算結(jié)果相關(guān)性分析Fig.6 Correlation analysis of calculated erosion-accumulation sediment volumes using different methods

表1 1991年10月～2012年10月孫口-艾山段主河道沖淤體積計算結(jié)果Table 1 Erosion-accumulation sediment volume in main channel from October 1991 to October 2012 between bathymetrical cross sections Sunkou and Aishan

注：1絕對誤差=斷面法或矩形網(wǎng)格法-本文方法；2相對誤差=(斷面法或矩形網(wǎng)格法-本文方法)/本文方法。
Note:1.Absolute error=Gross-section or Method using rectangulargrid-Method in this paper; 2.Relative error=(Gross-section or Method using rectangulargrid-Method in this paper)/Method in this paper.

4 討論

本文方法充分考慮了河床地形的復(fù)雜多變情況，從衛(wèi)星遙感圖片提取河道邊界，使用更好適應(yīng)性的曲線正交網(wǎng)格，能夠很好的重構(gòu)復(fù)雜多變的河道形態(tài)。對于河道地形變化橫向插值時僅使用插值點最近兩個點的斷面數(shù)據(jù)，保證了插值結(jié)果與斷面數(shù)據(jù)的高度一致，降低了斷面方位變化對插值結(jié)果的影響(見圖7)。此外，用曲線正交網(wǎng)格的河流中泓河勢方向作為模型的插值方向，符合河道地形的縱深變化，采用線性內(nèi)插進行縱向插值，保證了網(wǎng)格點數(shù)據(jù)只受上下兩斷面數(shù)據(jù)影響，能夠利用有限的斷面數(shù)據(jù)較為準確重構(gòu)河道地形(見圖8a)。

(a.垂直斷面(邵莊)Perpendicular section(Shaozhuang)； b.斜交斷面(孫口)Oblique section(Sunkou).)

4.1 矩形網(wǎng)格建模的誤差原因分析

矩形網(wǎng)格建模與本文方法誤差較大，其原因是河道斷面高程數(shù)據(jù)空間分布極不均勻，相鄰兩固定斷面間隔較大(2～9 km)，數(shù)據(jù)點對網(wǎng)格控制明顯不足(見圖8b)。直接插值計算地形高程與距離的關(guān)系時，使用的是點與點之間的直線距離L1(見圖8b)。實際上，河道地形高程是點與點沿平行或垂直主河道主流河勢方向距離河長L2的函數(shù)，對于河彎曲折的河段L1明顯小于L2，使插值結(jié)果與實際河道情況不符，得到的河道地形模型有明顯的不連續(xù)條帶。此外，在對河道邊界外部區(qū)域的插值數(shù)據(jù)進行剔除的過程中，只有部分DEM數(shù)據(jù)參與計算，會導(dǎo)致結(jié)果進一步偏小。因此，導(dǎo)致基于矩形網(wǎng)格的計算的河道沖淤量顯著小于斷面法和基于曲線正交網(wǎng)格法(見圖5)。

(a.本文方法建模This study； b.矩形網(wǎng)格建模Model using rectangulargrid.)圖8 不同方法建模結(jié)果高程圖Fig.8 Modeling resultsbased on different methods

4.2 斷面法的誤差原因分析

斷面法與本文方法計算結(jié)果總體較為接近，孫口-艾山河段斷面法計算河道沖淤量較本文曲線正交網(wǎng)格地形法平均高約9%，但個別年份相差可達71%(1997—1998年)(見表1)。斷面法主要基于斷面面積變化和斷面間距進行沖淤計算，并假設(shè)斷面間河道寬度無明顯變化以及斷面方位與河道基本垂直。然而，受自然和人為因素的影響，實際河道寬度沿程變化顯著且不均勻和隨機性較強，河道寬度最小約為250 m，最大可達650 m(見圖9)。斷面法僅使用有限的斷面測點高程數(shù)據(jù)以及斷面間距數(shù)據(jù)，將河道地形簡單的使用梯形體或者截錐表示，其結(jié)果強烈依賴于斷面代表性，存在一定的偶然性，無法準確精細的刻畫河道地形的演變，導(dǎo)致斷面法與本文的曲線正交網(wǎng)格法存在一定的差異。此外，由于黃河徑流量和泥沙含量年際和季節(jié)性變化較大，自然狀態(tài)下的河道易發(fā)生擺動，使斷面方位與河流流向的夾角不斷變化，這一現(xiàn)象在研究區(qū)普遍存在(見表2)。當觀測斷面與河流流向存在一定夾角時(，則觀測斷面長度、面積與實際斷面長度和面積存在如下關(guān)系：

圖9 研究區(qū)河道寬度變化Fig.9 Variation in the width of the river channel

(4)

其中：L為觀測斷面長度；l為實際斷面長度，對應(yīng)實際斷面面積也為觀測斷面面積的sinα倍。斷面法計算時沒有改正角度引起的斷面面積偏大，使得計算河槽體積也相應(yīng)偏大，沖淤量也相應(yīng)偏大(見表3)，這是導(dǎo)致斷面法較本文方法偏大的另一因素。

隨著河道治理與監(jiān)測的不斷發(fā)發(fā)展，河道斷面密度不斷增加。1998年以后黃河水利委員會對黃河下游淤積斷面觀測普遍加密，斷面間距進一步減小，至2012年西霞院以下河段共有觀測斷面約370個，最小斷面間距僅450 m。此外，隨著科學(xué)技術(shù)的發(fā)展，野外測量的精度和遙感影像分辨率也在不斷提高。這將會進一步提高構(gòu)建河道DEM的精度，推進河道地形地貌演化研究及其機理研究和河道綜合治理。

表2 研究區(qū)觀測斷面與河流流向夾角表Table 2 Average angle between bathymetrical cross section and the river channel in study area during 1991—2012

表3 1991年10月～2012年10月孫口-艾山斷面與流向夾角校正結(jié)果Table 3 Angle correction resultsfrom 1991 to 2012 between cross sections Sunkou and Aishan

5 結(jié)論

(1)本文提出了一種應(yīng)用河道固定斷面測深資料和衛(wèi)星遙感數(shù)據(jù)，基于曲線正交網(wǎng)格的河道數(shù)字高程模型(DEM)定量計算河道沖淤的普適方法，該方法不受河道形態(tài)和斷面方位等因素的限制，能夠利用有限的斷面數(shù)據(jù)較為準確的構(gòu)建河道DEM進行沖淤計算。

(2)以黃河孫口-艾山河段為例，分別應(yīng)用本文方法、基于矩形網(wǎng)格的網(wǎng)格地形法和斷面法進行沖淤計算，結(jié)果表明：基于矩形網(wǎng)格建模和斷面法計算結(jié)果與本文結(jié)果趨勢基本一致，但均存在明顯誤差。斷面法計算結(jié)果與本文方法結(jié)果較為接近，平均相對誤差約為9%，個別年份可達70%以上；基于矩形網(wǎng)格建模計算結(jié)果與本文方法結(jié)果相差較大，平均相對誤差可達40%以上，最大誤差超過157%。

(3)現(xiàn)有斷面數(shù)據(jù)往往分布不均勻，對于彎曲河道的控制明顯不足，是矩形網(wǎng)格建模計算誤差的主要原因。而斷面法計算誤差主要由于未考慮河道寬度的沿程變化，僅基于斷面面積和斷面距離的梯形體或截錐計算河道沖淤變化。另外，由于河道長期擺動，斷面與河流流向不垂直導(dǎo)致測量斷面面積偏大，也是導(dǎo)致計算誤差的主要原因。