無人機遙感技術在北京首云鐵礦儲量動態(tài)監(jiān)測中的應用

向 杰， 陳建平， 李 詩， 賴自力， 黃浩中， 劉 靜， 謝 帥

(1.中國地質科學院礦產資源研究所自然資源部成礦作用與資源評價重點實驗室，北京 100037； 2.中國地質大學(北京)地球科學與資源學院，北京 100083)

0 引言

無人機(unmanned aerial vehicle，UAV)遙感技術是無人機駕駛飛行器技術、遙感傳感器技術、遙測遙控技術、通信技術、POS定位定姿技術和GPS差分定位技術等一系列技術的組合[1]。近年來，UAV遙感技術廣泛應用于國土資源環(huán)境領域，成為國土資源遙感領域的研究熱點之一。例如，利用UAV遙感開展土地資源調查與分類； 進行林業(yè)資源分布調查與資源總量估算； 開展環(huán)境監(jiān)測與污染追蹤調查等[2-3]。由于UAV遙感和礦山開采條件的限制，在礦山儲量動態(tài)監(jiān)測方面應用較少[4-5]。如何快速、準確、廉價地進行礦山監(jiān)測一直是國土資源遙感領域的研究難點。

本文以北京市密云區(qū)首云鐵礦為研究區(qū)，利用無人機攝影測量技術(UAV structure from motion，SfM)獲取不同時相的高精度數字地表模型(digital surface model，DSM)，快速準確計算礦山動用儲量，從而實現對儲量的動態(tài)監(jiān)控。與傳統(tǒng)方法相比，該方法具有準確、快速、性價比高的特點，對于中小范圍地區(qū)的礦山監(jiān)測具有重要的推廣價值。

1 研究區(qū)概況

首云鐵礦屬于密云鐵礦的一部分，位于北京市密云區(qū)巨各莊鎮(zhèn)境內，距北京80 km，屬于較典型的丘陵地帶。首云鐵礦始建于1959年，1970年建成投產，是集采礦、選礦、豎爐焙燒連續(xù)生產工藝為一體的國有中型冶金礦山企業(yè)，曾隸屬于北京市冶金局和首鋼總公司，是原首鋼鐵礦冶金部重點礦山企業(yè)之一。

首云鐵礦為露天開采，該礦區(qū)開采面積大約為1.5 km2，平均海拔為200 m，礦區(qū)植被稀少(圖1)，這些條件都為UAV遙感調查提供了良好的環(huán)境。在2014年8月—2016年10月期間，該礦山一直處于露天開采狀態(tài)，并且該礦山礦石體重和礦石品位比較均勻，通過計算地表體積變化即可計算出礦山動用儲量，這為UAV進行礦山儲量監(jiān)測提供了前提條件。

圖1 北京市密云區(qū)首云鐵礦位置Fig.1 Location of Shouyun iron mine in Miyun district, Beijing

2 UAV數據獲取與處理

2.1 UAV數據獲取

本次工作分別于2014年8月和2016年10月對首云鐵礦進行了UAV野外調查。2次UAV野外調查期間的天氣晴朗無風，非常適合UAV飛行。另外，盡量選擇正午期間進行飛行，減少了山體影像對結果的影響。2次UAV調查都采用臺灣觀天科技生產的X5型小型電動固定翼UAV，該UAV翼展為1.2 m，起飛重量為1.5 kg，有效載荷為300 g，續(xù)航時間可達40 min。UAV搭載的自駕儀和GPS可以實現按照規(guī)劃航線自動飛行和定時、定距拍照； 機載POS系統(tǒng)可同步記錄拍照時UAV姿態(tài)以用于數據后續(xù)處理； 搭載的相機為Sony QX100鏡頭，其基本參數如表1。

表1 UAV搭載數碼相機基本參數Tab.1 Parameters of camera equipped on UAV

2次UAV調查飛行高度都設定為海拔450 m，設定航向重疊率80%，旁向重疊率60%。2014年8月的調查中獲取了499張帶POS信息的高清影像，實際航向重疊率81%，旁向重疊率62%； 2016年10月的調查過程中獲取了674張帶POS信息的高清影像，實際航向重疊率80%，旁向重疊率45%，并利用DGPS在特征明顯且無變化的區(qū)域測定了7個控制點(ground control point，GCP)。UAV飛行航線如圖2所示。

圖2 UAV野外調查與數據獲取航線Fig.2 Field surveys and data acquisition lines of UAV

2.2 UAV數據處理

UAV數據處理主要是為了得到高精度的DSM和高空間分辨率的數字正射影像圖(digital orthophoto map，DOM)。近年來市場上涌現了一大批UAV遙感影像處理軟件，如國內的DPGrid，PixelGrid和DPMatrix等，國外的PhotoScan，Pix4Dmapper和Inpho等[6]。

本次DOM的獲取采用PhotoScan V1.2.5進行流程化處理，生成了高清的DOM。為方便后續(xù)計算，DSM獲取采用如下步驟： ①在PhotoScan V1.2.5中預處理，其過程包括添加影像—加載POS點信息—對齊影像—生成網格—添加GCP—優(yōu)化相機對齊方式—建立密集點云—生成點云數據； ②將導出的點云數據導入到cloud compare中，去除掉植被或一些噪聲點； ③將處理好的點云數據導入到ArcGIS中，采用自然鄰域法進行插值，生成空間分辨率為1 m的DSM[7]。

要進行地表變化的計算，2期DSM之間的配準極為重要，因為它會直接影響最后的研究精度。為此，本文首先利用DGPS測量的7個GCP數據，在ArcGIS中對2016年數據進行配準； 其次，對比2014年與2016年的DOM，在已配準完成的2016年數據上，選取特征明顯且沒有發(fā)生改變的點作為2014年的GCP； 最后，利用選取出的10個GCP配準2014年數據。采用這種方法不僅很好地解決了2014年控制點缺失的問題，也使得2期數據互相匹配，提高了后期地表變化計算的精度。

2014年和2016年2期UAV攝影測量獲取的點云及DSM如圖3所示。

(a) 2014年點云 (b) 2016年點云

(c) 2014年DSM (d) 2016年DSM

圖32期UAV攝影測量獲取的點云與DSM

Fig.3CloudpointandDSMoftwoperiodsusingUAVdata

通過以上處理過程，獲得了2期高清DOM與配準過后的1 m空間分辨率DSM。2016年DSM誤差采用已經測定的7個GCP來計算均方根誤差(root mean square error，RMSE)得到。由于2014年沒有測量GCP，因此基于配準完成的2016年數據，在沒有發(fā)現變化的穩(wěn)定區(qū)域挑選出10個GCP來計算DSM誤差。2014年8月DOM的空間分辨率可達到0.081 m, 獲得的DSM誤差為0.44 m； 2016年10月DOM的空間分辨率達到0.072 m，獲得的DSM誤差為0.24 m。2期DSM的誤差都屬于亞m級，能滿足進行礦山儲量估算的精度要求。

3 動用儲量計算

地形變化算法(DSM of Difference，DoD)最早運用于地形地貌領域。20世紀90年代早期，瑞士地貌學家Lane等便開始利用多期DSM相減來定量描述地形的變化[8]。近年來，獲取DSM (或DEM)的技術手段不斷豐富，如： 衛(wèi)星立體像對、合成孔徑雷達、SfM和激光雷達等，DSM精度在不斷提高，DoD算法也不斷改進完善[9]。本研究通過SfM，分別獲取了2014年與2016年首云露天礦山的DSM，采用2期DSM相減的方式即可表示礦山的開采量。由于該礦山礦石密度與礦石品位都比較均勻，與礦山開采量相乘即可得到儲量變化，進而動態(tài)監(jiān)測露天礦山的儲量變化。

但是，當用DSM來表示實際的地表模型時存在不確定性δZ，即

Z真實=ZDSM±δZ,

(1)

式中：Z真實表示實際的高程；ZDSM表示DSM所表示的高程；δZ表示對于單幅DSM的誤差，而要計算地形變化必須要考慮到2期DSM的傳播誤差，通常δZ受到很多因素的影響，如測量誤差、采樣密度、類型和插值方法等。

2003年英國地貌學家Brasington等對于傳播誤差δDoD提出的公式[10]被廣為接受，即

(2)

式中：δznew表示較新DSM的誤差；δzold表示較老DSM的誤差。

為了更加準確地計算礦山動用儲量，消除傳播誤差對于DoD結果的影響顯得尤為重要，本文采用了最小檢出限法和轉換概率法2種方法進行計算。

3.1 最小檢出限法

最小檢出限法是將傳播誤差δDoD作為最小檢出限，即2期DSM相減，高程值變化小于δDoD視為是由于誤差造成，不是由于礦山開采造成的，只有當高程值變化大于δDoD時，才視其為真正的開采所致[8]。2014年與2016年這2期DSM直接相減得到原始DoD如圖4(a)所示。如前所述，2014年DSM誤差為0.44 m，2016年DSM誤差為0.24 m，根據公式(2)可求得傳播誤差為0.5 m。將傳播誤差0.5 m作為最小檢出限法的DoD結果如圖4(b)所示。

(a) 原始DoD (b) 最小檢出限法DoD

圖4最小檢出限法計算地形變化結果

Fig.4ResultsofDSMofdifferencebyusingthelimitofdetectionmethod

通過對比圖4(a)和(b)可以看出，最小檢出限法所摒棄的信息主要分布在礦區(qū)中部和北部。從DOM上可以清晰地看到，該區(qū)域主要是道路與廠房，確實屬于穩(wěn)定非開采區(qū)域。礦山的主采區(qū)分布在礦坑南緣，尤其是中南部區(qū)域，開采造成的高程變化達到了80余m。礦山的堆積區(qū)主要分布在礦坑中部以及西北部，高程變化在20 m左右。

3.2 轉換概率法

同樣也可以將DSM的不確定性轉化成概率值，通過用戶自主定義置信區(qū)間的方式來去除不確定性對最終結果的影響[10-11]。假設傳播誤差的臨界誤差為U，則式(2)可以轉換為式(3)， 基于t檢驗來確定一個置信區(qū)間。

(3)

(4)

式中|Znew-Zold|可以由原始DoD的絕對值得到。

因此，可以將原始DoD轉化成如圖5(a)所示的T值分布圖。然后，可以通過將t統(tǒng)計量與其累積分布函數相關聯(lián)來計算純粹由于誤差而引起的高程變化，利用t的累計分布函數將其轉化為高程變化概率分布，通過設定合適的概率作為閾值來更加準確地估計動用儲量。本研究按照慣例選取95%作為置信區(qū)間，來求取轉換概率法DoD，結果如圖5(b)所示。

(a) T值分布 (b) 轉換概率法DoD

圖5轉換概率法計算地形變化結果

Fig.5ResultsofDSMofdifferencebyusingtheconvertedprobabilitymethod

圖5結果可以看出，其主采區(qū)和堆積區(qū)分布與最小檢出限法結果相同，摒棄的信息也主要分布在道路與廠房區(qū)域。但是，轉換概率法摒棄了更多的信息量，其空值區(qū)域范圍較大。

為了定量對比2種改進方法與原始DoD的計算結果，除了對比體積變化，還通過查詢首云鐵礦2014年8月—2016年10月間的開采報告，該地區(qū)鐵礦石平均密度大約為4.8 t/m3，平均含礦率約為10%，礦石平均品位約為30%。將2種方法所計算得到的總體積變化分別乘以平均密度、平均含礦率和礦石平均品位即可估算出動用儲量，計算結果如表2所示。

表2不同方法計算儲量變化結果對比

Tab.2Comparisonofresultsbydifferentmethods

方法堆積/m3開挖/m3總體積變化/m3儲量變化/t原始DoD1 325 479-14 633 968-13 308 489-1 916 422最小檢出限法(0.5 m)1 264 146-14 559 008-13 294 862-1 914 460最小檢出限法摒棄量61 333-74 960-13 627-1 962轉換概率法(95%)1 262 525-14 556 857-13 294 332-1 914 384轉換概率法摒棄量62 954-77 111-14 157-2 039

通過表2可以看到，儲量體積變化比較接近(13 294 332～13 294 862 m3)，并且都在一定程度上排除了誤差對于體積變化的干擾(13 627～14 157 m3)； 轉換概率法采用95%作為置信區(qū)間所摒棄的信息，要大于最小檢出限法以0.5 m為檢出限所摒棄的信息。查詢實際開采資料發(fā)現，2014年8月—2016年10月期間，首云鐵礦的實際開采量約為180萬t，本次研究2種方法的計算結果與礦山實際開采量都非常接近，其誤差約為6.4%。

4 結論

本文將無人機攝影測量的方法運用于露天礦的儲量動態(tài)監(jiān)測中，以獲取多時相、高空間分辨率的DSM為基礎，利用最新的地形變化算法快速計算出礦山的體積變化，進而估算出礦山動用儲量，為實現礦山科學管理提供了新的思路。主要結論如下：

1)通過2種地形變化算法計算出的動用儲量體積，通過乘以該礦山平均含礦率、礦石平均體重以及礦石平均品位，估算的該期間動用儲量與礦山實際開采量較為接近，誤差僅為6.4%，證實了該方法的實用性。

2)由于金屬礦產的品位、含礦率和密度具有一定的差異性，該方法的精度還有待進一步提高，目前更加適用于露天開采的煤礦或各類非金屬礦產。另外，由于無人機續(xù)航能力的限制，該方法比較適合中小范圍礦山的持續(xù)監(jiān)測。

3)無人機不僅可以獲得高空間分辨的DSM，還能方便快速地獲取高清的DOM，基于DOM與DSM，還能進行礦山環(huán)境監(jiān)測、災害巡查、礦山施工管理和邊坡設計與監(jiān)測等。無人機攝影測量技術的運用為礦山監(jiān)測和科學管理提供了全新的思路。

志謝： 感謝北京市國土局和密云區(qū)首云鐵礦的大力協(xié)助，以及中國冶金地質總局遙感院王乾工程師對DGPS測量的指導。