數字孿生基礎數據底板DEM修正方法研究

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

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

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

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

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

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

1 材料與方法

1.1 研究區概況和數據來源

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

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

1.2 計算原理

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

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

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

a)像元位置

(1)

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

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

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

(2)

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

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

a)1個像元無值

1.3 研究技術路線及步驟

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

圖3 研究技術路線

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

步驟二將處理好的等高線和高程點數據作為創建TIN的輸入數據,其中等高線要素作為contour,高程點要素作為masspoint。

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

步驟四重復步驟一、二完成所有單幅地形圖轉為DEM后將所有DEM進行鑲嵌。

步驟五對鑲嵌后DEM進行NoData值修正處理。

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

1.4 分析指標及評定標準

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

(3)

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

當m

參照GB 50026—2020《工程測量標準》條文說明5.1.3中,地形圖的基本等高距,是以等高線的高程中誤差的經驗公式驗算,數字高程模型格網間距的選取及格網點高程中誤差應符合表1的規定。

表1 數字高程模型格網間距的選取及格網點高程中誤差 單位:m

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

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

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

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

2 結果與討論

2.1 修正方法的適用性

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

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

表3 空間自相關分析結果數據

圖4 比較值DEM空間自相關分析圖莫蘭指數

2.2 修正結果

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

a)NoData像元空值

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

圖6 修正值與比較值數據點比較

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

2.3 中誤差分析

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

圖8 中誤差分布

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

(4)

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

表4 單位成果質量評定等級

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

3 結論

本文通過Python調用ArcGIS相應的工具來批量修正某市某研究區域的地形轉DEM后接縫處存在的高程NoData值,其主要結論如下。

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

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

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