朱 濤

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

0 引言

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

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

1 無人機(jī)DEM測點(diǎn)插值算法

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

1.1 DEM測點(diǎn)獲取

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

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

(1)

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

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

(2)

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

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

(3)

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

1.2 分維值測定

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

(4)

1.3 插值拐點(diǎn)

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

圖1 插值拐點(diǎn)取樣原則

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

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

(5)

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

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

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

2.1 像控點(diǎn)布設(shè)

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

圖2 像控點(diǎn)取樣

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

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

2.2 空中三角

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

(6)

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

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

(7)

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

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

(8)

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

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

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

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

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

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

(9)

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

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

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

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

(10)

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

3.3 像片傾角與旋偏角

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

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

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

3.4 測量精度評價(jià)指標(biāo)

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

(11)

4 實(shí)例分析

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

4.1 礦區(qū)情況

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

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

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

4.2 沉陷測繪點(diǎn)取樣

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

表1 測繪點(diǎn)沉陷深度 m

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

4.3 實(shí)驗(yàn)結(jié)果分析

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

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

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

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

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

(12)

其中，M1表示實(shí)驗(yàn)組差值，M2表示對照A組差值，M3表示對照B組差值，單位均為m。

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

由上述實(shí)驗(yàn)結(jié)果，可得出以下結(jié)論：

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

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

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

5 結(jié)束語

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