非連續(xù)河道地形插值方法的比選

華祖林,王海燕,汪 靚,王玉琳（1.河海大學(xué)淺水湖泊綜合治理與資源開發(fā)教育部重點(diǎn)實(shí)驗(yàn)室,江蘇南京 210098；2.河海大學(xué)環(huán)境學(xué)院,江蘇南京 210098；3.河海大學(xué)水資源高效利用與工程安全國(guó)家工程研究中心,江蘇南京 210098）

?

非連續(xù)河道地形插值方法的比選

華祖林1,2,3,王海燕1,2,汪 靚1,2,王玉琳1,2
（1.河海大學(xué)淺水湖泊綜合治理與資源開發(fā)教育部重點(diǎn)實(shí)驗(yàn)室,江蘇南京 210098；2.河海大學(xué)環(huán)境學(xué)院,江蘇南京 210098；3.河海大學(xué)水資源高效利用與工程安全國(guó)家工程研究中心,江蘇南京 210098）

摘要:針對(duì)在許多河道水動(dòng)力及水環(huán)境模擬計(jì)算過(guò)程中僅有一定間距下的斷面地形數(shù)據(jù)的實(shí)際情形,以黃浦江為例,分別采用三角網(wǎng)格法、反距離權(quán)重法和克里金法3種方法對(duì)強(qiáng)彎段、干支流交匯處及復(fù)式分汊河道的水下地形進(jìn)行插值,以對(duì)比分析3種方法的適應(yīng)性。結(jié)果表明,對(duì)于強(qiáng)彎段、干支流交匯處及復(fù)式分汊河道等地形起伏變化較大區(qū)域,三角網(wǎng)格法會(huì)造成插值結(jié)果不穩(wěn)定,誤差及波動(dòng)大,未能反映部分地形起伏實(shí)際情況,反距離權(quán)重法可能出現(xiàn)“牛眼”現(xiàn)象,克里金法表現(xiàn)出較好的適應(yīng)性;在減少測(cè)量斷面時(shí),克里金法的插值結(jié)果誤差最小,強(qiáng)健性較好;將3種插值方法獲得的黃浦江全河段插值地形用于整體潮位計(jì)算,克里金法計(jì)算精度最高。

關(guān)鍵詞:非連續(xù)河道;地形插值;三角網(wǎng)格法;反距離權(quán)重法;克里金法;黃浦江

在水動(dòng)力及水環(huán)境數(shù)值模擬中,計(jì)算網(wǎng)格上地形值能否反映真實(shí)地形非常重要,會(huì)直接影響計(jì)算精度。對(duì)河道來(lái)講,在許多實(shí)際情況下沒有整個(gè)河道的連續(xù)地形數(shù)據(jù),只有一定間距下的斷面地形數(shù)據(jù),數(shù)值模擬計(jì)算時(shí)需要通過(guò)插值方法由斷面數(shù)據(jù)得到網(wǎng)格地形值,尤其河道轉(zhuǎn)彎或河流交匯等起伏變化很大處的地形,因此插值方法是否適用極為重要。

插值方法有很多,在進(jìn)行空間地形插值時(shí),構(gòu)造地形插值曲面的常用方法主要有三角插值和距離加權(quán)插值兩類,目前使用最為廣泛的是三角網(wǎng)格法、反距離權(quán)重法及克里金法,其中三角網(wǎng)格法為三角插值法,反距離權(quán)重法和克里金法為距離加權(quán)插值法。很多水動(dòng)力計(jì)算軟件如Mike使用的地形插值方法為三角網(wǎng)格法,三角網(wǎng)格法插值速度快,更適合于大面積且均勻化分布數(shù)據(jù)的網(wǎng)格化；Delft3D使用的是三角網(wǎng)格法及反距離權(quán)重法,EFDC使用的是反距離權(quán)重法,反距離權(quán)重法可以通過(guò)權(quán)重調(diào)整空間插值等值線的結(jié)構(gòu)。克里金法也是一種常用的空間插值方法,具有很好的空間相關(guān)性,其結(jié)果是一種無(wú)偏最優(yōu)估計(jì)量,適于各種類型的離散數(shù)據(jù)。這3種方法在地形插值中得到了廣泛的應(yīng)用,如何磊等[1]利用三角形網(wǎng)格法得到了珠海的三維地形圖；陳中等[2]利用改進(jìn)的三角形網(wǎng)格法對(duì)某河道地形進(jìn)行了地形插值；劉兆平等[3]用克里金法、三角形網(wǎng)格法和反距離權(quán)重法分別繪制了河北省某市郊區(qū)垃圾填埋場(chǎng)的等值線斷面圖,得到了垃圾填埋廠的垃圾分布情況；牛文杰[4]運(yùn)用克里金法和改進(jìn)的反距離權(quán)重法對(duì)某礦區(qū)煤層數(shù)據(jù)及某油田地震測(cè)線數(shù)據(jù)分別進(jìn)行了插值,得到煤層分布及地層深度等值線圖；Detweiler等[5]利用反距離權(quán)重法和克里金法計(jì)算得到了三維地形；李曉軍等[6]借鑒地質(zhì)統(tǒng)計(jì)學(xué)的Kriging方法給出了一種距離加權(quán)插值算法；張殿新等[7]用改進(jìn)的克里金法得到了珠海港的網(wǎng)格數(shù)字地形；包世泰等[8]利用克里金法對(duì)廣州市南沙區(qū)進(jìn)港大道進(jìn)行了地形高程插值；杜國(guó)明等[9]用克里金法得到了珠江河道地形陰影圖；曾懷恩等[10]利用改進(jìn)變異函數(shù)估計(jì)的方法對(duì)某次地形測(cè)量的數(shù)據(jù)進(jìn)行了插值。

但是,有關(guān)不同方法的比較和適合性研究成果還較少,由于通過(guò)插值得到的網(wǎng)格地形數(shù)據(jù)對(duì)數(shù)值計(jì)算結(jié)果有重要影響,尤其是河道地形起伏較大的地方,可通過(guò)比較選擇適應(yīng)性更好的插值方法來(lái)提高數(shù)值模擬的精度。本文采用上海黃浦江干流地區(qū)的實(shí)測(cè)地形數(shù)據(jù),選取強(qiáng)彎段、干支流交匯處及復(fù)式分汊河道3種類型的地形,分別用三角網(wǎng)格法、反距離權(quán)重法和克里金法進(jìn)行插值和對(duì)比分析。

1 空間插值方法及數(shù)據(jù)來(lái)源

1.1 空間插值方法的選取

三角網(wǎng)格法、反距離權(quán)重法和克里金法是水動(dòng)力模型中常用的插值方法,本文選用這3種方法進(jìn)行對(duì)比分析。

三角網(wǎng)格法是通過(guò)空間坐標(biāo)與水下地形高程的關(guān)系式求得未知點(diǎn)高程值的,即求解如下空間上i個(gè)點(diǎn)組成的方程式得到任意未知點(diǎn)高程:

反距離權(quán)重法和克里金法均是通過(guò)權(quán)重估計(jì)待測(cè)點(diǎn)的值（計(jì)算公式見式（2）[11]）,但不同的是,反距離權(quán)重法的權(quán)重僅與距離有關(guān),而克里金法則同時(shí)考慮了距離及空間相關(guān)性,通過(guò)變異函數(shù)[12]來(lái)刻畫區(qū)域化變量的空間結(jié)構(gòu)性變化和描述其隨機(jī)性變化。

式中:z為估計(jì)值；zi為已知點(diǎn)高程值；wi為權(quán)重；n為空間內(nèi)已知坐標(biāo)點(diǎn)的個(gè)數(shù)。

1.2 數(shù)據(jù)來(lái)源

采用的黃浦江地形數(shù)據(jù)為斷面數(shù)據(jù),干流斷面間距約0.5 km,每個(gè)斷面上有150～200個(gè)實(shí)測(cè)地形數(shù)據(jù),支流斷面間距約0.4km,每個(gè)斷面上約有50～100個(gè)實(shí)測(cè)地形數(shù)據(jù)。黃浦江潮位數(shù)據(jù)為河祝、金澤、東團(tuán)及吳淞口共4個(gè)站點(diǎn)的實(shí)測(cè)數(shù)據(jù),每5min一個(gè)值,數(shù)據(jù)均由上海市水文總站提供。

2 不同河道形態(tài)的地形插值結(jié)果比較

對(duì)于水動(dòng)力及水環(huán)境模型計(jì)算而言,對(duì)河道斷面地形量測(cè)的間距是需要根據(jù)國(guó)家規(guī)范要求布設(shè)的,根據(jù)SL257—2000《水道觀測(cè)規(guī)范》,黃浦江實(shí)測(cè)地形斷面間距為0.4～0.5 km。本文采用黃浦江干流的實(shí)測(cè)地形數(shù)據(jù),分別選取強(qiáng)彎段、干支流交匯處（吳淞口）及上游分汊河段3種不同類型地形（位置關(guān)系見圖1）,用三角網(wǎng)格法、反距離權(quán)重法和克里金法進(jìn)行插值。

圖1 研究河段位置示意圖

2.1 強(qiáng)彎段

考慮到在轉(zhuǎn)彎處河道受凹凸不同的水流作用,會(huì)產(chǎn)生不同于順直河道的特殊地形,選取強(qiáng)彎段7個(gè)斷面進(jìn)行插值,3種方法等值線插值結(jié)果如圖2所示。

從圖2可以看出,三角網(wǎng)格法插值結(jié)果從某個(gè)平穩(wěn)高程陡變?yōu)榱硪粋€(gè)高程,結(jié)果不夠光滑,且忽略了地形的部分起伏變化。反距離權(quán)重法出現(xiàn)“牛眼”現(xiàn)象,即很多計(jì)算值把實(shí)測(cè)值包圍在中間的情況,這表明反距離權(quán)重法的權(quán)重僅與距離有關(guān),沒有空間方向相關(guān)性,不能判斷待測(cè)點(diǎn)的方位,對(duì)于河道等狹長(zhǎng)型地形來(lái)說(shuō),極容易將邊界處與河道中心等權(quán)處理,尤其在轉(zhuǎn)彎處無(wú)法判別有效待測(cè)點(diǎn),使插值結(jié)果失真?？死锝鸱▌t克服了“牛眼”這一缺點(diǎn),其最大的特點(diǎn)就是可以通過(guò)改變變異函數(shù)的結(jié)構(gòu)及參數(shù)從不同的角度反映空間變異性,解決空間連續(xù)變化的屬性非常不規(guī)則時(shí)數(shù)據(jù)的插值問(wèn)題,從而對(duì)空間變量的連續(xù)性、相關(guān)性、影響范圍、間斷性等要素進(jìn)行描述,因而具有很好的適應(yīng)性,對(duì)于不同復(fù)雜程度的地形,插值結(jié)果都較好。

圖2 強(qiáng)彎段等值線插值結(jié)果

2.2 干支流交匯處

吳淞口位于黃浦江與長(zhǎng)江匯流處,干支流在交匯處發(fā)生強(qiáng)烈交互,水流能量損失大,容易在不同時(shí)期表現(xiàn)出不同的地形形態(tài)。因此,選取吳淞口處3個(gè)斷面數(shù)據(jù),用相同方法對(duì)地形數(shù)據(jù)插值,3種方法插值得到的三維圖如圖3所示。

圖3 干支流交匯處地形插值三維結(jié)果

從圖3可以看出,反距離權(quán)重法和克里金法在河口位置的插值結(jié)果沒有很大區(qū)別,在一定程度上反映出地形明顯起伏情況；但三角網(wǎng)格法插值結(jié)果平坦化,忽略了地形起伏較劇烈的部分,未能很好地反映出吳淞口的地形變化,使地形插值結(jié)果產(chǎn)生部分失真。

2.3 復(fù)式分汊河道

對(duì)于復(fù)式分汊河道,選取黃浦江上游一段如圖4所示的分汊河道進(jìn)行對(duì)比分析,中間部分為凸起的小島,將該河道分割為兩個(gè)支汊。圖4所示斷面為實(shí)測(cè)斷面位置,是根據(jù)復(fù)式河道測(cè)量的規(guī)范要求布設(shè)的,斷面間距平均0.4 km左右。采用3種方法對(duì)分汊河道進(jìn)行插值得到的等值線如圖5所示。

圖4 黃浦江復(fù)式分汊河道示意圖

圖5 分汊河道等值線插值結(jié)果

由圖5可知,在斷面滿足規(guī)范要求的情況下,將小島邊界作為固定邊界,相鄰汊道的水下地形插值過(guò)程互不影響,3種方法均能實(shí)現(xiàn)對(duì)分汊河道的地形插值,但三角網(wǎng)格法仍忽略了部分地形實(shí)際起伏情形,反距離權(quán)重法仍然出現(xiàn)了“牛眼”現(xiàn)象,相對(duì)于三角網(wǎng)格法和反距離權(quán)重法,克里金法得到了較好的插值結(jié)果。可見,在斷面合理布設(shè)條件下,3種方法均可以應(yīng)用到復(fù)式河道中進(jìn)行插值,且相較于其他兩種方法而言,克里金法對(duì)復(fù)式河道仍然表現(xiàn)出較好的適應(yīng)性及強(qiáng)健性。

3 斷面減少情況下的插值結(jié)果比較

為了考量克里金法插值效果的強(qiáng)健性,開展了減少斷面數(shù)量測(cè)試檢驗(yàn),并對(duì)其精度進(jìn)行了定量分析。

3.1 干支流交匯處

選取吳淞口處的兩個(gè)間隔斷面,在斷面量減少的情況下,用同樣方式進(jìn)行插值,與原實(shí)測(cè)斷面數(shù)據(jù)進(jìn)行對(duì)比,并對(duì)其插值結(jié)果進(jìn)行誤差分析,結(jié)果見表1。

表1 干支流交匯處插值結(jié)果誤差分析

根據(jù)插值結(jié)果,3種插值方法得到的地形數(shù)據(jù)與實(shí)測(cè)值有一定差異,三角網(wǎng)格法的地形數(shù)據(jù)變化較劇烈,與實(shí)測(cè)值相差較大,反距離權(quán)重法和克里金法的插值結(jié)果變化梯度則較小。從誤差分析結(jié)果來(lái)看,三角形網(wǎng)格法誤差最大,反距離權(quán)重法次之,克里金法誤差最小。實(shí)測(cè)值與計(jì)算值的方差結(jié)果表明,三角網(wǎng)格法方差最大,克里金法最小。

3.2 復(fù)式分汊河道

用相同方法對(duì)復(fù)式分汊河道進(jìn)行插值及誤差分析,結(jié)果見表2。

表2 復(fù)式分汊河道插值結(jié)果誤差分析

根據(jù)插值結(jié)果,克里金法在原斷面處的插值結(jié)果與實(shí)測(cè)結(jié)果比較接近,表明在斷面減少情況下,該方法仍能得到理想的插值結(jié)果,與不減少斷面情況下插值效果最為接近。同時(shí)由表2可知,復(fù)式分汊河道的誤差分析結(jié)果仍是三角形網(wǎng)格法誤差最大,反距離權(quán)重法次之,克里金法誤差最小,表明克里金法具有一定的穩(wěn)定性。

克里金法通過(guò)設(shè)計(jì)變差函數(shù)來(lái)實(shí)現(xiàn)局部加權(quán)插值,克服了其他方法對(duì)起伏地形插值不穩(wěn)定的問(wèn)題,克里金法能夠保持實(shí)測(cè)值與估計(jì)值的偏差期望為零,其結(jié)果是一種最優(yōu)無(wú)偏估計(jì)量,在實(shí)際應(yīng)用中保持插值誤差小是其一個(gè)顯著特點(diǎn),因此在斷面減少的情況下仍能保持較強(qiáng)的健壯性和較高的精度。

4 黃浦江潮位模擬結(jié)果比較

用3種方法對(duì)整個(gè)黃浦江干流地形進(jìn)行插值,用得到的地形數(shù)據(jù)采用EFDC模型對(duì)黃浦江潮位進(jìn)行模擬,計(jì)算時(shí)段為2014年7月1—5日,除地形外,在其他計(jì)算條件相同情況下,潮位模擬結(jié)果如圖6所示。

圖6 潮位模擬結(jié)果對(duì)比

由圖6可知,以2014年7月1日的首個(gè)最高潮位誤差來(lái)看,實(shí)測(cè)潮位值為3.80m,三角網(wǎng)格法地形插值結(jié)果得到的模擬潮位值為3.62m,誤差18cm,反距離權(quán)重法地形插值結(jié)果得到的模擬潮位值為3.73m,誤差7cm,克里金法插值得到的模擬潮位值為3.83m,誤差3cm,后續(xù)計(jì)算時(shí)段高潮位誤差與首個(gè)高潮位誤差接近。在其他條件保持不變的情況下,克里金法插值地形數(shù)據(jù)得到的潮位模擬結(jié)果誤差最小,反距離權(quán)重法次之,三角網(wǎng)格法誤差最大。由此可知,不同地形插值方法對(duì)數(shù)值模擬的計(jì)算精度有顯著影響,插值結(jié)果越接近真實(shí)地形,其模擬精度越高,克里金法對(duì)黃浦江地形插值效果最好,反距離權(quán)重法次之,三角形網(wǎng)格法最差。

5 結(jié) 論

a.對(duì)河道強(qiáng)彎段、干支流交匯處及復(fù)式分汊河道區(qū)域,克里金法較好地避免了反距離權(quán)重法出現(xiàn)的“牛眼”現(xiàn)象,且比三角網(wǎng)格法能更好表現(xiàn)地形實(shí)際起伏變化。

b.克里金法通過(guò)設(shè)計(jì)變差函數(shù)來(lái)實(shí)現(xiàn)局部加權(quán)插值,其結(jié)果是一種最優(yōu)無(wú)偏估計(jì)量,在斷面數(shù)量減少的條件下仍能表現(xiàn)出其較強(qiáng)的健壯性及較高的插值精度。

c.其他計(jì)算條件相同的情況下,3種方法對(duì)黃浦江全段地形插值并進(jìn)行潮位模擬計(jì)算的精度仍是克里金法最高。

參考文獻(xiàn):

[1]何磊,袁德奎,陶建華.一種適用波浪潮流計(jì)算的地形插值算法[J].港工技術(shù),2004（4）:4-6.（HE Lei,YUAN Dekui, TAO Jianhua.An interpolationmethod forgenerating topography for computation of wave and tide [J].Port Engineering Technology,2004（4）:4-6.（in Chinese））

[2]陳中,夏云峰,聞云呈.一種適用于河道水流數(shù)值模擬的地形插值方法[C] //第十四屆中國(guó)海洋（岸）工程學(xué)術(shù)討論會(huì)論文集:下冊(cè).北京:海洋出版社,2009:1230-1233.

[3]劉兆平,楊進(jìn),武煒.地球物理數(shù)據(jù)網(wǎng)格化方法的選取[J].物探與化探,2010,34（1）:93-97.（LIU Zhaoping, YANG Jin,WU Wei.The choice ofgriddingmethods forgeophysical data [J].Geophysical andgeochemical Exploration,2010,34（1）:93-97.（in Chinese））

[4]牛文杰.基于變異函數(shù)的距離加權(quán)反比法[J].東華大學(xué)學(xué)報(bào)（自然科學(xué)版）,2011,37（3）:362-367.（NIU Wenjie.Inverse distance weightingmethod based on variogram [J].Journal of Donghua University（Natural Science）,2011,37（3）:362-367.（in Chinese））

[5] DETWEILER Z R,FERRIS J B.Interpolationmethods for high-fidelity three-dimensional terrain surfaces [J].Journal of Terramechanics, 2010,47（4）:209-217.

[6]李曉軍,王長(zhǎng)虹,朱合華.Kriging插值方法在地層模型生成中的應(yīng)用[J].巖土力學(xué),2009（1）:157-162.（LI Xiaojun, WANG Changhong, ZHU Hehua.Kriging interpolation and its application togenerating stratummodel[J].Rock and Soilmechanics,2009（1）:157-162.（in Chinese））

[7]張殿新,袁德奎,陶建華.波浪、潮流數(shù)值模擬中的網(wǎng)格數(shù)字地形的生成及可視化[J].港工技術(shù),2005（4）:5-8.（ZHANG Dianxin,YUAN Dekui,TAO Jianhua.Digital terrain is produced and visualized in wave and tide[J].Port Engineering Technology, 2005（4）: 5-8.（in Chinese））

[8]包世泰,廖衍旋,胡月明,等.基于Kriging的地形高程插值[J].地理與地理信息科學(xué),2007（3）:28-32.（BAO Shitai, LIAO Yanxuan, HU Yueming, et al.Terrain interpolation based on Krigingmethod[J].Geography andgeo-Information Science,2007（3）:28-32.（in Chinese））

[9]杜國(guó)明,汪光松,吳超羽,等.克里金在珠江河道地形空間數(shù)據(jù)內(nèi)插中的應(yīng)用[J].中山大學(xué)學(xué)報(bào)（自然科學(xué)版）, 2007, 46（1）: 119-122.（DUguoming, WANGguangsong, WU Chaoyu, et al.Application of Kriging interpolation to spatial data of Pearl River riverway[J].Acta Scientiarum Naturalium Universitatis Sunyatseni, 2007,46（1）:119-122.（in Chinese））

[10]曾懷恩,黃聲享.基于Kriging方法的空間數(shù)據(jù)插值研究[J].測(cè)繪工程,2007,16（5）:5-8.（ZENG Huaien, HUANG Shengxiang.Research on spatial data interpolation based on Kriging interpolation [J].Engineering of Surveying andmapping,2007,16（5）:5-8.（in Chinese））

[11]蔡迪花,郭鈮,李崇偉.基于DEM的氣溫插值方法研究[J].干旱氣象,2009,27（1）:10-17.（CAI Dihua,GUO Ni,LI Chongwei.Interpolation of air temperature based on DEM over eastern region ofgansu[J].Journal of Aridmeteorology,2009,27（1）:10-17.（in Chinese））

[12]李新,程國(guó)棟,盧玲.空間內(nèi)插方法比較[J].地球科學(xué)進(jìn)展,2000,15（3）:260-265.（LI Xin,CHENGguodong, LU Ling.Comparison of spatial interpolationmethods[J].Advances in Earth Science,2000,15（3）:260-265.（in Chinese））

中圖分類號(hào):P217

文獻(xiàn)標(biāo)志碼:A

文章編號(hào):1006- 7647（2016）03- 0016- 04

DOI:10.3880/j.issn.1006- 7647.2016.03.004

基金項(xiàng)目:國(guó)家自然科學(xué)基金（51379060,51179052）；上海市水務(wù)局科研項(xiàng)目（2014-10）

作者簡(jiǎn)介:華祖林（1965—）,男,教授,博士,主要從事水動(dòng)力與水環(huán)境模擬研究。E-mail:zulinhua@ hhu.edu.cn

收稿日期:（2015- 10 16 編輯:熊水斌）

Comparison of differentmethods for interpolation of topography of discrete rivers

HUA Zulin1,2,3, WANGHaiyan1,2, WANG Liang1,2, WANG Yulin1,2（1.Key Laboratory of Integrated Regulation and Resource Development on Shallow Lakes,ministry of Education, Hohai University, Nanjing 210098, China；2.College of Environment, Hohai University, Nanjing 210098, China；3.National Engineering Research Center of Water Resources Efficient Utilization and Engineering Safety, Hohai University, Nanjing 210098, China）

Abstract:According to the actual situation in which the cross-sectional topographic data are discontinuous for hydrodynamic and aquatic environmental simulations, topographic interpolations for strongly curved reaches, river confluence, and braided rivers were conducted with the triangularmeshmethod, inverse distance weightingmethod, and Krigingmethod using the Huangpu River as an example.The reliability of these three interpolationmethods was analyzed and compared.The results show that, for strongly curved reaches, river confluence, and braided rivers where topography varies significantly, the interpolation results of the triangularmeshmethod are unstable, with high errors and significant fluctuations, indicating that thismethod failed to reflect the actual topographic situation.The“buphthalmos”phenomenonmay appear when the inverse distance weightingmethod is used.Compared with othermethods, the Krigingmethod shows higher reliability.When fewer cross-sectional data are used, the results of the Krigingmethod show lower error and higher robustness than the other twomethods.Based on the interpolation results of the topography of the Huangpu River with the above-mentioned three interpolationmethods, amathematicalmodel was used to compare the tidal level processes, and the results show that the computational precision of the Krigingmethod is the highest.

Key words:discrete river；topographic interpolation；triangularmeshmethod；inverse distance weightingmethod；Krigingmethod；Huangpu River