洪明海,嚴(yán) 濤,劉 輝,王興建,敬 娜,汪仕偉,魏文杰,李析男,馮楚橋

(貴州省水利水電勘測設(shè)計研究院有限公司,貴州 貴陽 550002)

為全面推進(jìn)算據(jù)、算法、算力建設(shè),加快構(gòu)建具有“四預(yù)”(預(yù)報、預(yù)警、預(yù)演、預(yù)案)功能的智慧水利體系,水利部發(fā)布的《“十四五”智慧水利建設(shè)實施方案》《數(shù)字孿生流域建設(shè)技術(shù)大綱(試行)》《數(shù)字孿生水利工程建設(shè)技術(shù)導(dǎo)則(試行)》《水利業(yè)務(wù)“四預(yù)”功能基本技術(shù)要求(試行)》等指導(dǎo)文件和技術(shù)要求[1],明確提出數(shù)據(jù)底板是數(shù)字孿生流域建設(shè)的重點任務(wù)之一,數(shù)字高程模型(DEM)則是數(shù)據(jù)底板中最重要的數(shù)據(jù)之一。

DEM于1958年由國外提出,基于DEM的數(shù)字地形分析在遙感測繪[2]、資源調(diào)查、環(huán)境保護(hù)、工程建設(shè)、災(zāi)害防治、水利規(guī)劃設(shè)計[3]以及地學(xué)研究等各方面都展現(xiàn)出越來越重要的作用。通常,DEM主要來源于攝影測量、已有地形圖數(shù)字化、已有的DEM庫中提取、地面測量等。通過地面測量的等高線生成DEM不僅精度可控,操作方便,而且現(xiàn)有技術(shù)相對成熟穩(wěn)定,是工程實踐中最常用的DEM獲取方式之一。然而,當(dāng)研究區(qū)域較大(如一個地級市或者一個省)且等高線為大比例尺(萬分之一或者五千分之一甚至更大)時,常常會因為電腦等硬件設(shè)施的限制使得DEM的獲取不能一次性完成,需要分幅生成后再鑲嵌,這樣在分幅之間的縫隙往往就會產(chǎn)生誤差甚至是Nodata值,給工程應(yīng)用帶來不小的困擾。

目前對于DEM修復(fù)的方法主要有3種:一是利用局部修正后進(jìn)行反距離加權(quán)融合填補[4];二是建立誤差分布模型來判斷系統(tǒng)誤差等指標(biāo)后采用自適應(yīng)平滑[5]、智能濾波[6]等對DEM噪音進(jìn)行修復(fù);三是以別的DEM來做參考校正缺失的區(qū)域[7]。事實上,DEM數(shù)據(jù)集的誤差分布與土地利用類型、坡度、坡向、高程地形、地表植被類型等有關(guān)[8],也就是缺失的像元高程值與周邊的像元是有關(guān)系的。當(dāng)處理大面積大批量的DEM修復(fù)時,高效快速且精度可控的方法顯得尤為重要[9]。

本文擬采用ESRI公司的ArcGIS軟件,通過萬分之一等高線來獲取某市DEM,同時通過Python對生成的DEM分幅處進(jìn)行適當(dāng)?shù)男拚?將修正的結(jié)果和已有的數(shù)據(jù)進(jìn)行比較分析,為大范圍大尺度的等高線生成DEM提供一種借鑒方法。

1 材料與方法

1.1 研究區(qū)概況和數(shù)據(jù)來源

本次研究的區(qū)域位于貴州省某市某局部區(qū)域,研究區(qū)面積約為57 km2,海拔高度為935～1 345 m,境內(nèi)地形起伏大,地貌類型較為復(fù)雜。

本文使用的研究數(shù)據(jù)為測量部門提供的1∶10 000的等高線和高程點數(shù)據(jù),為shp格式,坐標(biāo)系為CGCS2000投影坐標(biāo),其中等高線共4 716段,高程點共415個;對比的DEM為自然資源部門提供的同一區(qū)域數(shù)據(jù),GDB格式,CGCS2000投影坐標(biāo)。

1.2 計算原理

等高線生產(chǎn)DEM[10],主要涉及等高線/高程點轉(zhuǎn)TIN,然后由TIN生成DEM,對DEM進(jìn)行鑲嵌后再修正。

不管是地形轉(zhuǎn)TIN還是TIN轉(zhuǎn)DEM工具,都是一個占用大量內(nèi)存的應(yīng)用程序,因此不能創(chuàng)建較大的輸出柵格。按照ArcGIS官方(https://resources.arcgis.com/zh-cn/help/main/)給出的處理方案“創(chuàng)建多個DEM后,最好使用鑲嵌地理處理工具的BLEND選項或MEAN選項將它們合并”。因此在處理全市乃至全省上萬副大比例尺的等高線生成DEM后最終需要將它們拼接起來[11-12],然后對分幅的縫隙采用MEAN方法修正。

從概念上講,該修正方法會將要修正的像元或中心像元周圍一個3×3的像元鄰域的高程值參與計算。圖1所示,如果鄰域內(nèi)某個像元位置的高程值為NoData,則將中心像元周邊8個有值像元值的平均值指定給該位置,計算見式(1):

a)像元位置

(1)

式中Dx——待修正的DEM像元;Di——周邊8個像元有值的柵格像元值;n——周邊8個像元中有值的個數(shù),最大取8。

例如圖1a中NoData的Dx,當(dāng)其周邊的8個像元值都不為NoData時,

Dx=(D1+D2+D3+D4+D5+D6+D7+D8)/8

(2)

當(dāng)NoData的Dx周邊的8個像元值均為NoData時,則該像元值賦值也為NoData。

當(dāng)修正像元周邊未讀取到3×3的像元鄰域或者讀取到的3×3的像元鄰域中的8個像元中可能會有多個NoData時,則參與計算的周邊像元數(shù)是不為NoData的個數(shù),見圖2。圖2a、2b、2c為周邊鄰域為2×2的像元,且3個像元值分別有1個、2個、3個NoData的像元,當(dāng)2×2的鄰域像元中有3個NoData時,則待修正的像元柵格也為NoData,見圖2c。此外,當(dāng)處理DEM邊界時,還可能讀取到2×3的像元或3×2的像元,見圖2d、2e;其處理方法同上?；蛘咦x取到3×3的像元鄰域中有多個NoData(圖2f),則參與計算的為有值部分。

a)1個像元無值

1.3 研究技術(shù)路線及步驟

本文的研究技術(shù)路線見圖3,其主要操作步驟如下。

圖3 研究技術(shù)路線

步驟一對研究區(qū)范圍內(nèi)的等高線和高程點數(shù)據(jù)處理。將CAD格式的地形圖分別轉(zhuǎn)為等高線shp和高程點shp,并定義正確的投影坐標(biāo)系后,對等高線和高程點數(shù)據(jù)中高程值為0、無高程值、高程值異常的數(shù)據(jù)進(jìn)行細(xì)部識別并處理(刪除或重新賦值)為正常值,然后去除等高線或高程點重合的要素。

步驟二將處理好的等高線和高程點數(shù)據(jù)作為創(chuàng)建TIN的輸入數(shù)據(jù),其中等高線要素作為contour,高程點要素作為masspoint。

步驟三使用前文生成的TIN轉(zhuǎn)DEM,根據(jù)輸入地形圖的比例尺設(shè)置合適的容差和像元大小。

步驟四重復(fù)步驟一、二完成所有單幅地形圖轉(zhuǎn)為DEM后將所有DEM進(jìn)行鑲嵌。

步驟五對鑲嵌后DEM進(jìn)行NoData值修正處理。

步驟六對修正后的DEM進(jìn)行精度評定,符合要求則導(dǎo)出DEM,不符合則回到地形圖的處理重新設(shè)置像元大小、容差等。

1.4 分析指標(biāo)及評定標(biāo)準(zhǔn)

常用的高程評定[13]方法有檢查點法[14-15]、剖面法、分形法、影像分析法、等高線回放法[16-17]等,本文采用中國國家標(biāo)準(zhǔn)化管理委員會發(fā)布的GB/T 18316—2001《數(shù)字測繪產(chǎn)品檢查驗收規(guī)定和質(zhì)量評定》[14]中的檢查點法進(jìn)行分析,即對生成的DEM采用修正方法修正后NoData像元值和已知的DEM值進(jìn)行比較,計算中誤差。中誤差按照GB 50026—2020《工程測量標(biāo)準(zhǔn)》[18]5.9.6計算,見式(3):

(3)

式中m——數(shù)字高程模型(DEM)高程中誤差;n——檢測點數(shù),按照GB/T 18316—2008《數(shù)字測繪成果質(zhì)量檢查與驗收》[15],要求檢測點分布均勻,位置易于辨認(rèn),不少于50個;Ri——檢測點的修正高程值;Zi——檢測點的原始高程值。

當(dāng)m

參照GB 50026—2020《工程測量標(biāo)準(zhǔn)》條文說明5.1.3中,地形圖的基本等高距,是以等高線的高程中誤差的經(jīng)驗公式驗算,數(shù)字高程模型格網(wǎng)間距的選取及格網(wǎng)點高程中誤差應(yīng)符合表1的規(guī)定。

表1 數(shù)字高程模型格網(wǎng)間距的選取及格網(wǎng)點高程中誤差 單位:m

參照5.9.6條數(shù)字高程模型建立后應(yīng)進(jìn)行檢查,并應(yīng)符合下列條件:對于實測數(shù)據(jù)所建立的數(shù)字高程模型,應(yīng)進(jìn)行外業(yè)實測檢查并統(tǒng)計精度。每個圖幅的檢查點數(shù),不應(yīng)少于20點,檢查點與模型插值點的高程較差不應(yīng)大于本標(biāo)準(zhǔn)第5.1.7條相應(yīng)格網(wǎng)點高程中誤差的2倍。等高線高程中誤差mb的取值,對于常用的設(shè)計坡度,均不能大于基本等高距的1/2;對于較大的設(shè)計坡段,也不能大于基本等高距。

此外,參照DL/T 5001—2014《火力發(fā)電廠工程測量技術(shù)規(guī)程》[19]中第6.1.4條等高線插值點或相對與臨近圖根點的高程中誤差的規(guī)定,見表2。

表2 等高線插值求點的高程中誤差 單位:m

綜上,本研究區(qū)地形類別為山地和高山地,為1∶10 000比例尺,數(shù)字高程模型的插值點高程中誤差限值m0取1倍基本等高距(格網(wǎng)尺寸),為5 m。

2 結(jié)果與討論

2.1 修正方法的適用性

地理學(xué)第一定律指出,任何事物都是與其他事物相關(guān)的,相近的事物關(guān)聯(lián)更緊密[20],空間自相關(guān)分析對具有地理坐標(biāo)的要素進(jìn)行相關(guān)程度表征[21],即空間分布模式(spatial distribution pattern),一般有聚集、離散和隨機3種模式,而空間自相關(guān)中用來表征空間分布模式的量化指數(shù)一般使用莫蘭指數(shù)[22]。DEM是典型的具有坐標(biāo)信息的要素,其高程值常常和地形的變化存在很大的聯(lián)系,即空間相關(guān)性。

圖4是比較值DEM轉(zhuǎn)換為高程像元方格后做的空間自相關(guān)分析結(jié)果,參與計算的78 018個單元值(一般莫蘭指數(shù)分析樣本大于30個即可)表明DEM的高程分布存在空間上的高度聚集(圖4 Clustered)。分析結(jié)果數(shù)據(jù)見表3,其中莫蘭指數(shù)為0.999 796,莫蘭指數(shù)代表數(shù)據(jù)的相關(guān)程度,其值范圍為[-1,1],本次自相關(guān)分析表明高程值和位置關(guān)系相關(guān)性非常大,即高值與高值發(fā)生聚集,低值與低值聚集,空間上呈現(xiàn)正相關(guān)模式;表3中的Z得分和p值則代表分析的置信度,當(dāng)Z>2.58且p<0.01時其置信度為99%,本次分析結(jié)果Z得分為393.212 387 369,則隨機產(chǎn)生此聚類模式的可能性小于1%,因此高程值和地理位置是高度自相關(guān)的,這也印證了ArcGIS官方推薦的鑲嵌DEM柵格時為處理接幅的縫隙采用“BLEND”或者“MEAN”方法是適宜的,也即是表明本文對NoData采用周邊9個有值像元修正的方法是合適的。

表3 空間自相關(guān)分析結(jié)果數(shù)據(jù)

圖4 比較值DEM空間自相關(guān)分析圖莫蘭指數(shù)

2.2 修正結(jié)果

圖5是接幅縫隙DEM的NoData值修正前后的局部截圖比較,其中圖5a為修正前,其接幅處存在空白的NoData值,會影響DEM的后續(xù)數(shù)據(jù)分析,比如集雨面積量算;圖5b是修正后的像元,圖5c是比較值的像元,可以看到采用本文提出的方法修正的像元高程值和比較值的差異非常小,圖5中展示的39個像元值中37個的修正值和比較值相差均在1 m以內(nèi),遠(yuǎn)遠(yuǎn)低于規(guī)范規(guī)定的誤差限5 m和采樣間距5 m,最大值也僅為1.11 m,也滿足規(guī)范要求的誤差限值。

a)NoData像元空值

圖6是本次研究區(qū)域內(nèi)接幅縫隙中306個NoData像元的修正值和比較值的點繪圖,可以看出,306對數(shù)據(jù)幾乎均勻分布在1∶1線附近,表明采用本文提出的方法修正的DEM高程值和真值(比較值)非常接近,數(shù)據(jù)計算結(jié)果也表明修正值和比較值的最大差值僅為3.55 m(圖7),也低于規(guī)范規(guī)定的誤差限5 m和采樣間距5 m。

圖6 修正值與比較值數(shù)據(jù)點比較

圖7 修正值與比較值高程差

2.3 中誤差分析

將研究區(qū)域306個NoData像元的修正值和比較值按照各自所在的空間位置關(guān)系,自然形成11組,分別計算中誤差,結(jié)果見圖8,中誤差不等于真誤差,它僅是一組真誤差的代表值。中誤差的大小反映了該組觀測值精度的高低,因此,通常稱中誤差為觀測值的中誤差。其中,中誤差最大值為1.22,是第三組,最小值是第7組的0.37。此外,還可以看出11組的中誤差值和總體的中誤差值0.75都在0.3倍誤差限的下方,低于規(guī)范要求的中誤差限值。

圖8 中誤差分布

按照GB/T 18316—2008《數(shù)字測繪成果質(zhì)量檢查與驗收》[15],當(dāng)成果質(zhì)量符合要求后還可以計算質(zhì)量元素分值S,并對單位成果質(zhì)量進(jìn)行等級評定。計算見式(4):

(4)

式中S——質(zhì)量元素分值,其分值對應(yīng)的等級分類見表4;m0——允許中誤差值,參照表1取值;m——數(shù)字高程模型(DEM)高程中誤差,按式(3)計算。

表4 單位成果質(zhì)量評定等級

根據(jù)GB/T 18316—2008《數(shù)字測繪成果質(zhì)量檢查與驗收》規(guī)范,采用本方法修正的研究區(qū)域的DEM高程值的中誤差均小于0.3倍中誤差,其質(zhì)量元素分值S為100分,按照表4的等級評定結(jié)果為優(yōu)級品。

3 結(jié)論

本文通過Python調(diào)用ArcGIS相應(yīng)的工具來批量修正某市某研究區(qū)域的地形轉(zhuǎn)DEM后接縫處存在的高程NoData值,其主要結(jié)論如下。

a)方法適用性強。參與計算的78 018個單元值的自相關(guān)分析表明DEM的高程值存在高度的自相關(guān)性,空間分布模式結(jié)果為聚類,其中莫蘭指數(shù)為0.999 796,Z得分為393.212 387 369且p小于0.1,則隨機產(chǎn)生此聚類模式的可能性小于1%。此外ESRI官方推薦的鑲嵌DEM柵格時為處理接幅的縫隙采用“BLEND”或者“MEAN”方法。因此本文提出DEM的NoData值采用周邊8個有值像元的均值來修正的方法是可行的。

b)修正結(jié)果可靠。局部計算39個像元值中37個的修正值和比較值相差均在1 m以內(nèi),最大值也僅為1.11 m;本次研究區(qū)域內(nèi)接幅縫隙中306個NoData像元的修正值和比較值的數(shù)據(jù)幾乎均勻分布在1∶1線附近,修正值和比較值的最大差值僅為3.55 m,低于規(guī)范規(guī)定的誤差限5 m和采樣間距5 m。11組中誤差最大值為1.22,最小值是0.37,且總體的中誤差值0.75都在0.3倍誤差限內(nèi),低于規(guī)范要求的中誤差限值,其質(zhì)量元素分值S為100分,評定結(jié)果為優(yōu)級品。

c)計算速度快,節(jié)約時間成本。從柵格計算的原理出發(fā),采用Python調(diào)用ArcGIS將處理圖像問題轉(zhuǎn)化為矩陣運算問題,而矩陣運算速度極快則是Numpy最大的優(yōu)勢,本研究區(qū)域的修正計算用時僅為10 s(ArcGIS暫時沒有直接的修正工具,采用鑲嵌柵格也不一定能覆蓋所有的接縫處NoData值)。此外,本方法還可以批量處理,通過矩陣分塊快速完成任務(wù),不存在內(nèi)存限制處理失敗等情形,這些都是ArcGIS無法比擬的,生產(chǎn)實踐中可極大地提高工作效率,為DEM后續(xù)的計算分析奠定基礎(chǔ)。