葉春陽 許傳華 孫國權 聶 聞,4
(1.中國科學院海西研究院泉州裝備制造研究所,福建泉州362000;2.中北大學電氣與控制工程學院,山西太原030000;3.中鋼集團馬鞍山礦山研究總院股份有限公司,安徽馬鞍山243000;4.江西理工大學資源與環境工程學院,江西贛州341000)
采礦區域的不斷擴張,伴隨而來的是裸露邊坡數量的增加,礦山邊坡一旦發生形變失穩將會給礦山生產以及人員生命安全帶來巨大的威脅[1-4]。因此通過采取一定的技術措施降低事故發生率對礦山生產安全具有重要意義[5-6]。邊坡形變識別方法按照監測形式可分為接觸式監測和非接觸式監測兩類。以接觸式監測為代表的邊坡形變識別方法主要是在邊坡中埋設各種傳感器設備,然后通過傳感器數據來識別邊坡形變區域。這種方式主要存在兩方面的不足:一是需要耗費大量的人力物力,且無法對設備安裝人員生命安全提供保障;二是獲取的信息有限,通常一個傳感器只能反映一個點的變化情況,因此這種方法獲取的數據難以反映邊坡整體變化情況。機器視覺是一種近幾年來發展較快且比較成熟的技術,它作為一種非接觸式的技術能夠有效減少人為參與并提升生產效率。通常機器視覺所用的監測儀器主要有深度相機、雷達等,許多學者利用這些設備提出了許多經典的邊坡數據處理方法。例如YANG等[7]利用超像素算法處理裁剪的邊坡視頻數據,獲取了邊坡形變區域的輪廓;MWANIKI等[8]利用圖像增強方法處理圖像數據提升邊坡的識別精度;OHNISHI等[9]研究了數字相機在邊坡形變位移上的應用,實現了用數字相機測量邊坡位移量;劉昌軍等[10]利用激光數據實現了高邊坡危巖體的識別;李洪梁等[11]利用三維激光掃描技術實現了堰塞湖的應急測繪。以上方法剔除了作業人員與邊坡的近距離接觸環節,在一定程度上避免了人員受到生命安全威脅;此外,相機和雷達獲取的數據范圍更廣、維度更高,從而更能反映邊坡的整體趨勢。
本研究通過室內模型試驗,首先利用深度相機實時獲取邊坡形變圖像數據和點云數據;然后分析相機坐標系、像素坐標系、世界坐標系、圖像坐標系之間的轉換關系,通過相機標定最終建立圖像與點云之間的聯系;最后利用背景差分法實現邊坡形變區域的識別,通過前后時刻的邊坡表面位移量估算邊坡形變的體積。本研究旨在探索邊坡形變二維數據與三維數據相結合實現邊坡連續形變區域的體積估算,為邊坡形變特征信息的自動提取提供參考。
本研究分別建立世界坐標系、相機坐標系(點云坐標系)、圖像坐標系以及像素坐標系,并對相機進行標定,確定點云坐標與像素坐標之間的參數關系。
世界坐標系用于描述相機的空間位置,真實反映相機在三維空間的坐標關系。如圖1所示,Ow、Xw、Yw、Zw分別代表世界坐標系的原點、X軸、Y軸和Z軸。

相機坐標系是聯系世界坐標系與圖像坐標系的過渡坐標系,是空間數據過渡到平面坐標之間的橋梁。從世界坐標系變換到相機坐標系屬于剛體變換,即物體不會發生形變,只需要進行平移和旋轉。通常以 Oc、Xc、Yc、Zc表示相機坐標系,本研究使用的點云數據所在坐標系即為相機坐標。
圖像坐標系是相機坐標系沿Z軸投影而來的,常用O、X、Y表示圖像坐標系的原點、X軸和Y軸,如圖2所示。

如圖2所示,坐標系原點O與相機坐標系原點Oc之間的距離為f,其中ΔABOc與ΔONOc相似,ΔPcBOc與ΔPNOc相似,則可得以下關系式:

進而求得相機坐標系與圖像坐標系的關系為

像素坐標系由圖像坐標系沿X軸負方向和Y軸負方向分別平移u0和v0個單位得到,常用O0、U、V來分別表示像素坐標系的原點、X軸和Y軸。如圖3所示,若矩形O0QWE為一幅圖像,則點(u0,v0)為該圖像的中心。

由圖3可得圖像坐標系與像素坐標系轉換關系:

式中dx,dy分別表示圖像坐標系下沿X軸和Y軸的長度系數。
由相機坐標系、圖像坐標系以及像素坐標系三者之間的轉換關系,可以推導出相機坐標系與像素坐標系有以下關系:

式中,a,b表示相機內參。
如圖4所示,選取15張不同拍攝角度下的標定板照片,根據張正友[12]標定法對相機標定。相機內參計算得:a=624.19,b=635.99。

光照不均會導致部分點云數據缺失。數據的缺失一方面會影響數據的完整性,另一方面可能導致計算結果與真實結果相差較大。插值技術能夠有效解決因數據缺失造成的誤差較大問題。因此本研究采用最近鄰插值算法[13-14]估計缺失數據。最近鄰插值算法的原理是根據距離缺失值最近的k個點的值估算缺失值,其定義如下:

式中,P為缺失值的估計值;pi為缺失值附近第i點的值。
本研究首先由式(4)確定點云數據中缺失值的位置,然后由式(5)對缺失值進行插值估算,其中k設置為5。
將坡體數據與背景數據分離有利于坡體中形變區域的識別。因此在進行邊坡形變區域識別前需要對數據進行預處理。首先在點云3個坐標軸方向設置閾值,獲取粗提取的邊坡數據;然后利用坡體與背景的顏色差異對坡體進行精細提取;最后保留的數據即為最終的坡體數據。依次處理每個時刻的邊坡點云數據以獲得完整時間序列下的坡體數據。
邊坡發生形變時會有兩個特征發生顯著性變化:一是邊坡形變區域顏色會發生改變;二是邊坡表面深度會發生改變。合理地利用邊坡形變時的兩個特征能夠有效地識別邊坡形變區域。本研究首先用顏色特征處理邊坡數據獲取帶有噪聲的邊坡形變區域;然后利用邊坡表面深度特征細化上述結果,獲取更加精細的邊坡形變區域。
在外界光源穩定的情況下,邊坡形變區域顏色特征不會發生變化或者發生微小變化,而形變區域顏色變化明顯。基于以上信息,本研究采用背景差分算法[15-17]剝離邊坡非形變區域數據,得到帶有噪聲的邊坡形變區域數據。
利用邊坡深度特征細化帶有噪聲的邊坡形變區域。隨著時間推移,邊坡表面某一區域并未發生形變,但是該區域在后續時刻的深度會發生細微變化。造成這種細微變化原因是邊坡表面土壤顆粒發生輕微移動。為防止邊坡表面土壤顆粒發生輕微移動的區域被誤識別為形變區域,設置一個閾值來避免該種情況發生。閾值選取規則為:選取初始時刻處理的點云數據,記該時刻點云平均深度為h0;選取視覺可見的坡體發生形變的第一張點云數據,記該時刻點云平均深度為h;變化量d=h-h0則為閾值。當邊坡點云變化量大于d,則視為發生形變;反之,則未發生形變。經過上述兩個步驟處理后的區域即為最終識別的邊坡形變區域。
為獲取邊坡形變數據,在室內進行模擬降雨導致邊坡形變破壞的試驗,并用深度相機對邊坡形變過程進行拍攝記錄。圖5為試驗現場圖,試驗模型箱長1.6 m、寬0.8 m、高0.8 m。試驗共分為兩組,A組降雨為65 mm/h,B組降雨為30 mm/h,兩組坡角均為45°。
全部程序由C++編程語言編寫。如圖6所示,第1列為原始點云,第2列為粗提取的邊坡形變區域結果,第3列為粗提取后細化的邊坡形變區域結果。從第1列數據可以看出邊坡發生形變后,坡體向內坍塌,此時該形變區域的Z坐標將發生變化。第2列數據顯示基于閾值分割后坡體上部未發生顯著形變的部分區域仍然被標記為形變區域,造成該種情況的原因是模擬的降雨影響了相機對該區域坐標的計算。噪聲剔除后,如第3列所示,可以看見優化后的區域只剩下坡底被標記的形變區域。


為了定量分析本研究方法的可靠性,以目標分割系數Dice[18]作為評價指標定量分析所提方法的性能。Dice是描述識別結果與原始破壞區域的重合一致性程度,其公式為

其中,G表示原始破壞區域,由人工分割獲取;P為本研究方法識別結果。Dice取值范圍為0~1。當Dice值越接近1時,表明識別結果越好;當Dice值越接近0時,識別結果越差。
由式(6)計算的圖6中A組、B組數據Dice值表明,A組數據Dice得分為91.24%,B組數據Dice分值為93.53%,兩者平均Dice為92.39%。計算結果充分證明本研究方法能夠精確提取邊坡形變區域。
通過上述方法對每組數據進行處理,獲取每個時刻標記的邊坡形變區域,然后利用識別的邊坡形變區域的坐標與初始時刻坐標進行作差后積分即可獲得隨時間變化的累計邊坡體積。如圖7所示,圖中兩條曲線為本研究方法求取的兩組試驗中邊坡累計體積隨時間的變化曲線。在相同的坡腳情況下,降雨大小為65 mm/h的A組試驗邊坡比降雨大小為30 mm/h的B組試驗邊坡更早產生形變。

針對邊坡連續形變信息難以精確提取的問題,提出了一種基于深度信息的邊坡形變區域識別方法。該方法基于識別的形變區域前后時刻深度信息的變化能夠準確估算邊坡形變的體積。由于該方法仍以試驗數據為主,因此采用真實數據對該方法進行進一步優化,是下一步的工作方向。