基于傾斜攝影測量的露天礦數值模型重構與計算實例

劉陽曉，朱萬成，劉文勝，劉溪鴿，王江梅

( 1. 東北大學 資源與土木工程學院，遼寧 沈陽 110819；2. 鞍鋼集團礦業(yè)有限公司，遼寧 鞍山 114001 )

數值模擬作為常用的露天礦穩(wěn)定性分析方法，其模型是否精準直接影響計算結果的準確性。傳統(tǒng)的數值計算模型，是由地質工作者現場實地采集有限的地形數據建立的。然而，由于大型露天礦山地形、地表構造較為復雜，與剝離開采、降水及植被生長等多因素關聯，以實際采集有限離散數據構建簡化地質模型的方法亟需改進。目前，隨著工業(yè)無人機技術的日益普及，傾斜攝影測量及配套三維實景建模軟件發(fā)展迅速，人們可以快速高效地獲得視覺良好的高精度露天礦三維實景模型。然而，如何進一步通過三維實景模型快速重構露天礦數值計算模型，為高效準確的數值建模提供新方法，更具現實意義。

傾斜攝影測量技術通過同一飛行平臺搭載多鏡頭，從垂直、傾斜多角度采集影像數據，利用三維實景建模軟件經過空中三角測量、影像密集匹配、構建TIN模型和紋理映射后[1]生成精度為厘米級[2-3]的三維實景模型。利用無人機靈活、安全、精度高、周期短的特點，傾斜攝影測量在露天礦山的地質調查[4-5]、資源監(jiān)管[6]、運輸監(jiān)測[7]、工程量計算[8-9]和邊坡變形監(jiān)測[10-11]方面均有研究和應用。針對通過傾斜攝影三維實景模型重構數值計算模型的問題，國內外學者也有一定的研究。高相波[12]等對具有地質災害隱患的危巖體進行高精度三維激光掃描獲得點云模型后，對邊坡發(fā)育的植被部分點云進行降噪、網格清理等人工處理后，構造邊坡NURBS曲面模型，實體化并劃分網格后可導入FLAC3D軟件；金愛兵[13]等選取植被覆蓋較少的邊坡進行攝影測量得到數字高程模型，通過Surfer、Rhino軟件轉化為網格模型，實現無人機測量到三維數值模型的快速建立；鄧小龍[14]等利用逆向工程軟件，由點云模型重構三維地質模型，用Hypermesh劃分網格后用于數值計算?？傮w而言，國內外學者的研究多集中在傾斜攝影測量精度、模型建立及露天礦監(jiān)測方面，針對傾斜攝影模型重構數值計算模型的研究，多選取沒有植被干擾的邊坡進行攝影測量，或對采集到的數據進行人工去噪后，利用軟件完成傾斜攝影模型到數值計算模型的轉換。

露天礦傾斜攝影測量，受植被、礦山設備和人員等的影響，建立的三維實景模型包含大量干擾 物，致使在進一步的模型重構與網格劃分中頻繁出錯，無法直接建立有限元數值計算模型。并且，對傾斜攝影測量模型人工去噪又極為費時費力，模型重構效率低下。針對上述問題，本文采用數字圖像處理的方法清除干擾物，生成規(guī)則化等高線，再構建實體模型并劃分網格，可高效完成由三維實景模型重構數值計算模型。

1 模型重構原理及網格劃分

本文數值計算模型重構方法分為兩部分，第1部分是編程實現的等高線提取及基于數字圖像處理的等高線規(guī)則化方法；第2部分是借助第三方軟件實現的由規(guī)則化等高線建立實體模型并劃分網格，生成網格模型用于數值計算。數值計算模型重構流程如圖1所示。

圖1 數值計算模型重構流程 Fig. 1 Reconstruction flow of numerical calculation model

1.1 等高線提取及規(guī)則化實現

1.1.1 模型預處理及等高線提取

由傾斜攝影采集數據建立的三維實景模型，多用于在場景中展示而無法導入其他三維建模軟 件。為解決三維實景模型難以二次開發(fā)的問題，本文通過遍歷模型三角網，提取三角網模型節(jié)點三維坐標構建點云模型。點云模型為已知坐標點的集合，模型可保存為文本，易于二次開發(fā)。

本文采用直通濾波的方法在點云模型中提取等高點云。直通濾波是對點云模型指定維度指定區(qū)間進行剪裁，保留區(qū)間內的點。將濾波維度設置在露天礦深度方向，將濾波區(qū)間設置在較小范圍，則直通濾波得到的點近似為等高點。

1.1.2 等高線規(guī)則化方法

受礦山情況復雜及實景模型存在離散三角面的影響，由點云模型提取的等高點存在噪聲點而無法直接導入三維建模軟件實現實體模型重構及網格劃分。本文將三維點云映射至二維圖像[15]，利用數字圖像處理方法剔除噪聲點，規(guī)則化等高線。采用點柵格化的方法將三維點云映射至二維圖像，以點云的點平均距離為像素點邊長，構造值為0的圖像矩陣。遍歷點云，利用式( 1 )確定三維點在二維圖像矩陣中的位置，并將矩陣值設置為1，即得到三維點云對應的二維圖像。

式中，round( )為四舍五入運算；(X，Y)為點云坐標； (I，J)為點在圖像矩陣中的坐標；L為圖像像素邊 長。

將三維點云映射至二維圖像后，通過填充圖像孔洞并統(tǒng)計連通區(qū)域面積，識別篩除噪聲點。剔除噪聲點后通過檢測孔洞邊緣獲得規(guī)則化后的等高點，將等高點按逆時針順序連接即得到規(guī)則化等高線。

( 1 ) 形態(tài)學開操作與閉操作。形態(tài)學操作利用結構元在圖像中遍歷，提取特定圖像分量。形態(tài) 學[16]的基本操作包括膨脹和腐蝕，膨脹操作使圖像變粗，腐蝕使圖像變細。膨脹和腐蝕的定義分別為

式中，A為圖像；B為結構元；為結構元相對原點的反射；z為結構元的原點；?為空集。

由式( 2 )可知，結構元B對圖像A膨脹的結果是結構元B的反射與圖像A有重疊時z的集合。結構元B對圖像A進行腐蝕的結果是結構元B中包含于圖像A中時z的集合。膨脹和腐蝕的組合使用構成開操作和閉操作，開操作和閉操作的定義為

開操作在圖像中的作用是斷開狹窄的間斷和消除細的突出物，閉操作則是為了填補輪廓線中的斷裂和消除小的孔洞。在本文研究中，將點云圖像填充孔洞前，用閉操作連接細的斷開部分，以避免填充溢出。在填充孔洞后，統(tǒng)計連通區(qū)域前，用開操作斷開細的連接以保證區(qū)分等高點云與噪聲點部分。利用形態(tài)學閉操作和開操作，使規(guī)則化等高點的程序更加健壯。

( 2 ) 填充孔洞及連通分量統(tǒng)計。填充孔洞在本文有2個作用，一是對原點云圖像進行邊界提取，會提取到2條邊緣線，填充孔洞后可以只提取到1條邊緣線；二是填充孔洞后可以根據連通區(qū)域的面積剔除離群點?？锥刺畛湟钥锥催吔鐬榧s束條件，以孔洞內1點為種子點，通過對孔洞循環(huán)膨脹操作直至圖像穩(wěn)定，孔洞全部填充?？锥刺畛涞墓綖?/p>

式中，Xk為填充孔洞后的圖像；為圖像補集；k＝0時，X0為初始種子點圖像。

對點云圖像填充孔洞后，由于噪聲點與等高點孔洞屬于不同連通分量，通過統(tǒng)計各個連通分量像素點個數，可以區(qū)分等高點和噪聲點，剔除噪聲 點。

( 3 ) Sobel邊緣檢測。邊緣檢測最常用的方法是檢測亮度的不連續(xù)性，其由梯度的幅值衡量。邊緣檢測的準則是在圖像中尋找梯度的幅值大于閾值的地方。梯度的幅值為

Sobel算子[17]采用水平和垂直2個方向的卷積模板與圖像進行卷積運算，以3×3鄰域的行和列之間的離散差計算梯度。圖像3×3鄰域示意圖如圖2( a )所示，Sobel算子水平和垂直梯度模板如圖2( b )，( c )所示。

圖2 3×3圖像鄰域及Sobel算子示意 Fig. 2 Neighborhood of image and Sobel operator

Sobel算子2個方向卷積模板與圖像進行卷積運算后，可得水平和垂直方向一階導數為

式中，Z為亮度。

由Sobel算子卷積得到水平和垂直梯度后，可得Z5處梯度公式為

如果?f大于等于閾值T，則該位置像素為邊緣像素。

( 4 ) 改進的區(qū)域生長算法連接無序點成線。傳統(tǒng)區(qū)域生長算法是基于灰度相似性的圖像分割方法，其基本思想是將具有相似性質的像素集合起來構成區(qū)域。具體地，確定種子點后，將在種子點像素周圍滿足相似規(guī)則的種子像素合并到區(qū)域中，直到沒有滿足條件的像素點可被包括，則1個區(qū)域生長完成，尋找下一個像素點進行生長[18]。本文在種子選取方法上進行改進，傳統(tǒng)的種子隨機生成，本文中種子的選取，是在1個區(qū)域生長的最后1個點按照8鄰域擴散的方式查找，由此可以保證2個鄰域首尾按順序連接。

1.2 實體模型構建及網格劃分

Rhino作為專業(yè)的三維建模軟件，可快速實現由線生成面及由面重構實體模型。本文由實體模型提取等高線后，通過對等高線的規(guī)則化處理，很大程度上避免了在Rhino軟件計算中的布爾運算錯誤和劃分網格時存在自相交面的錯誤。通過建模軟件Rhino可以由等高線生成露天礦表面模型，由表面模型邊框曲線擠出閉合曲面可快速建立三維實體模型。對實體模型利用Griddle插件劃分網格后，即生成可導入數值模擬的網格模型。

2 模型重構方法應用實例

2.1 三維實景模型建立

大孤山礦區(qū)總面積176.98 hm2，開采深度至-350 m。采用科比特M6無人機搭載5鏡頭從垂直方向及前后左右與地面45°方向進行一定重疊率的拍攝，采集得到圖像共31 110張及對應無人機POS( Positioning and Orientation System定位定姿系統(tǒng) )數據。外業(yè)獲得傾斜攝影的影像數據、POS信息和像控點坐標后，導入三維建模軟件ContextCapture進行空中三角測量，多視角影像密集匹配，構建TIN模型和3D紋理映射后，生成高精度的三維場景模型，如圖3所示。

圖3 大孤山露天礦三維實景模型 Fig. 3 Real scene terrain model of Dagushan open-pit mine

2.2 等高線提取及規(guī)則化處理

2.2.1 大孤山模型預處理及等高點提取

通過分析三維實景模型的數據結構，利用三維渲染引擎OSG庫的節(jié)點訪問器NodeVisitor遍歷實景模型的三角面節(jié)點，提取節(jié)點三維坐標構建點云模型。大孤山點云模型點數量8 000余萬，為減少運算成本，提高計算效率，采用體素化網格的方法對模型下采樣，在保持點云形態(tài)不變的情況下，將模型點數量由8 000余萬減少至532余萬，下采樣后，大孤山露天礦點云模型如圖4所示。

圖4 大孤山露天礦點云模型 Fig. 4 Point cloud model of Dagushan open-pit mine

大孤山露天礦深度方向即為Z軸方向，最小值為-351.92 m，最大值為50 m。點平均間距為 0.95 m。以Z軸為濾波字段，濾波區(qū)間設置1.5 m，相鄰濾波區(qū)間設置為3 m，對模型直通濾波處理后獲得等高點云模型，如圖5所示。

圖5 等高點云示意 Fig. 5 Contour line of point cloud model

2.2.2 等高點規(guī)則化

受點云模型中冗余點與離群點的影響，未規(guī)則化的等高線在Rhino?軟件中構建實體模型出現布爾運算錯誤及網格劃分存在自相交面錯誤，無法完成網格模型重構。本文通過點柵格化的方法，將三維點云映射至二維圖像后，通過數字圖像處理的方法規(guī)則化等高線。

在Beijing1954坐標系下，大孤山點云模型X范圍[-608.36，1 192.35]，Y范圍[-792.4，858.43]。以大孤山模型-269 m高程等高點( 圖6( a ) )為例，構建大小為1 892×1 735的全0圖像矩陣，通過式( 1 )，令L＝0.95，Xmin＝-608.36，Ymin＝-792.4，將三維點云映射至二維圖像，如圖6( b )所示。對于等高點圖像，填充孔洞后( 圖6( c ) )，規(guī)則化等高點與噪聲點區(qū)域不連 通，屬于不同連通分量。通過統(tǒng)計圖像中連通區(qū)域分布，生成連通分量標記矩陣( 圖6( d ) )來區(qū)分噪聲 點。

-269 m點云圖像，連通分量標記結果如圖6( d )所示。點云圖像共包含8個連通分量，等高點構成的連通分量( 圖6( d )粉色區(qū)域 )面積最大，7個噪聲點構成的連通分量( 圖6( d )紅、藍、黃、綠、粉紅、紫、青區(qū)域 )面積小。統(tǒng)計各連通分量像素點個數，區(qū)分并剔除噪聲點，得到無噪聲點云圖像( 圖6( e ) )。經過邊緣檢測，提取邊緣點即得到規(guī)則化等高點圖像( 圖 6( f ) )。為檢驗規(guī)則化等高線的精度，將規(guī)則化等高線與等高點云對比( 圖6( g ) )可以看出，規(guī)則化等高線保持原點云模型形態(tài)。將無序的等高點規(guī)則化處理后，通過改進的區(qū)域生長算法，將點按逆時針方向連接為線即為規(guī)則化等高線( 圖7 )。

圖6 -269 m等高線規(guī)則化處理流程 Fig. 6 Flow diagram of regulating contour line at -269 m

圖7 大孤山露天礦規(guī)則化等高線 Fig. 7 Regularized contour line of Dagushan open-pit mine

2.3 實體模型建立重構及網格劃分

規(guī)則化的等高線，利用Rhino?軟件，可快速由線生成礦山殼體模型，拉伸生成封閉模型后，利用Griddle插件劃分網格后，生成FLAC3D格式的網格模型，如圖8所示。模型中面體單元數為687 052，節(jié)點數為132 713。

圖8 大孤山露天礦三維網格模型 Fig. 8 Three dimensional grid model of Dagushan open-pit mine

3 模型精度驗證及數值計算實例

3.1 模型精度驗證

完成數值計算模型的重構后，對模型精度 進行分析。以三維實景模型為準確地形模型參 照，對比本文重構網格模型與傳統(tǒng)方法建立的網 格模型精度。分別對三維實景模型、重構網格模 型及傳統(tǒng)網格模型在X與Y方向的1/4，1/2，3/4處 取6條剖線，如圖9( a )所示，剖線5的對比如圖9( b )所示。

圖9 大孤山露天礦代表性剖面選取與對比 Fig. 9 Selection and comparison of representative sections in Dagushan open-pit mine

每個模型獲得6個剖面的高程向量，通過計算網格模型與實體模型剖面高程向量的歐式距離，判斷剖面相似性。歐式距離定義為

式中，z1，z2為剖面高程向量；dist(z1，z2)為2向量的歐式距離。

歐式距離越小，代表剖面相似度越高，距離大則剖面相似度低。令z1為三維實景模型剖面高程向量，令z2分別為本文重構網格模型和傳統(tǒng)網格模型的高程向量，計算模型剖面間歐式距離( 表1 )。由表1可知，本文重構模型與實景模型剖面歐式距離小于傳統(tǒng)網格模型與實景模型剖面的歐式距離，即本文重構模型與三維實景模型剖面更加相似。這是因為，傳統(tǒng)數值計算模型的建立由于人工采集地形數據有限而被大大簡化，與實際礦山地形不符。本文基于無人機采集的大量礦山地形數據重構網格模型與實際礦山地形更加符合，準確度高于傳統(tǒng)數值計算模型。

表1 模型剖面歐式距離 Table 1 Euclidean distance of profiles

3.2 數值計算

為驗證重構的網格模型在數值模擬中的三維應用，將模型導入FLAC3D進行力學分析計算。計算模型大小為2 000 m×2 000 m×890 m。數值模擬采用線彈性本構模型，根據對露天礦取巖樣后進行室內物理力學試驗結果，采用計算模型參數為：巖體體積模量(K)為4.44 GPa；剪切模量(G)為3.33 GPa；密度(ρ)為2 300 kg/m3；重力加速度為9.8 m/s2。在模型前后邊界施加沿X方向固定水平位移約束，左右邊界施加沿Y方向的固定水平位移約束，底部施加沿Z方向的固定垂直位移約束。通過迭代計算，當體系最大不平衡力與典型內力的比率小于給定值( 1.0×10-5)時計算收斂，得到在重力作用下生成的初始地應力場。其中，垂直地應力分量計算結果如圖10所示。

圖10 大孤山露天礦垂直地應力分布規(guī)律 Fig. 10 Distribution of vertical stress in Dagushan open-pit mine

沿該露天礦最深標高處的點以下、巖體內部一點P為研究對象，其中，此處P恰取為最深標高點與整體模型底部距離的中點。通過數值計算后處理，得到該點的垂直地應力分量σv＝5.0 MPa；而通過巖體容重( 或密度 )進行地應力理論計算：σv＝ρgh＝2 300×9.8×222×10-6＝5.004 MPa；二者之間的相對誤差不超過0.8%。

4 結 論

( 1 ) 提出了針對三維實景模型的等高線提取及規(guī)則化方法。采用直通濾波的方法提取等高點，通過數字圖像處理的方法剔除噪聲點及冗余點，生成規(guī)則化等高線。

( 2 ) 借助第三方軟件，由規(guī)則化等高線可快速實現實體模型重構及網格劃分，生成數值計算模型，提高數值建模效率。

( 3 ) 通過與實測的三維地質模型進行剖面的對比分析，發(fā)現本文方法相比傳統(tǒng)數值建模，模型更加符合露天礦實際地形。