史雨川，李 浩，楊 彪，劉慶群

（1. 河海大學(xué) 地球科學(xué)與工程學(xué)院，江蘇 南京 210098）

基于普通數(shù)碼影像的露采礦山三維可視化分析

史雨川1，李 浩1，楊 彪1，劉慶群1

（1. 河海大學(xué) 地球科學(xué)與工程學(xué)院，江蘇 南京 210098）

研究了基于普通數(shù)碼影像的礦山三維可視化技術(shù)。首先將檢校后的數(shù)碼相機(jī)以物方控制的單模型攝影方式對(duì)目標(biāo)拍攝，取得原始影像并逐一校正；然后自動(dòng)進(jìn)行像對(duì)定向、立體匹配，采集影像點(diǎn)云；最后根據(jù)影像匹配點(diǎn)云三維建模，并依據(jù)前期地形數(shù)據(jù)、終采境界設(shè)計(jì)圖、地質(zhì)勘察數(shù)據(jù)對(duì)露采礦山進(jìn)行三維可視化分析，實(shí)現(xiàn)動(dòng)態(tài)檢測(cè)。

露采礦山；普通數(shù)碼影像；三維可視化；近景攝影測(cè)量

本文采用近景攝影測(cè)量方法，基于普通數(shù)碼影像研究露采礦山三維可視化技術(shù)，通過(guò)量測(cè)化檢校的相機(jī)對(duì)礦山開(kāi)采宕口進(jìn)行拍攝，并對(duì)原始影像畸變校正，自動(dòng)進(jìn)行相對(duì)定向、絕對(duì)定向以及影像立體匹配，利用采集的匹配點(diǎn)云重建礦山三維模型，并根據(jù)前期地形數(shù)據(jù)、終采境界設(shè)計(jì)圖進(jìn)行可視化分析[1-3]，最終實(shí)現(xiàn)礦山影像數(shù)據(jù)采集、處理、三維重建和可視化的一體化實(shí)時(shí)動(dòng)態(tài)檢測(cè)。

1 數(shù)碼相機(jī)量測(cè)化檢校

由于普通數(shù)碼相機(jī)不提供畸變系數(shù)和內(nèi)方位元素，因此必須進(jìn)行量測(cè)化檢校，使之獲取的數(shù)字影像具有可控的精度。為求得畸變系數(shù)k1、k2、k3、p1、p2，應(yīng)用精密平面控制場(chǎng)，將構(gòu)像畸變差改正模型納入影像透視變換，并用未知數(shù)交替解算[4]。再根據(jù)三維控制場(chǎng)的物方坐標(biāo)及經(jīng)過(guò)畸變改正的控制點(diǎn)像片量測(cè)坐標(biāo)，利用單張像片后方交會(huì)算法解算內(nèi)方位元素f、x0、y0[5,6]。將解算后的內(nèi)方位元素重新代入上述畸變系數(shù)求解過(guò)程，利用新的畸變系數(shù)繼續(xù)解算內(nèi)方位元素，如此交替迭代直到畸變系數(shù)、內(nèi)方位元素變化較小時(shí)停止，流程如圖1。

實(shí)驗(yàn)顯示，普通數(shù)碼影像畸變最大可達(dá)幾十個(gè)像元，校正后中誤差小于1個(gè)像元，且主距固定時(shí)畸變系數(shù)和內(nèi)方位元素比較穩(wěn)定。

2 礦山數(shù)據(jù)快速采集與三維建模

2.1 影像快速采集

圖1 畸變系數(shù)和內(nèi)方位元素檢校交替計(jì)算流程

從露采礦山需求考慮，由于單幅影像已經(jīng)可以覆蓋目標(biāo)，所以采用外業(yè)工作量較小的物方控制的單模型攝影方式。拍攝前，在礦山宕口布置若干個(gè)控制點(diǎn)，精確測(cè)定控制點(diǎn)的坐標(biāo)，通過(guò)單張像片空間后方交會(huì)計(jì)算影像的外方位元素。采取立體拍攝的辦法分別在左站和右站對(duì)區(qū)域目標(biāo)進(jìn)行外業(yè)拍攝，且交向角不宜太大，基線長(zhǎng)度與拍攝距離比為1/10左右，既能保證精度，又便于影像處理，相鄰像對(duì)的拍攝覆蓋范圍保證銜接即可，如圖2。

圖2 礦山攝影方法示意圖

為滿足像對(duì)絕對(duì)定向中模型變形改正要求，每個(gè)像對(duì)至少要有4個(gè)控制點(diǎn)，且均勻分布在像片四角，相鄰像對(duì)共用2個(gè)像控點(diǎn)?？刂泣c(diǎn)可以采用標(biāo)志點(diǎn)布設(shè)在礦山上，也可以選取特征明顯的點(diǎn)。使用免棱鏡全站儀以自由設(shè)站方式精確測(cè)量其空間坐標(biāo)，利用GPS測(cè)量其中任意2個(gè)以上控制點(diǎn)的大地坐標(biāo)，使其在數(shù)據(jù)處理時(shí)控制點(diǎn)坐標(biāo)能夠轉(zhuǎn)換為大地坐標(biāo)系。礦

山單像對(duì)拍攝原始影像如圖3所示。

圖3 礦山原始影像

根據(jù)測(cè)定得到的相機(jī)畸變系數(shù)和內(nèi)方位元素，利用畸變差改正模型，對(duì)原始數(shù)碼影像像素重排，得到畸變校正后的影像。

2.2 像對(duì)自動(dòng)定向參數(shù)解算

傳統(tǒng)的相對(duì)定向元素求解是采用迭代法，優(yōu)點(diǎn)是只需較少的相對(duì)定向點(diǎn)，且定向精度高；缺點(diǎn)是必須提供較準(zhǔn)確的未知數(shù)初值，一旦初值設(shè)置不當(dāng)，求解方程迭代就可能不收斂。因露采礦山采用徒手拍攝，條件多變且目標(biāo)復(fù)雜，所以相對(duì)定向元素初值很難估算[7,8]。為此，利用控制點(diǎn)像片坐標(biāo)和空間坐標(biāo)，通過(guò)角錐體法單像空間后方交會(huì)計(jì)算得到左右像片外方位元素，分別設(shè)為和并依據(jù)外方位元素自動(dòng)設(shè)定相對(duì)定向初值。在攝影測(cè)量坐標(biāo)系中，左右片旋轉(zhuǎn)矩陣為：

設(shè)R1、R2分別為左右片的旋轉(zhuǎn)矩陣，則右片相對(duì)于左片的旋轉(zhuǎn)矩陣為：

經(jīng)解算得出連續(xù)法相對(duì)定向初值：

初值確定之后，還需要若干對(duì)同名像點(diǎn)。利用Forstner算子自動(dòng)提取左右片中的特征點(diǎn)，利用相對(duì)定向初值建立的概略核線作為約束條件，然后進(jìn)行特征點(diǎn)匹配，自動(dòng)剔除小區(qū)域內(nèi)視差異常的錯(cuò)誤匹配點(diǎn)，得出正確的同名點(diǎn)。將同名像點(diǎn)坐標(biāo)與初值代入迭代公式，自動(dòng)進(jìn)行相對(duì)定向元素解算。

為確定相對(duì)定向所得幾何模型在物方坐標(biāo)系中的絕對(duì)位置，需進(jìn)行空間相似變換，即像對(duì)絕對(duì)定向?？臻g相似變換是7個(gè)絕對(duì)定向元素的非線性函數(shù)，仍采用迭代法進(jìn)行解算。由于連續(xù)法相對(duì)定向中采用左片的像空間坐標(biāo)系為基準(zhǔn)，所以模型空間方位由左片外方位元素確定，故絕對(duì)定向元素初值為取2個(gè)控制點(diǎn)的實(shí)際邊長(zhǎng)和對(duì)應(yīng)模型長(zhǎng)度比值作為模型比例系數(shù)的初值，絕對(duì)定向點(diǎn)自動(dòng)選取之前在計(jì)算影像外方位元素時(shí)已經(jīng)人工選取的實(shí)測(cè)控制點(diǎn)，最后自動(dòng)進(jìn)行絕對(duì)定向元素的迭代解算。

2.3 礦山影像立體匹配

由于露采礦山攝影距離較近，且開(kāi)采表面崎嶇變化較大，所以拍攝位置的變化容易引起較大的目標(biāo)構(gòu)像變形。其中包括攝影比例尺的差異、影像斷裂、同名子影像形狀變化和位移，而且不同視點(diǎn)會(huì)導(dǎo)致目標(biāo)的背景不同，因此常規(guī)的單基線立體匹配的可靠性和精度較低[9]。為提高匹配準(zhǔn)確性，設(shè)計(jì)了一種基于種子點(diǎn)約束的金字塔多級(jí)概率松弛匹配方法。

獲取種子點(diǎn)的目的是用其約束搜索范圍以提高匹配精度及降低匹配代價(jià)。方法是對(duì)于待匹配的格網(wǎng)點(diǎn)根據(jù)其周?chē)姆N子點(diǎn)估計(jì)其右片的對(duì)應(yīng)位置，將左右影像的種子點(diǎn)對(duì)看作匹配的線段，種子點(diǎn)對(duì)之間的格網(wǎng)點(diǎn)根據(jù)右片相匹配的線段及格網(wǎng)點(diǎn)個(gè)數(shù)均分估計(jì)待匹配點(diǎn)的初始位置，并由前一點(diǎn)的相關(guān)系數(shù)修改估計(jì)的初始值位置，以適應(yīng)視差的變化。依次遍歷每個(gè)格網(wǎng)點(diǎn)，得到每個(gè)格網(wǎng)點(diǎn)的候選匹配點(diǎn)[10,11]。流程如圖4，匹配結(jié)果如圖5。

圖4 影像匹配流程圖

圖5 匹配效果圖

圖6 紋理映射后的單像對(duì)模型

2.4 坐標(biāo)轉(zhuǎn)換與三維建模

匹配完成后，即根據(jù)相對(duì)定向、絕對(duì)定向參數(shù)和匹配結(jié)果生成影像點(diǎn)云，對(duì)目標(biāo)進(jìn)行三維重建。考慮到控制點(diǎn)坐標(biāo)是通過(guò)全站儀自由設(shè)站獲得的，因此需要將采集點(diǎn)坐標(biāo)從自由坐標(biāo)系轉(zhuǎn)換到大地坐標(biāo)系，轉(zhuǎn)換時(shí)只需考慮水平方向的旋轉(zhuǎn)和原點(diǎn)平移。根據(jù)實(shí)測(cè)的2個(gè)以上控制點(diǎn)的大地坐標(biāo)、自由坐標(biāo)，計(jì)算坐標(biāo)系旋轉(zhuǎn)參數(shù)和平移參數(shù)，然后根據(jù)式（5）將采集點(diǎn)坐標(biāo)（X，Y，Z）轉(zhuǎn)換到大地坐標(biāo)（Xt，Yt，Zt），通過(guò)插值生成DEM。紋理映射后的單像對(duì)模型如圖6所示，礦山開(kāi)采宕口三維模型如圖7所示。

圖7 分級(jí)渲染后的礦山開(kāi)采宕口三維模型

3 礦山三維可視化分析

3.1 精度評(píng)價(jià)

為評(píng)價(jià)精度，精確測(cè)量布置的8個(gè)檢查點(diǎn)坐標(biāo)，將攝影測(cè)量計(jì)算的空間坐標(biāo)與外業(yè)實(shí)測(cè)坐標(biāo)進(jìn)行對(duì)比統(tǒng)計(jì)。使用Canon G5數(shù)碼相機(jī)對(duì)南京浦口某礦山進(jìn)行拍攝，拍攝距離約100 m，對(duì)比結(jié)果如表1所示。

表1 檢查點(diǎn)精度/m

統(tǒng)計(jì)檢查點(diǎn)坐標(biāo)中誤差為：mx=0．080 m，my=0．059 m，mz=0．075 m。實(shí)驗(yàn)表明，使用本文方法采集點(diǎn)云坐標(biāo)相對(duì)精度高于1‰，滿足一般礦山精度要求。

3.2 儲(chǔ)量分析

礦山三維模型底面積計(jì)算根據(jù)終采邊界拐點(diǎn)的XY平面坐標(biāo)計(jì)算得到，即計(jì)算終采邊界拐點(diǎn)所構(gòu)成平面區(qū)域的面積，亦為終采范圍內(nèi)模型在水平面上的投影面積。礦山三維模型表面積計(jì)算是基于DEM進(jìn)行的，將DEM每個(gè)格網(wǎng)分成兩個(gè)對(duì)角三角形，遍歷每個(gè)三角形，其面積之和即為格網(wǎng)的表面積，亦為模型的表面積。體積計(jì)算時(shí)，將礦山三維模型分解成若干個(gè)三棱柱體，計(jì)算每個(gè)三棱柱體在最低開(kāi)采高程以上的體積，然后對(duì)這若干個(gè)體積累加求和，即得到模型體積。

儲(chǔ)量計(jì)算是基于高精度礦山三維模型自動(dòng)進(jìn)行的，因此依次建立以下礦山DEM：①將礦山開(kāi)采前的原始柵格地形圖矢量化，插值得到礦區(qū)原始DEM0；②將設(shè)計(jì)好的礦區(qū)終采境界圖數(shù)字化，生成終采境界DEM1；③根據(jù)數(shù)碼影像采集的點(diǎn)云數(shù)據(jù)生成礦山開(kāi)采宕口DEM2，獲取目前已開(kāi)采區(qū)域R；④根據(jù)平面位置將DEM0與DEM2疊加，用區(qū)域R內(nèi)的DEM0數(shù)據(jù)替換為DEM對(duì)應(yīng)區(qū)域的數(shù)據(jù)，得到礦區(qū)當(dāng)前DEM[12]。

23

儲(chǔ)量計(jì)算具體步驟如下：①分別計(jì)算礦界內(nèi)DEM0、DEM1、DEM3的體積V0、V1、V3；②計(jì)算總浮土層表面積S0，即最低開(kāi)采高程以上DEM0的表面積；③計(jì)算已剝離浮土層表面積S1，即最低開(kāi)采高程以上、開(kāi)采區(qū)域R內(nèi)的DEM0的表面積；④總儲(chǔ)量T0=（V0-V1-S0×h）×a×t，已采量T1=（V0-V3-S1×h）×a×t，剩余儲(chǔ)量T2=T0-T1。其中，a為含礦率，t為礦石比重，h為浮土層厚度。利用每期采集的礦山開(kāi)采宕口DEM不斷更新，可以清晰地反映出礦山儲(chǔ)量動(dòng)態(tài)，精確分析超開(kāi)采情況。

以南京浦口某礦山為例，對(duì)開(kāi)挖宕口共拍攝5個(gè)像對(duì)，使用免棱鏡全站儀測(cè)量23個(gè)控制點(diǎn)，并通過(guò)GPS測(cè)量其中2個(gè)點(diǎn)的大地坐標(biāo)作為轉(zhuǎn)換點(diǎn)，外業(yè)影像獲取及數(shù)據(jù)采集約45 min，內(nèi)業(yè)數(shù)據(jù)處理耗時(shí)1 h左右。礦區(qū)面積102 900 m2，礦石比重2．8 t/m3，含礦率78%，浮土層厚度2 m。儲(chǔ)量計(jì)算結(jié)果如表2所示，對(duì)比外業(yè)實(shí)測(cè)計(jì)算已采儲(chǔ)量2 123 806 t，相對(duì)精度達(dá)到1．5‰，完全滿足礦山要求。

表2 儲(chǔ)量計(jì)算結(jié)果/t

[1] 朱煜峰． 概述我國(guó)礦山測(cè)量技術(shù)的新進(jìn)展[J]． 中國(guó)礦業(yè), 2004, 13(4): 7-8

[2] 段奇三． 徠卡HDS 8800三維激光掃描儀在露天礦中的應(yīng)用[J]．測(cè)繪通報(bào), 2011(12): 79-80

[3] 楊彪, 嚴(yán)麗英, 李浩． 基于普通數(shù)碼相機(jī)的露采礦儲(chǔ)量動(dòng)態(tài)檢測(cè)系統(tǒng)[J]． 工程勘察, 2013(3): 54-59

[4] 陳舒, 李浩, 黃河． 普通數(shù)碼相機(jī)三種檢校方法的量測(cè)精度評(píng)價(jià)[J]． 遙感信息, 2012, 27(5): 73-76

[5] Lucchese L, Mitra S K． Using Saddle Points for Subpixel Feature Dectecyion in Camerapore, a Calibration Targets[C]． 2002 Asia Pacific Conference on Circutis and Systems, Singapore, 2002

[6] Lucchese L． Geometric Calibration of Digital Cameras Through Multi-view Rectification[J]． Image and Vision Computing, 2005, 23(5):517-539

[7] 楊立君, 馬明棟, 苗立志． 非量測(cè)相機(jī)近景數(shù)字影像相對(duì)定向方法研究[J]． 大地測(cè)量與地球動(dòng)力學(xué), 2012, 32(4): 135-138

[8] Tsai R Y, A Versatile Cemera Calibration Technique for Highaccurancy 3D Machine Vision Metrology Using Off-theshelf Cameras and Lenses[J]． IEEE Journal of Robotics and Automation, 1987, 3(4): 323-344

[9] 陳新璽． 多基線普通數(shù)碼影像近景攝影測(cè)量技術(shù)研究[D]． 南京: 河海大學(xué), 2006

[10] Zhang L． Automatic Digital Surface Model(DSM) Generation from Linear Array Images[D]． Zurich: Swiss Federal Institute of Technology, 2005

[11] Lowe D G． Distinctive Image Features from Scale-Invariant Interest Points[J]． International Journal of Computer Vision, 2004, 60(2): 91-110

[12] 蔣銳, 宋煥斌． 基于三維柵格數(shù)據(jù)的露采礦山儲(chǔ)量動(dòng)態(tài)監(jiān)測(cè)研究與應(yīng)用[J]． 礦產(chǎn)與地質(zhì), 2009, 23(5): 469-472

P234.1

B

1672-4623（2014）04-0017-04

10.11709/j.issn.1672-4623.2014.04.006

史雨川，碩士，研究方向?yàn)閿?shù)字?jǐn)z影測(cè)量。

2013-08-27。

項(xiàng)目來(lái)源：國(guó)家自然科學(xué)基金資助項(xiàng)目（51079053）。