河道二維數(shù)值模擬地形數(shù)據(jù)前處理精細(xì)化研究

趙科鋒，曹慧群，羅平安

(長(zhǎng)江科學(xué)院 流域水環(huán)境研究所，武漢 430010)

1 研究背景

隨著經(jīng)濟(jì)及城鎮(zhèn)化發(fā)展，地表水環(huán)境問(wèn)題越來(lái)越突出，尤其是城鎮(zhèn)河流、湖泊等地表水體受到越來(lái)越嚴(yán)重的污染。為了研究地表水水環(huán)境變化狀況，許多學(xué)者通過(guò)建立數(shù)值模型，分析地表水水動(dòng)力水質(zhì)變化規(guī)律。國(guó)內(nèi)外針對(duì)地表水水動(dòng)力、水質(zhì)、泥沙的數(shù)值模擬軟件主要有EFDC、MIKE、Delft3D、CJK3D等[1]，而地表水?dāng)?shù)值模擬的一個(gè)重要內(nèi)容是地形模擬，地形模擬的好壞直接影響模擬結(jié)果的準(zhǔn)確性與可信度[2]。針對(duì)地形模擬，部分學(xué)者開(kāi)展了相關(guān)研究，如：任勇等[3]對(duì)金沙江河道不規(guī)則三角網(wǎng)進(jìn)行人工干預(yù)和優(yōu)化，達(dá)到快速、合理、準(zhǔn)確生成等高線的目的。侯明瑞[4]將元胞方法引入地形重構(gòu)過(guò)程中，通過(guò)將復(fù)雜的河流動(dòng)力學(xué)過(guò)程進(jìn)行簡(jiǎn)化，建立基于元胞方法的地形重構(gòu)模型，此種模型能較好避免地形構(gòu)造生成的非自然條帶狀地形。陳中[2]利用SMS軟件，提出了一種基于三角網(wǎng)格的內(nèi)插技術(shù)。而在眾多的研究當(dāng)中，針對(duì)地形數(shù)據(jù)優(yōu)化的方法卻較少。地形模擬的優(yōu)劣極大影響模型的準(zhǔn)確性，而受人力、物力、財(cái)力及技術(shù)等條件的限制，無(wú)法對(duì)河流、湖泊等地表水體進(jìn)行全覆蓋測(cè)量，實(shí)際操作中，僅能進(jìn)行有限的斷面測(cè)量，導(dǎo)致地形模擬嚴(yán)重不符合實(shí)際情況，影響模型準(zhǔn)確性，甚至造成模型不收斂，導(dǎo)致無(wú)法計(jì)算。

MIKE是由丹麥的DHI公司開(kāi)發(fā)的一款非常優(yōu)秀的應(yīng)用于水動(dòng)力、水質(zhì)、泥沙輸移等方面的軟件，MIKE可用于河流、湖泊、水庫(kù)、河口、海灣等的一維、二維、三維模擬[5]，水底地形插值是MIKE建模必不可少的重要環(huán)節(jié)[6]。目前，研究者使用較多的斷面地形細(xì)化方法，是通過(guò)第三方模擬軟件進(jìn)行插值后，導(dǎo)出插值數(shù)據(jù)，而這種方法操作復(fù)雜，依然無(wú)法解決任意斷面數(shù)及任意距離的插值計(jì)算，且常常導(dǎo)致地形數(shù)據(jù)出現(xiàn)較大的偏差[7-8]。為了解決因?qū)崪y(cè)斷面地形數(shù)據(jù)不足導(dǎo)致的模型插值與實(shí)際情況偏差較大的問(wèn)題，本文通過(guò)對(duì)研究區(qū)河道邊界節(jié)點(diǎn)編碼，利用實(shí)測(cè)斷面坐標(biāo)對(duì)編碼進(jìn)行識(shí)別，可沿河道彎曲方向插值，實(shí)現(xiàn)縱向任意斷面數(shù)、橫向任意距離的河道地形插值優(yōu)化，避免因地形數(shù)據(jù)量少而生成異常地形的問(wèn)題，達(dá)到減少模型運(yùn)算錯(cuò)誤率及模型更接近真實(shí)情況的目的。

2 地形數(shù)據(jù)前處理精細(xì)化方法

2.1 邊界及斷面數(shù)據(jù)編碼

MIKE21模型邊界基礎(chǔ)數(shù)據(jù)是由連續(xù)坐標(biāo)點(diǎn)組成的，描述連續(xù)坐標(biāo)點(diǎn)的方法應(yīng)用較多的是樣條曲線法[9]，然而對(duì)于彎曲較大的河道，樣條曲線無(wú)法準(zhǔn)確描述河道邊界形狀。為此，本文提出對(duì)河道邊界節(jié)點(diǎn)進(jìn)行編碼(編碼與坐標(biāo)一一對(duì)應(yīng))，通過(guò)遍歷所有邊界節(jié)點(diǎn)，識(shí)別出距實(shí)測(cè)斷面最近的邊界節(jié)點(diǎn)，根據(jù)定位的邊界節(jié)點(diǎn)及設(shè)置的插值斷面距離(如圖1中x)，計(jì)算內(nèi)插斷面的編碼(如圖1中A1、A2等)，由此可識(shí)別出所有插值斷面的邊界節(jié)點(diǎn)編碼，進(jìn)而確定新增斷面坐標(biāo)。節(jié)點(diǎn)編碼及識(shí)別方法為：利用ArcGIS對(duì)河道左、右岸邊界進(jìn)行等間距取節(jié)點(diǎn)(間距設(shè)置成1 m為宜)，并對(duì)每個(gè)節(jié)點(diǎn)進(jìn)行編碼，通過(guò)節(jié)點(diǎn)遍歷，識(shí)別出距實(shí)測(cè)斷面最近的節(jié)點(diǎn)。需要注意，邊界文件(MIKE21能識(shí)別擴(kuò)展名為xyz的格式)分左、右岸2個(gè)節(jié)點(diǎn)坐標(biāo)文件，斷面上實(shí)測(cè)點(diǎn)需按上游至下游、左岸至右岸的規(guī)則編寫文件。河道左、右岸分別為ABC、abc，實(shí)測(cè)斷面為Aa、Bb、Cc，新增內(nèi)插橫向斷面為A1a1、A2a2等，新增內(nèi)插縱向斷面為Ⅰ、Ⅱ、Ⅲ，如圖1所示。

圖1 河道邊界編碼計(jì)算示意圖Fig.1 Schematic diagram of coding of river boundary

2.2 內(nèi)插節(jié)點(diǎn)坐標(biāo)及高程計(jì)算

內(nèi)插節(jié)點(diǎn)坐標(biāo)計(jì)算包括橫向斷面及縱向斷面節(jié)點(diǎn)坐標(biāo)計(jì)算。橫斷面節(jié)點(diǎn)坐標(biāo)計(jì)算：根據(jù)河道左、右岸邊界上確定的斷面節(jié)點(diǎn)及設(shè)置的插值斷面間距，計(jì)算內(nèi)插斷面的編碼，由編碼確定相應(yīng)的坐標(biāo)。縱斷面節(jié)點(diǎn)坐標(biāo)計(jì)算：根據(jù)每條橫斷面兩端節(jié)點(diǎn)，確定等間距縱向斷面節(jié)點(diǎn)，擴(kuò)展方式如圖2所示。

圖2 節(jié)點(diǎn)坐標(biāo)內(nèi)插示意圖Fig.2 Schematic diagram of node coordinate interpolation

主要步驟有：①?gòu)牡貓D上矢量化河道邊界(河道左岸及右岸)，得到KML或SHP等格式的文件；②利用ArcGIS對(duì)河道左、右岸分別進(jìn)行等間距取點(diǎn)，得到左、右岸多組邊界節(jié)點(diǎn)坐標(biāo)文件；③斷面上實(shí)測(cè)點(diǎn)按上游至下游、左岸至右岸的規(guī)則編寫文件；④利用實(shí)測(cè)河道橫斷面距河道最近一個(gè)節(jié)點(diǎn)坐標(biāo)，遍歷對(duì)比河道邊界所有節(jié)點(diǎn)坐標(biāo)，獲取實(shí)測(cè)斷面與河岸邊界最近節(jié)點(diǎn)(此節(jié)點(diǎn)與交叉點(diǎn)相距極小，近似為交叉點(diǎn)，如圖1中點(diǎn)A、B、C)的編碼；⑤根據(jù)設(shè)置的橫向斷面剖分間距x(如圖1中所示)，計(jì)算兩實(shí)測(cè)斷面間河道邊界等間距內(nèi)插節(jié)點(diǎn)編碼，如圖1中點(diǎn)A1、B1等；⑥根據(jù)縱向剖分?jǐn)嗝鏀?shù)量，計(jì)算實(shí)測(cè)斷面(圖1中Aa、Bb等)及內(nèi)插斷面(A1a1、B1b1等)相應(yīng)的河道內(nèi)等間距節(jié)點(diǎn)編碼；⑦根據(jù)2條實(shí)測(cè)斷面河底高程，插值計(jì)算兩實(shí)測(cè)斷面間河道內(nèi)節(jié)點(diǎn)河底高程值及相應(yīng)編碼；⑧對(duì)各剖分節(jié)點(diǎn)坐標(biāo)與河底高程進(jìn)行編碼匹配，最終得到多組河底精細(xì)化坐標(biāo)點(diǎn)及對(duì)應(yīng)的河底高程。

2.3 插值計(jì)算關(guān)鍵代碼模塊

由于計(jì)算量較大，手動(dòng)計(jì)算費(fèi)時(shí)費(fèi)力，且出錯(cuò)率高，通過(guò)編寫程序，可便捷實(shí)現(xiàn)河道地形精細(xì)化內(nèi)插，獲取大量地形數(shù)據(jù)，達(dá)到快速、準(zhǔn)確建立數(shù)值模型的目的。關(guān)鍵代碼模塊主要包括：①河道左、右岸邊界距斷面最近點(diǎn)求解模塊；②實(shí)測(cè)橫向斷面間內(nèi)插斷面編碼計(jì)算模塊；③縱向斷面內(nèi)插節(jié)點(diǎn)編碼計(jì)算模塊；④內(nèi)插節(jié)點(diǎn)河底高程插值模塊；⑤內(nèi)插節(jié)點(diǎn)坐標(biāo)與高程匹配模塊；⑥數(shù)據(jù)顯示、輸出模塊。

3 實(shí)例應(yīng)用

為了檢驗(yàn)河道地形內(nèi)插生成精細(xì)化地形數(shù)據(jù)的準(zhǔn)確性與適用性，本文選取荊門城區(qū)某河流部分河段作為研究對(duì)象，根據(jù)河道實(shí)測(cè)斷面地形數(shù)據(jù)，沿河道彎曲方向進(jìn)行精細(xì)化插值計(jì)算，并選取部分?jǐn)嗝鎸?shí)測(cè)數(shù)據(jù)進(jìn)行對(duì)比驗(yàn)證。

3.1 河道基本情況

研究河道長(zhǎng)約1.0 km，寬約18～60 m，河道彎曲變化較大。為研究河岸污水處理廠排污對(duì)下游監(jiān)測(cè)斷面水質(zhì)的影響，擬建立二維水動(dòng)力水質(zhì)模型進(jìn)行預(yù)測(cè)，為此，開(kāi)展了河道地形測(cè)量，現(xiàn)場(chǎng)共測(cè)量斷面10組，每個(gè)斷面等間距測(cè)量7組高程點(diǎn)，共計(jì)70個(gè)測(cè)量點(diǎn)，斷面位置如圖3所示。

圖3 河道地形測(cè)量斷面位置Fig.3 Location of sections for river topographic survey

3.2 河道地形精細(xì)化插值及檢驗(yàn)

利用6個(gè)實(shí)測(cè)斷面地形數(shù)據(jù)進(jìn)行地形精細(xì)化內(nèi)插計(jì)算，4個(gè)實(shí)測(cè)斷面地形數(shù)據(jù)作為檢驗(yàn)。橫向斷面間距設(shè)置為20 m，縱向斷面數(shù)設(shè)置為7，利用Python編寫的程序自動(dòng)生成約60組斷面，共約420個(gè)坐標(biāo)點(diǎn)，如圖4所示。

圖4 經(jīng)程序插值后生成內(nèi)插節(jié)點(diǎn)Fig.4 Nodes generated by programming afterinterpolation

根據(jù)生成的內(nèi)插節(jié)點(diǎn)坐標(biāo)，利用MIKE鄰近點(diǎn)法插值，生成研究區(qū)河道地形，如圖5(a)所示。未經(jīng)斷面內(nèi)插計(jì)算，直接通過(guò)10組實(shí)測(cè)地形數(shù)據(jù)生成的河道地形如圖5(b)所示。從圖5(a)、圖5(b)的對(duì)比可以明顯看出，經(jīng)河道斷面精細(xì)化內(nèi)插計(jì)算生成的河道地形較為合理，河道地形較為平滑；而未經(jīng)河道斷面精細(xì)化插值計(jì)算，直接利用實(shí)測(cè)斷面生成的河道地形明顯不準(zhǔn)確，甚至出現(xiàn)河道中泓地形高于河道兩岸地形的情況，這主要是由于河道較長(zhǎng)，河道寬度較小，而實(shí)測(cè)斷面相距較遠(yuǎn)，數(shù)據(jù)量較小，導(dǎo)致MIKE軟件插值不合理。

圖5 河道地形精細(xì)化插值和未精細(xì)化插值計(jì)算地形Fig.5 Topography maps before and after refined riverchannel interpolation

本文利用4組實(shí)測(cè)斷面數(shù)據(jù)進(jìn)一步檢驗(yàn)精細(xì)化插值計(jì)算后生成的河道地形的準(zhǔn)確性。通過(guò)對(duì)比河道斷面對(duì)應(yīng)高程實(shí)測(cè)值與精細(xì)化計(jì)算值及未精細(xì)化計(jì)算值(通過(guò)第三方軟件插值計(jì)算)，分析河道地形精細(xì)化內(nèi)插后的準(zhǔn)確性，如圖6所示。28組地形高程實(shí)測(cè)值與精細(xì)化計(jì)算值中，相差最大為0.148 5 m，最小為0.012 1 m，平均為0.076 6 m，相對(duì)誤差均不到1%；而實(shí)測(cè)值與未精細(xì)化計(jì)算值中，相差最大為1.510 5 m，最小為0.011 0 m，平均為0.511 4 m。從分析結(jié)果來(lái)看，經(jīng)程序精細(xì)化內(nèi)插計(jì)算后的河道地形能較好地反映原始河道地形，為研究區(qū)建立二維水動(dòng)力水質(zhì)模型提供了大量可靠的基礎(chǔ)數(shù)據(jù)，同時(shí)，減少了因數(shù)據(jù)量少導(dǎo)致的模擬地形不合理的情況，提高了數(shù)值建模的效率及準(zhǔn)確性。

圖6 斷面地形高程實(shí)測(cè)值與計(jì)算值對(duì)比Fig.6 Comparison of terrain elevation between measuredand calculated values

4 結(jié) 論

(1)地形模擬的優(yōu)劣極大影響了模型的準(zhǔn)確性，而受人力、物力、財(cái)力及技術(shù)等條件的限制，無(wú)法對(duì)河流、湖泊等地表水體進(jìn)行全覆蓋地形測(cè)量，實(shí)際多數(shù)操作中，僅能進(jìn)行有限的斷面測(cè)量，導(dǎo)致地形模擬嚴(yán)重不符合實(shí)際情況，甚至導(dǎo)致模型不收斂。

(2)樣條曲線法無(wú)法準(zhǔn)確描述彎曲較大的河道，而通過(guò)對(duì)河道連續(xù)點(diǎn)編碼定位坐標(biāo)點(diǎn)，結(jié)合設(shè)置的插值橫向斷面間距及插值縱向斷面數(shù)，可以沿河道方向計(jì)算內(nèi)插節(jié)點(diǎn)坐標(biāo)及高程，此方法能較好地精細(xì)化河道地形。

(3)選取荊門城區(qū)某河流部分河段進(jìn)行檢驗(yàn)，利用編寫的程序結(jié)合實(shí)測(cè)斷面精細(xì)化河道地形，對(duì)比地形實(shí)測(cè)值與計(jì)算值，可以看出，計(jì)算河道斷面與實(shí)測(cè)斷面基本一致，測(cè)量點(diǎn)高程相差最大為0.148 5 m，平均為0.076 6 m，精細(xì)化內(nèi)插生成的地形基本能反映原始河道地形。

(4)通過(guò)編寫程序，能實(shí)現(xiàn)任意橫向斷面間距及任意縱向斷面數(shù)的精細(xì)化坐標(biāo)計(jì)算，此方法生成的河道地形數(shù)據(jù)，能應(yīng)用于MIKE、EFDC、SMS、Delf3D等地表水?dāng)?shù)值模擬軟件，并減少模型出錯(cuò)率，提高模型精度。