無人機DEM測點插值算法的礦區沉陷測繪技術

朱 濤

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

0 引言

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

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

1 無人機DEM測點插值算法

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

1.1 DEM測點獲取

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

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

(1)

oX、oY分別表示橫軸、縱軸方向上的礦區地質特征，δX表示橫軸方向上的測點參數，δY表示縱軸方向上的測點參數。

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

(2)

oZ表示空間軸方向上的礦區地質特征，δZ表示空間軸方向上的測點參數。

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

(3)

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

1.2 分維值測定

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

(4)

1.3 插值拐點

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

圖1 插值拐點取樣原則

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

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

(5)

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

2 礦區沉陷區域地理模型

在無人機DEM測點插值算法的作用下，構建礦區沉陷區域地理模型，還需根據像控點布設情況，定義空中三角，再聯合相關數據參量，求解徑向插值基函數表達式。

2.1 像控點布設

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

圖2 像控點取樣

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

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

2.2 空中三角

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

(6)

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

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

(7)

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

2.3 徑向插值基函數

(8)

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

3 測繪實施方法設計

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

3.1 空間坐標系轉換

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

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

(9)

3.2 矢量數據疊加

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

圖4 矢量數據疊加原則

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

(10)

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

3.3 像片傾角與旋偏角

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

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

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

3.4 測量精度評價指標

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

(11)

4 實例分析

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

4.1 礦區情況

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

圖6 礦區地質遙感圖像

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

4.2 沉陷測繪點取樣

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

表1 測繪點沉陷深度 m

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

4.3 實驗結果分析

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

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

圖7 礦區地表沉陷深度實驗數值

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

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

(12)

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

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

由上述實驗結果，可得出以下結論：

1)基于無人機激光雷達技術的沉陷監測方法的應用，對于沉陷深度測量結果與礦區真實沉陷深度之間差值水平的控制能力有限，不足以實現對礦區沉陷深度的精準測量；

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

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

5 結束語

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