楊雪峰，徐漢超，范 鵬

(1.遼寧江海水利工程有限公司，遼寧 沈陽(yáng) 110006；2.沈陽(yáng)市勘察測(cè)繪研究院有限公司，遼寧 沈陽(yáng) 110004；3.遼寧省水利水電勘測(cè)設(shè)計(jì)研究院有限責(zé)任公司，遼寧 沈陽(yáng) 110006)

無(wú)人機(jī)測(cè)量獲取的DSM數(shù)據(jù)為包含了地表建筑物、植被等高度的地表高程模型，需通過(guò)優(yōu)化處理后方可作為真實(shí)數(shù)據(jù)使用，目前尚無(wú)有效技術(shù)手段將DSM數(shù)據(jù)直接轉(zhuǎn)化為DEM數(shù)據(jù)。本文采用雙線性內(nèi)插計(jì)算原理，通過(guò)lisp程序?qū)o(wú)人機(jī)獲取的DSM、DOM數(shù)據(jù)相結(jié)合，建立可靠數(shù)字地形模型，真實(shí)還原目標(biāo)區(qū)域的地形情況，為水土流失治理項(xiàng)目設(shè)計(jì)工作提供真實(shí)可靠的基礎(chǔ)數(shù)據(jù)，保證水土保持措施的精準(zhǔn)實(shí)施，解決以往項(xiàng)目中使用投影面積乘系數(shù)所產(chǎn)生的誤差。

1 基于無(wú)人機(jī)影像的坡度分析與項(xiàng)目區(qū)域劃定

目前水土保持防治項(xiàng)目區(qū)域的選擇和面積的量算主要在地面上進(jìn)行，無(wú)法全面觀測(cè)地形地勢(shì)，難以精準(zhǔn)選擇治理區(qū)域。隨著無(wú)人機(jī)技術(shù)的發(fā)展，可以通過(guò)無(wú)人機(jī)拍照和視頻攝制定性、無(wú)人機(jī)測(cè)量功能定量的方法和手段，實(shí)現(xiàn)全面多樣、大面覆蓋的方式獲取目標(biāo)區(qū)域的數(shù)據(jù)，從而實(shí)現(xiàn)精準(zhǔn)實(shí)施水土流失治理方案[1]。

在水土流失治理工作中，無(wú)人機(jī)航空攝影得到的DOM和DSM數(shù)據(jù)經(jīng)初步坡度分析后，成果僅用作項(xiàng)目區(qū)位置選取。故本次以DSM數(shù)據(jù)作為分析樣本，樣本數(shù)據(jù)獲取于2021年5月15日的實(shí)地?zé)o人機(jī)測(cè)量，航攝面積為0.63km2，共239張照片。沿東西方向布設(shè)航線，沿航線方向等距間隔拍攝，航向重疊度80%，旁向重疊度75%，航高100m，選用1英寸2 000萬(wàn)像素CMOS傳感器，地面分辨率5.11cm/2.01in。獲取了分辨率為1m的數(shù)字地表模型(DSM)和分辨率為0.3 m的數(shù)字正射影像(DOM)，如圖1所示。

圖1 樣區(qū)數(shù)據(jù)情況

利用DSM數(shù)據(jù)成果，提取柵格數(shù)據(jù)中每個(gè)像元的坡度值，利用像元坡度值圖層結(jié)合無(wú)人機(jī)影像數(shù)據(jù)獲取的項(xiàng)目區(qū)坡度、坡向、植被條件和土壤侵蝕情況，對(duì)項(xiàng)目區(qū)進(jìn)行坡度分級(jí)分類，選取項(xiàng)目區(qū)實(shí)施區(qū)域[2]。區(qū)域選取情況如圖2所示。

圖2 項(xiàng)目區(qū)選取結(jié)果

在水土保持項(xiàng)目設(shè)計(jì)和施工的具體實(shí)施階段，精確的地形地貌數(shù)據(jù)可以保證水土保持設(shè)計(jì)和施工方案具有針對(duì)性，故需要將DSM數(shù)據(jù)轉(zhuǎn)化為更可靠的DEM數(shù)據(jù)。

2 影像數(shù)據(jù)特性分析

遙感數(shù)據(jù)具有一定的精度限制，而無(wú)人機(jī)采集的地形數(shù)據(jù)可以用高分辨率的DSM數(shù)據(jù)表達(dá)，但DSM數(shù)據(jù)區(qū)別于DEM數(shù)據(jù)，需分離地面上的建筑物和植被方可得到真實(shí)的DEM數(shù)據(jù)模型。

早期主流的數(shù)據(jù)處理辦法包括數(shù)學(xué)形態(tài)學(xué)處理、最小二乘內(nèi)插和TIN加密。ISPRS組織的實(shí)驗(yàn)結(jié)果表明，這3種方法對(duì)于比較光滑的地形，能夠較好的剔除建筑物和植被，但在地形粗糙復(fù)雜的地區(qū)存在一定問(wèn)題[3]。另一方面，大部分算法對(duì)濾波和各種闕值等參數(shù)的變化十分敏感，且地形變化本身具有復(fù)雜和不可預(yù)測(cè)性，當(dāng)全局點(diǎn)云數(shù)據(jù)采用同一套計(jì)算參數(shù)時(shí)，必然會(huì)產(chǎn)生較大誤差。因此，當(dāng)前項(xiàng)目，可以采用局部地區(qū)雙線性內(nèi)插計(jì)算內(nèi)插高程值的方式，在小范圍內(nèi)對(duì)地形變化進(jìn)行定量表達(dá)，使計(jì)算參數(shù)分區(qū)域適應(yīng)調(diào)整，解決DEM數(shù)字模型建立的精度問(wèn)題。

3 DEM數(shù)據(jù)獲取實(shí)現(xiàn)方法

在無(wú)人機(jī)拍攝的DSM數(shù)據(jù)中，地表裸露的區(qū)域可以大量提取真實(shí)地形的高程數(shù)據(jù)，而植被覆蓋的區(qū)域是影響DEM數(shù)據(jù)成果精度的主要原因。目前尚無(wú)直接剔除植被等數(shù)據(jù)的辦法，故研究采用局部小區(qū)域雙線性內(nèi)插計(jì)算內(nèi)插點(diǎn)的辦法。技術(shù)流程圖如圖3所示。

圖3 技術(shù)路線圖

3.1 不考慮地形特征的內(nèi)插方式

如圖4所示，局部?jī)?nèi)插高程區(qū)域有一條山脊線，需要在山脊線左側(cè)內(nèi)插點(diǎn)，當(dāng)采用不考慮地形特征的線性內(nèi)插，內(nèi)插高程值將大大小于真實(shí)值。因此，為了提高DEM數(shù)據(jù)的精度，在內(nèi)插高程點(diǎn)的時(shí)候應(yīng)加以考慮地形特征[4]。

圖4 不考慮地形特征的內(nèi)插

目前，工程中常用CASS軟件內(nèi)插高程值的解決辦法有4種，分別為不擬合的折線法、張力樣條擬合法、三次B樣條擬合法和SPLine擬合法。其中折線和SPLine方法使用環(huán)境局限性較大，一般不采用；張力樣條擬合法和三次B樣條擬合法精度較高[5]。

本次通過(guò)張力樣條擬合法和三次B樣條擬合法分別對(duì)樣本數(shù)據(jù)進(jìn)行全局解算內(nèi)插，并對(duì)內(nèi)插成果與真實(shí)值比對(duì)，高差最大差值分別為0.351m和0.499m，中誤差分別為0.35m和0.48m。

3.2 局部雙線性內(nèi)插方法

3.2.1局部雙線性內(nèi)插原理

首先在一維中理解雙線性內(nèi)插計(jì)算，將參考點(diǎn)構(gòu)成的地形趨勢(shì)線垂直投影到通過(guò)內(nèi)插點(diǎn)的線上，得到一個(gè)交點(diǎn)，再通過(guò)山脊線上參考點(diǎn)趨勢(shì)線的校正，得到的交點(diǎn)均值作為內(nèi)插點(diǎn)的高程值，如圖5所示。

圖5 雙線性內(nèi)插一維示意圖

雙線性內(nèi)插可以理解為兩次或多次的線性內(nèi)插，如圖6所示，已知4個(gè)參考點(diǎn)坐標(biāo)點(diǎn)分別為(a1，b1)、(a1，b1+1)、(a1+1，b1+1)、(a1+1，b1)，H軸表示數(shù)據(jù)高程。在A面部分，變量為a1，而b1為常數(shù)，連線H(a1，b1)和H(a1+1，b1)做一維線性內(nèi)插，求得H(a，b1)，同理在B面和C面部分可以分別求得H(a，b1+1)和H(a，b)。

圖6 雙線性內(nèi)插計(jì)算示意圖

其函數(shù)形式為：

Z=a0+a1x+a2x+a3xy

(1)

設(shè)4個(gè)參考點(diǎn)為p1(x1，y1，z1)，p2(x2，y2，z2)，p3(x3，y3，z3)和p4(x4，y4，z4)，代入下式：

(2)

當(dāng)求得4個(gè)參數(shù)后，將內(nèi)插點(diǎn)P(X，Y)的(x，y)代入公式(1)，便可求得p(z)。已知參數(shù)帶入越多，內(nèi)插高程值精度越高。

3.2.2局部雙線性內(nèi)插在lisp程序中的實(shí)現(xiàn)

以樣本DSM數(shù)據(jù)中提取的真實(shí)地面高程點(diǎn)作為原始數(shù)據(jù)，意在通過(guò)局部雙線性內(nèi)插的方法獲取被植被覆蓋處的地面高程數(shù)據(jù)。

在lisp程序中根據(jù)使用功能設(shè)計(jì)若干子程序，主要子程序包括光標(biāo)撲獲高程點(diǎn)的信息、消除字符串中的空格、對(duì)坐標(biāo)文件排序、判斷點(diǎn)與點(diǎn)之間的位置、角度法判斷點(diǎn)位、高程內(nèi)插、線性內(nèi)插、雙線性內(nèi)插、求雙線性內(nèi)插交點(diǎn)、判斷高程值h是否在h1和h2之間、追蹤和創(chuàng)建Cass高程點(diǎn)、檢查字體“HZ”是否存在、檢查是否存在高程點(diǎn)圖塊定義等。

主程序部分如下：

其中參數(shù)datasj為數(shù)據(jù)文件、gcsj為兩端點(diǎn)高程、dgxsj為兩條等高線、gcsj為單個(gè)高程點(diǎn)、xxjlsj為距離間隔、jlsj為距離。

為了提高內(nèi)插精度，在進(jìn)行內(nèi)插工作前，在DSM數(shù)據(jù)中標(biāo)繪出山脊線、山谷線等可以表達(dá)出地形驟變的地形線。

3.2.3插值精度分析

為了驗(yàn)證插值精度的準(zhǔn)確性，依據(jù)地形起伏條件的不同選取兩類樣本區(qū)域，分別抽取50個(gè)內(nèi)插點(diǎn)作為驗(yàn)證數(shù)據(jù)集，選取的已知高程點(diǎn)包括了平緩區(qū)域、地形起伏較大區(qū)域等具有不同代表性數(shù)據(jù)區(qū)域[6]。對(duì)比分析插值點(diǎn)和實(shí)測(cè)高程值得到不同區(qū)域數(shù)據(jù)精度對(duì)比結(jié)果，高程點(diǎn)位中誤差離散分布如圖7所示。

圖7 內(nèi)插點(diǎn)位中誤差離散圖

在地形較平坦的區(qū)域插值精度最高，樣本數(shù)據(jù)最大中誤差僅為0.16m，測(cè)區(qū)地形起伏較大的區(qū)域插值效果較為優(yōu)異，樣本數(shù)據(jù)最大中誤差僅為0.26m，優(yōu)于1m等高線內(nèi)插點(diǎn)高程精度小于1/3等高距、1/2等高距的規(guī)范要求。

3.3 精度對(duì)比

通過(guò)局部雙線性內(nèi)插高程值精度與張力樣條擬合法、三次B樣條擬合法全局解算內(nèi)插方式的成果對(duì)比，局部雙線性內(nèi)插成果精度高于另兩種方式獲取的高程值。

4 利用DEM數(shù)據(jù)精準(zhǔn)布局水土保持實(shí)施方案

根據(jù)DEM和DOM數(shù)據(jù)情況，本次治理片區(qū)內(nèi)坡耕地上游分布有林地、疏林地，水源涵養(yǎng)能力較強(qiáng)；部分坡耕地斑塊結(jié)合水平槽整地措施，增加了坡面蓄水能力[7]。同時(shí)，坡耕地上游雨水部分向兩側(cè)排入溝道，很大程度減少降雨徑流對(duì)坡耕地的沖刷，可以取消在坡耕地上游布設(shè)截排水溝設(shè)計(jì)。結(jié)合DEM的坡度分析結(jié)果，具體實(shí)施方案見(jiàn)表1。

表1 坡度分級(jí)與水土保持實(shí)施方案

對(duì)6°～15°坡耕地采取保土耕作和坡式梯田措施，無(wú)法修建梯田的坡耕地，耕作方向要基本沿等高線，以攔蓄部分降雨。

對(duì)15°～25°的荒山荒坡地進(jìn)行水平槽整地，減少荒坡地水土流失，并水平槽整地基礎(chǔ)上結(jié)合產(chǎn)業(yè)調(diào)整適時(shí)栽植樹(shù)木[8]。

對(duì)于>25°對(duì)現(xiàn)有的天然林、人工林、疏林地及荒地等進(jìn)行封山育林，提高林地覆被率，增強(qiáng)保水保土和抗蝕能力。明確專人管理，增設(shè)圍欄和標(biāo)志牌[9]。

5 結(jié)語(yǔ)

(1)以樣本數(shù)據(jù)為基礎(chǔ)，對(duì)DSM數(shù)據(jù)進(jìn)行局部雙線性內(nèi)插，并將內(nèi)插結(jié)果與真實(shí)測(cè)量值和全局解算內(nèi)插獲取高程方式的成果對(duì)比，該方法可以得到更趨近真實(shí)地表情況的高程值。

(2)通過(guò)局部雙線性內(nèi)插法優(yōu)化處理后的DEM數(shù)據(jù)對(duì)項(xiàng)目區(qū)坡度、坡向等因子提取精度更高，對(duì)水土保持方案的實(shí)施有更好的幫助作用。此外，通過(guò)三維數(shù)據(jù)模型量取項(xiàng)目區(qū)的表面積，可以準(zhǔn)確的計(jì)算項(xiàng)目實(shí)施面積。

(3)本文算法旨在解決局部區(qū)域的內(nèi)插精度，但在使用過(guò)程中，識(shí)別山脊線和山谷線等地性線時(shí)注意高程參考點(diǎn)的篩選，地形轉(zhuǎn)折處需要主觀判斷，今后還需進(jìn)一步優(yōu)化。