無人機DEM測點插值算法的礦區(qū)沉陷測繪技術(shù)

朱 濤

(河南省地質(zhì)局礦產(chǎn)資源勘查中心，鄭州 454000)

0 引言

DEM就是指數(shù)字高程模型，可以通過數(shù)字化模擬地形高程數(shù)據(jù)的方式，確定地面地形曲線的走向特征，簡單來說，就是利用有序數(shù)值陣列來表示實體地面模型。分辨率是DEM模型刻畫地形曲線的重要參考指標，其取值水平既影響數(shù)字化模擬手段的實施準確性，也影響有序數(shù)值陣列的實際排序形式[1]。無人機DEM測點插值算法是以DEM模型為基礎(chǔ)提出的測量數(shù)據(jù)處理思想，對于所得地面地形數(shù)據(jù)可以進行按需排序，并可以在不改變地形曲線走向特征的前提下，得到不同的參考坐標取值結(jié)果，從而在一次計算的過程中，得到更多的數(shù)據(jù)樣本坐標值。近年來，無人機DEM測點插值算法迅速發(fā)展，在工程建設(shè)、地質(zhì)勘探等多個領(lǐng)域中都得到了廣泛應(yīng)用[2]。此外，為保證數(shù)據(jù)坐標取值能夠準確表現(xiàn)出地面地形曲線的走向特征，還要求DEM模型不能在同一區(qū)域內(nèi)對數(shù)值陣列向量進行重復(fù)取值。

礦區(qū)沉陷就是指由煤礦地下采礦行為引起的地表塌陷現(xiàn)象，礦產(chǎn)資源被開采出來之后，礦區(qū)周圍巖體的力學平衡受到了破壞，而隨著這種破壞作用的持續(xù)累積，礦區(qū)地表就會表現(xiàn)出區(qū)域性塌陷的情況。在實際應(yīng)用過程中，如何根據(jù)礦區(qū)地表的實際沉陷程度，而制定具有針對性的測繪技術(shù)方案成為了一項亟待解決的難題。基于無人機激光雷達技術(shù)的沉陷監(jiān)測方法根據(jù)礦區(qū)沉陷特征的三維點云排列形式，建立多期數(shù)字高程模型，再通過DEM指標連續(xù)相減的方式，確定核心沉陷區(qū)地表的下陷程度[3]。基于HOG特征的檢測方法提取施工區(qū)域的邊緣特征，并借助訓練支持向量機原則，對相關(guān)特征參量進行分類，以便于實現(xiàn)對礦區(qū)沉陷程度的精準測量[4]。然而在上述兩種方法作用下，所得沉陷深度測量結(jié)果與真實沉陷深度之間的差值水平并不滿足實際應(yīng)用需求，故而按需所制定測繪技術(shù)的適用性也就相對較為局限。為解決上述問題，針對無人機DEM測點插值算法的礦區(qū)沉陷測繪技術(shù)展開研究。

1 無人機DEM測點插值算法

無人機DEM測點插值算法是完善礦區(qū)沉陷測繪技術(shù)實施方案的基礎(chǔ)，本章節(jié)將分別從DEM測點獲取、分維值測定、插值拐點求解3個方面，分析該算法的應(yīng)用流程。

1.1 DEM測點獲取

實施無人機測繪時，所得DEM測點數(shù)據(jù)包括平面位置、高程兩類信息參量[5]。其中，平面位置類的DEM測點數(shù)據(jù)對應(yīng)礦區(qū)地表及外部巖體的基本地質(zhì)特征，所得信息參量僅具有橫軸、縱軸兩個方向上的取值結(jié)果；高程類DEM測點數(shù)據(jù)對應(yīng)礦區(qū)無人機航攝圖像中的所有特征，所得信息參量具有橫軸、縱軸、空間軸3個方向上的取值結(jié)果[6]。

對于平面位置信息的定義式公為：

(1)

oX、oY分別表示橫軸、縱軸方向上的礦區(qū)地質(zhì)特征，δX表示橫軸方向上的測點參數(shù)，δY表示縱軸方向上的測點參數(shù)。

對于高程信息的定義式為：

(2)

oZ表示空間軸方向上的礦區(qū)地質(zhì)特征，δZ表示空間軸方向上的測點參數(shù)。

聯(lián)立公式(1)、公式(2)，可將無人機DEM測點獲取表達式定義為：

(3)

ΔO表示無人機航攝圖像中基本地質(zhì)特征的單位累積量，χ表示數(shù)據(jù)樣本測量參數(shù)。由于礦區(qū)的沉陷程度并不可控，所以oX、oY、oZ及δX、δY、δZ的取值存在相等的可能。

1.2 分維值測定

對于無人機DEM測點插值算法而言，分維值決定了測點參數(shù)之間的維數(shù)關(guān)系，且出于精確性考慮，分維值參數(shù)的取值應(yīng)屬于(0，1]的數(shù)值區(qū)間[7]。分維值等于“0”，表示所獲取DEM測點信息完全不能描述礦區(qū)的沉陷特征，在實際應(yīng)用過程中，即便是清晰度水平極低的無人機航攝圖像中也不可能完全不存在可用信息，所以這種情況并不會出現(xiàn)。分維值等于“1”，表示所獲取DEM測點信息能夠100%精準描述出礦區(qū)的沉陷特征，對于清晰度水平極高的無人機航攝圖像而言，這種表現(xiàn)情況是可能出現(xiàn)的，只是出現(xiàn)概率相對較小[8]。對于分維值測定表達式的求解滿足公式(4)。

(4)

1.3 插值拐點

插值拐點是指礦區(qū)無人機航攝圖像中DEM測點數(shù)據(jù)取值迅速發(fā)生改變的節(jié)點。所謂插值可以理解為安插在某一節(jié)點處的DEM測點數(shù)據(jù)，該數(shù)據(jù)樣本雖然具有當前節(jié)點處礦區(qū)沉陷地質(zhì)的全部特征，但由于其取值空間為無人機航攝圖像，所以測點數(shù)據(jù)結(jié)果僅具有運算作用，并不能完全代替礦區(qū)沉陷特征[9]。圖1反映了插值拐點的取樣原則。

圖1 插值拐點取樣原則

在無人機航攝圖像中，測繪方案實施方向并不一定與攝影方向保持一致，所以插值拐點之后DEM測點數(shù)據(jù)的取值既有增大也有縮小的可能(圖1中1曲線表示DEM測點數(shù)據(jù)取值的增大變化狀態(tài)，2曲線表示DEM測點數(shù)據(jù)取值的縮小變化狀態(tài))[10]。

設(shè)U0表示DEM測點數(shù)據(jù)的初始取值，Umax表示增大曲線上的DEM測點數(shù)據(jù)取值，Umin表示縮小曲線上的DEM測點數(shù)據(jù)取值，聯(lián)立公式(4)，可將插值拐點定義式表示為：

(5)

定義插值拐點時，如果不參考分維值測定條件，則有可能導致測繪數(shù)據(jù)與礦區(qū)地表真實沉陷情況出現(xiàn)偏差。

2 礦區(qū)沉陷區(qū)域地理模型

在無人機DEM測點插值算法的作用下，構(gòu)建礦區(qū)沉陷區(qū)域地理模型，還需根據(jù)像控點布設(shè)情況，定義空中三角，再聯(lián)合相關(guān)數(shù)據(jù)參量，求解徑向插值基函數(shù)表達式。

2.1 像控點布設(shè)

像控點布設(shè)就是在礦區(qū)沉陷航攝圖像中，按照無人機DEM測點插值算法所定義的測繪節(jié)點，一般來說，為避免插值拐點取樣結(jié)果無法真實反映出礦區(qū)地表的沉陷特征，在每一航攝方向上，都至少取兩個像控點作為實際布設(shè)位置[11]。如圖2為標準的像控點取樣結(jié)果。

圖2 像控點取樣

對圖2中的礦區(qū)沉陷航攝圖像進行背景去除處理，得到如圖3所示的像控點布設(shè)圖像。

根據(jù)圖3可知，對像控點進行布設(shè)時，應(yīng)在礦區(qū)沉陷航攝圖像中心定義一個標準控制區(qū)域，且該區(qū)域距離上、下兩端邊界的距離完全相等，距離左、右兩端邊界的距離也完全相等。航攝圖像上、下端存在兩組完全對稱的像控點，且這些節(jié)點的布設(shè)位置非常靠近圖像邊界[12]。左、右兩端像控點布設(shè)位置與標準控制區(qū)域的圓心等高，且這兩個節(jié)點的存在狀態(tài)也是完全對稱的。

2.2 空中三角

空中三角是輔助礦區(qū)沉陷測繪技術(shù)實施的必要條件，可以對已布設(shè)像控點對象進行分別取樣，再聯(lián)合所有數(shù)值結(jié)果，定義一個完整的空間數(shù)據(jù)集合，由于任意3個像控點都可以連接成一個三角形，所以該空間數(shù)據(jù)集合被簡稱為“空中三角”[13]。簡單來說，一個空中三角集合中所包含的全部數(shù)據(jù)樣本都描述礦區(qū)地表的沉陷特征，且這些數(shù)據(jù)的獲取完全依照無人機DEM測點插值算法。規(guī)定wX、wY、wZ分別表示橫軸、縱軸、空間中方向上的地表沉陷數(shù)據(jù)取值，其定義式如下：

(6)

其中，γX、γY、γZ分別表示橫軸、縱軸、空間中方向上的插值向量，RX、RY、RZ分別表示橫軸、縱軸、空間中方向上的測繪數(shù)據(jù)樣本。

在公式(6)的基礎(chǔ)上，設(shè)φ表示礦區(qū)航測圖像沉陷特征的全域測定參數(shù)，φ表示空間性標記參量，聯(lián)立上述物理量，可將空中三角表達式定義為：

(7)

為保證數(shù)據(jù)樣本空間的廣域性，一般不會對同一類型像控點中的沉陷特征進行重復(fù)取樣[14]。

2.3 徑向插值基函數(shù)

(8)

如果測繪區(qū)域的劃定范圍較大，利用無人機DEM測點插值算法對沉陷數(shù)據(jù)進行取值時，極有可能出現(xiàn)DEM測點數(shù)據(jù)與插值點高程信息不匹配的問題，此時空中三角集合中會出現(xiàn)不滿足徑向插值基函數(shù)表達式的數(shù)據(jù)樣本[16]。為避免上述情況的發(fā)生，在定義徑向插值基函數(shù)表達式之前，應(yīng)利用空中三角集合中的取樣信息進行重復(fù)計算，直至將非合理信息參量完全去除。

3 測繪實施方法設(shè)計

測繪實施方法設(shè)計就是按照空間坐標系轉(zhuǎn)換原則，對矢量數(shù)據(jù)進行疊加處理，再聯(lián)合像片傾角與旋偏角，確定測量精度評價指標的取值范圍，本章節(jié)將針對上述內(nèi)容展開研究。

3.1 空間坐標系轉(zhuǎn)換

空間坐標系轉(zhuǎn)換是按照無人機DEM測點插值算法定義原則，確定在測繪礦區(qū)發(fā)生沉陷時，橫向、縱向、空間向坐標軸之間的角度關(guān)系。轉(zhuǎn)換之前，橫向、縱向、空間向坐標軸之間的物理夾角均等于90°，這種定義模式只能適應(yīng)微傾斜情況下的測繪行為，當無人機航攝方向處于水平或豎直狀態(tài)時，至少有一個坐標軸在航攝圖像中的投影為一個點，而二維圖像中的點并不能描述出DEM測點數(shù)據(jù)的高程信息[17]。在轉(zhuǎn)換后的空間坐標系中，橫向、縱向、空間向坐標軸之間的夾角呈現(xiàn)出增大或縮小的數(shù)值變化狀態(tài)，在面對水平或豎直方向的測繪射線時，坐標軸在航攝圖像中的投影至少是具有一定程度的線段，而線段具有描述DEM測點數(shù)據(jù)高程信息的能力[18]。

對于空間坐標系轉(zhuǎn)換原則的定義滿足如下表達式：

(9)

3.2 矢量數(shù)據(jù)疊加

在實施測繪的過程中，矢量數(shù)據(jù)疊加就是將相似的DEM測點數(shù)據(jù)整合到一起，且處理后數(shù)據(jù)樣本的高程信息并不會發(fā)生改變[19]。對于矢量數(shù)據(jù)的疊加處理主要遵循如圖4所示的原則。

圖4 矢量數(shù)據(jù)疊加原則

無人機DEM測點插值算法規(guī)定，在空中三角集合中提取出的沉陷深度測量值保持矢量狀態(tài)。所謂矢量值是指同時具有方向特性與數(shù)值特性的數(shù)據(jù)樣本，在空間坐標系中，與無人機航測方向保持一致的數(shù)據(jù)樣本的取值方向為“+”，與無人機航測方向保持相反的數(shù)據(jù)樣本的取值方向為“-”[20]。兩個數(shù)值相同但取值方向相反的測繪數(shù)據(jù)樣本必須處于同一沉陷區(qū)域內(nèi)，但兩個取值方向相同、數(shù)值不同的測繪數(shù)據(jù)樣本有可能不屬于同一沉陷區(qū)域。設(shè)l1、l2、…、ln表示n個不相等的測繪數(shù)據(jù)樣本，κ1、κ2、…、κn分別表示與測繪數(shù)據(jù)樣本匹配的疊加系數(shù)，聯(lián)立公式(9)，可將矢量數(shù)據(jù)疊加原則表示為：

(10)

若存在數(shù)據(jù)樣本與DEM測點信息不完全匹配的情況，則表示疊加處理后，會剩余取值方向與數(shù)值水平都不固定的數(shù)據(jù)參量。

3.3 像片傾角與旋偏角

像片傾角是指無人機航攝圖像中測繪射線與沉陷深度所在地面水平線之間的物理夾角，如圖5中的ω。旋偏角是指測繪射線全旋角與像片傾角之間的差值，在圖5中可以表示為μ-ω。實施礦區(qū)沉陷測繪時，像片傾角的數(shù)值水平越大，就表示空間坐標系轉(zhuǎn)換處理過程中，坐標軸之間的物理夾角越大；而旋偏角越大，則表示DEM測點數(shù)據(jù)高程信息與真實沉陷數(shù)據(jù)之間的差值水平越大[21]。為實現(xiàn)對礦區(qū)沉陷深度的精準測量，應(yīng)同時控制像片傾角、旋偏角的實際取值。具體的像片傾角與旋偏角定義原則如圖5所示。

圖5 像片傾角與旋偏角定義

在無人機航攝圖像中，全旋角數(shù)值一般等于180°，而像片傾角數(shù)值則小于90°，因此旋偏角取值屬于90°-180°的數(shù)值區(qū)間[22]。

3.4 測量精度評價指標

測量精度評價指標影響無人機DEM測點插值算法對礦區(qū)沉陷程度的測繪處理能力，由于DEM測點數(shù)據(jù)高程信息始終保持定值狀態(tài)，所以即便在測繪方案不符合實際應(yīng)用需求的情況下，測量精度評價指標的計算數(shù)值也不會表現(xiàn)出波動變化行為[23-24]。如果沉陷區(qū)域面積保持不變，那么像片傾角與旋偏角數(shù)值也可以保持定值狀態(tài)，當前情況下求解所得的測量精度評價指標也就能夠符合測繪礦區(qū)的實際應(yīng)用需求。對于測量精度評價指標的計算滿足如下表達式：

(11)

4 實例分析

本次實驗的主要目的是分析無人機DEM測點插值算法的礦區(qū)沉陷測繪技術(shù)、基于無人機激光雷達技術(shù)的沉陷監(jiān)測方法、基于HOG特征的檢測方法的實用性能力，在不考慮其他干擾條件的情況下，設(shè)計如下對比實驗。

4.1 礦區(qū)情況

選擇位于中海拔地區(qū)的某土石礦區(qū)作為實驗環(huán)境，利用遙感技術(shù)對該區(qū)域的地質(zhì)特征進行檢測，詳情如圖6所示。

圖6 礦區(qū)地質(zhì)遙感圖像

圖6中部、下部兩塊小型區(qū)域的顏色明顯比其他遙感區(qū)域的顏色更深，根據(jù)遙感成像特點可知，這兩個區(qū)域為主要沉陷區(qū)，在上述兩區(qū)域外側(cè)劃定一片區(qū)域作為本次測繪的具體實驗區(qū)。為避免開采行為對礦區(qū)地質(zhì)造成嚴重影響，實驗過程中必須借助遙感技術(shù)對礦區(qū)沉陷程度進行實時監(jiān)測。

4.2 沉陷測繪點取樣

在實驗區(qū)域內(nèi)選擇8個不重復(fù)、深度水平也不相同的測繪節(jié)點作為實驗對象，并對每一節(jié)點處的地表沉陷深度進行測量，詳情見表1。

表1 測繪點沉陷深度 m

分析表1可知，該實驗區(qū)域內(nèi)，地表沉陷程度的差異性較大，4號測繪點的最大沉陷深度與7號測繪點的最小沉陷深度之間的數(shù)值差達到了22.51 m。

4.3 實驗結(jié)果分析

分別利用無人機DEM測點插值算法的礦區(qū)沉陷測繪技術(shù)、基于無人機激光雷達技術(shù)的沉陷監(jiān)測方法、基于HOG特征的檢測方法對8個測繪節(jié)點處的地表沉陷深度進行測量，其中第一組為實驗組方法、第二組為對照A組方法、第三組為對照B組方法。

圖7記錄了實驗組、對照組礦區(qū)地表沉陷深度的具體實驗數(shù)值。

圖7 礦區(qū)地表沉陷深度實驗數(shù)值

分析圖7可知，第7個測繪節(jié)點處，對照A組、對照B組沉陷深度測量值均為零，在精度分析方面的參考價值不大，因此在后續(xù)分析過程中，必須去除該測量結(jié)果。實驗組方法所測得的礦區(qū)地表沉陷深度均值水平相對較高，對照A組方法次之，而對照B組測量值最小。

取第4測繪節(jié)點處的最大實驗數(shù)值進行計算，可知沉陷深度測量結(jié)果與礦區(qū)真實沉陷深度之間的差值為：

(12)

其中，M1表示實驗組差值，M2表示對照A組差值，M3表示對照B組差值，單位均為m。

根據(jù)公式(12)可知，實驗組方法作用下，沉陷深度測量結(jié)果與礦區(qū)真實沉陷深度之間的差值最小，對照A組方法次之，對照B組方法最大。

由上述實驗結(jié)果，可得出以下結(jié)論：

1)基于無人機激光雷達技術(shù)的沉陷監(jiān)測方法的應(yīng)用，對于沉陷深度測量結(jié)果與礦區(qū)真實沉陷深度之間差值水平的控制能力有限，不足以實現(xiàn)對礦區(qū)沉陷深度的精準測量；

2)相較于基于無人機激光雷達技術(shù)的沉陷監(jiān)測方法，基于HOG特征的檢測方法的測量能力相對較強，能夠適當控制沉陷深度測量結(jié)果與礦區(qū)真實沉陷深度之間的差值水平，但依然達不到完善測繪技術(shù)實施方案的實際應(yīng)用需求；

3)無人機DEM測點插值算法的礦區(qū)沉陷測繪技術(shù)的應(yīng)用，可有效控制沉陷深度測量結(jié)果與礦區(qū)真實沉陷深度之間的差值水平，使其計算數(shù)值不超過1 m，不但實現(xiàn)了對礦區(qū)沉陷深度的精準測量，還可以對測繪技術(shù)具體實施方案進行完善，與無人機激光雷達技術(shù)的沉陷監(jiān)測方法、基于無人機激光雷達技術(shù)的沉陷監(jiān)測方法相比，更符合實際應(yīng)用需求。

5 結(jié)束語

新型礦區(qū)沉陷測繪技術(shù)在無人機DEM測點插值算法的基礎(chǔ)上，通過計算分維值測定結(jié)果的方式，確定插值拐點取樣條件，又聯(lián)合徑向插值基函數(shù)表達式，定義空間坐標系轉(zhuǎn)換原則，從而推導出測量精度評價指標的實際取值范圍。實用性方面，無人機DEM測點插值算法的礦區(qū)沉陷測繪技術(shù)的作用下，沉陷深度測量結(jié)果與礦區(qū)真實沉陷深度之間的差值水平得到了有效控制，有利于實現(xiàn)對礦區(qū)沉陷深度的精準測量，對于完善測繪技術(shù)的具體實施方案，可以起到一定的促進性影響作用。