張慶斌，趙之星，劉志奇

(中晉環(huán)境科技有限公司，山西 太原 030002)

0 引 言

滑坡監(jiān)測歷來是地質(zhì)災害研究領(lǐng)域中的重點，大量的開采、施工或長期惡劣氣候的影響，使得原始地貌形成了多種不穩(wěn)定邊坡，這些邊坡由于受力失衡會產(chǎn)生定向移動[1]。大面積滑坡會對礦區(qū)工作人員、附近居民以及其他相關(guān)作業(yè)人員的生命和財產(chǎn)造成嚴重危害。其中突發(fā)性的滑坡現(xiàn)象危險性最大，其隱蔽不易辨別的征兆和災害發(fā)生時無法應對的緊迫性，是迄今為止災害監(jiān)測預警的難題[2]。為此大量學者和機構(gòu)開發(fā)了各種監(jiān)測手段和預警系統(tǒng)。上海華測導航科技有限公司推出了基于綜合技術(shù)的自動實時監(jiān)測系統(tǒng)，該系統(tǒng)應用了GNSS接收機、固定式測斜儀、裂縫計、雨量計等傳感器，以及無線傳輸裝置和太陽能電池等輔助設(shè)備[3-6]。該系統(tǒng)利用多種手段獲取信息，綜合評估目標變化，得出結(jié)果，但是其設(shè)備價格昂貴，一旦發(fā)生滑坡，系統(tǒng)設(shè)備全部損壞，監(jiān)測成本較高。成都理工大學團隊[7-9]利用機載雷達和傾斜相機對災害體進行采集，通過內(nèi)業(yè)處理生成三維激光點云實景三維模型，并對多期數(shù)據(jù)進行比較，該方法能直觀顯示出測區(qū)隱患點，為決策提供準確的依據(jù)，且對植被茂密災害點隱蔽的測區(qū)尤為勝任。但該方法處理過程較長，特別是三維模型的生產(chǎn)占用大量的時間，因此監(jiān)測實時性略顯滯后。文獻[10-12]利用近景攝影測量、計算機視覺等手段，通過固定雙相機，利用影像匹配技術(shù)，結(jié)合內(nèi)外方位元素計算目標的三維坐標，并分析多期的坐標變化以達到目標位移的計算。該方法不需要在災害體上安裝設(shè)備，作業(yè)成本低[13-15]。但為了得出空間坐標，需要進行較為復雜的計算，且受制于攝影測量本身精度的限制，對突發(fā)性應急存在一定的不適應[16-18]。針對滑坡的突發(fā)性和應急性，通過固定相機對滑坡體拍照，對相片的像素進行計算，設(shè)定多層篩選條件和閾值，令計算機準確分析出實際滑坡和誤差影響或分散式局部移動的情況，做出相應決策。該方法延續(xù)了近景攝影非接觸式的優(yōu)勢，同時舍棄求取三維坐標的繁雜計算過程，僅通過影像進行判別，響應速度快，自動化程度高，對滑坡災害的應急較為適用[19-20]。

1 滑坡預警算法推導

1.1 初步判斷標準

固定相機按照一定時間間隔對滑坡體進行連續(xù)拍照。相機固定即像片不動，像片不動即像平面坐標系不變，當滑坡目標不發(fā)生移動時，理論上兩幀影像提取的同名點像素坐標應相等,即

受到噪聲影響，兩幀影像同名點坐標會產(chǎn)生一定的偏差，由于噪聲服從各向同性高斯分布，使得總體上點坐標偏差方向各異，偏差均值為0，即

(1)

事實上，除了隨機噪聲影響外，特征提取中存在粗差以及樣本數(shù)量有限等影響因素，導致上述坐標偏差的均值與理想值0存在偏離，因此需要定義置信閾值，當該偏離在閾值范圍內(nèi)，則同名點坐標的變化僅受到噪聲的干擾，沒有發(fā)生實際的移動。結(jié)合坐標誤差的理論分布，令

式中：Δui、Δvi分別為第k張影像與第1張影像的同名點橫、縱坐標差。

則樣本方差的估值為

(2)

式中：S為標準差。

構(gòu)建樣本樞軸量

(3)

圖1 誤差服從t分布的α分位Fig.1 Alpha quantiles with errors following t distribution

反之，該情況下特征點坐標偏差的均值和理論期望之間的差異是顯著的，因此認為像點受到了噪聲以外的干擾，需進行進一步分析。

1.2 位移情況甄別

若像點的偏差均值過大，則需深入分析，以判斷變化是由粗差影響還是本身發(fā)生移動造成的，具體方法如下。

計算特征點的梯度方向：

(4)

式中：θi為第1張影像相較于第k張影像同名點發(fā)生位移的梯度方向。

繪制特征點梯度方向的直方圖，橫軸為梯度方向的角度區(qū)間，在45°～135°內(nèi)，每間隔5°為1個區(qū)間，共18個，縱軸為滿足每個區(qū)間點的數(shù)量，即N1、N2、N3，…，N18，將點數(shù)量按大小排列，并提取出最多的2組進行分析。假設(shè)點數(shù)量最大的一組為Ni，將Ni個像點按照橫坐標從小到大依次計算相鄰點的距離

式中：di為同一影像上相鄰像點間的距離。

以第1個點與第2個點的相鄰距離d1為基準，定義相鄰間距的權(quán)重pi為

pi=d1/di

(5)

式中：mi為第1張影像相較于第k張影像第i個同名點發(fā)生的位移。

(6)

同時根據(jù)相似關(guān)系得出一個像素對應物方的實際距離為

GSD=f/(Lδ)

(7)

式中：GSD為1個像素對應的實際距離；f為焦距；L為拍攝距離；δ為單個像素大小。

因此發(fā)生移動的距離為D=GSD(u2+v2)1/2，其中，u、v為目標點像素坐標。

1.3 決策階段

根據(jù)以上分析，總結(jié)本次影像應急監(jiān)測預警的條件，如圖2所示。

圖2 滑坡應急預警決策流程

2 滑坡預警模擬試驗

2.1 試驗布置和數(shù)據(jù)采集

由于數(shù)量稀少，危險性等原因，對理想突發(fā)性滑坡測區(qū)的選定較為困難，考慮時間成本等因素，選擇進行模擬試驗。試驗區(qū)域選擇太原市中國冶金地質(zhì)總局新辦公樓施工區(qū)西北角的人工邊坡，邊坡坡度正切值約0.7(坡面角35°)，土層15 cm以下布置1.2 cm×1.2 cm的襯墊，以便模擬滑坡移動過程。為盡可能接近滑坡體結(jié)構(gòu)，在邊坡上摻雜了碎石組成的松散堆積體。采用SONY_A7Ⅱ相機，將其固定在改裝后的三腳架上，在距目標3 m處粗略整平，并保持不動。首先在坡面未動的情況下拍攝首張照片(圖3a)，然后分3種情況拍攝后續(xù)照片：保持坡面不動(圖3b)，人為使坡面上某點或較為分散的零星點位發(fā)生移動(圖3c)，人工拉拽襯墊，坡面滑移(圖3d)。分別將圖3b、圖3c、圖3d與圖3a進行對比分析，驗證方法的正確性。

2.2 數(shù)據(jù)處理

將圖3b中影像首先通過優(yōu)化后的SIFT算法進行同名點提取與匹配，利用基于本方法開發(fā)的程序進行處理分析，計算每組像片間同名點的像素坐標變化情況，驗證本方法對3種位移情況的判斷及響應效果。

2.2.1 未移動影像組處理分析

第1對影像選擇初始影像和滑坡未動影像進行軟件分析，通過同名點匹配、統(tǒng)計匹配點樣本數(shù)量，計算同名點像素坐標差值、平均值及標準差估計值，依據(jù)假設(shè)檢驗求取檢驗統(tǒng)計量的絕對值，并與其服從的分布拒絕閾的臨界點進行比較。處理結(jié)果如圖4所示。共提取特征點1 858個，匹配成功的同名點1 553個，紅色箭頭為像點移動的方向，經(jīng)比較計算得出檢驗統(tǒng)計量為1.78，α取0.05時，t分布的拒絕域臨界點為1.96，顯然統(tǒng)計量小于臨界點。

2.2.2 局部分散式移動影像組處理分析

第2對影像為初始影像與分散局部移動個別點位的影像之間的對比，如圖5所示。經(jīng)計算得出像點坐標差值組成的檢驗統(tǒng)計量已超過臨界值位于拒絕域內(nèi)，需要進一步分析，在分析階段，統(tǒng)計梯度方向角度45°～135°內(nèi)每5°一個區(qū)間的方向數(shù)量，如圖6所示。

圖3 試驗序列影像Fig.3 Experimental sequence images

圖4 未移動試驗序列影像Fig.4 Unmoved experimental sequence images

圖5 分散局部移動試驗序列影像Fig.5 Scattered local moving experimental sequence images

圖6 同名點梯度方向數(shù)目Fig.6 Number of gradient directions with same name

其中選定數(shù)量最多的前2組，進行第2步分析，依據(jù)程序計算2組中特征點的平均相鄰距離和影像間同名點的移動量，計算得出各自數(shù)值見表1。

表1 梯度方向數(shù)量前2組的距離Table 1 Top two distance of gradient directions

2組方向中發(fā)生滑移的特征點相鄰距離加權(quán)平均值遠大于所設(shè)定的3 pix，而影像間發(fā)生的位移也都遠大于5 pix。

2.2.3 移動影像組處理分析

第3組影像為初始影像和人工模擬移動影像對比，匹配點如圖7所示。

經(jīng)程序計算，像點坐標差值組成的檢驗統(tǒng)計量已遠大于t分布的拒絕閾臨界點，需要進一步分析，在分析階段，統(tǒng)計得出45°～135°內(nèi)按照每5°一個區(qū)間的方向數(shù)量(圖8)。

圖8 同名點方向梯度方向數(shù)目Fig.8 Number of gradient directions with same name

選定數(shù)量最多的前2組，進行第2步分析，依據(jù)程序計算2組中特征點的平均相鄰距離和影像間同名點的移動量，計算得出各自數(shù)值見表2。

表2 梯度方向數(shù)量前2組的距離Table 2 Top two distance of gradient directions

2組方向中發(fā)生滑移的特征點相鄰距離加權(quán)平均值小于方法中規(guī)定的3 pix，而影像間發(fā)生的滑移距離也都遠大于5 pix，同時根據(jù)相似幾何關(guān)系，按照拍攝距離3 m計算得出該滑移量大約為5.1 m。

2.3 結(jié)果分析

1)第1組試驗計算結(jié)果表明：影像間同名點像素坐標差平均值，其絕對值小于臨界值，說明不存在顯著性差異，造成的位移很大概率上是由于受到高斯噪聲的影響，因此可以判定未發(fā)生移動，并終止分析。

2)第2組試驗中，僅在較為分散的局部區(qū)域內(nèi)發(fā)生變動，其像素坐標差絕對值超出臨界值，經(jīng)計算得出兩影像間像點位移加權(quán)平均值超過設(shè)定閾值5 pix，而發(fā)生變化的點位相鄰間距加權(quán)平均值卻大于設(shè)定的3 pix，表明盡管點位發(fā)生了較為明顯變化，但點位之間的距離較遠，位置分散，區(qū)別于真實滑坡情況下的連片移動，故判定疑似人為所致，或滑坡前兆，發(fā)送提醒信號、不報警。

3)第3組試驗為已發(fā)生移動情況下的分析，該像素坐標差平均值遠大于臨界值，在分析階段，篩選出梯度方向數(shù)量最多的2組，計算得出變化點相鄰間距加權(quán)平均值小于設(shè)定值3 pix，同時影像間像點滑移量的加權(quán)平均值遠大于設(shè)定閾值5 pix，證明發(fā)生變化的點分布緊密且連接成片，目標位移也遠超過高斯噪聲和粗差的影響，故可判定目標發(fā)生滑坡移動，并報警。

3 結(jié) 論

1)提出了依靠單目序列影像進行滑坡位移的判斷算法，利用像素變化的統(tǒng)計規(guī)律，加入梯度方向、目標點相鄰距離和影像間目標點位移等參數(shù)，并設(shè)定初始篩選、數(shù)據(jù)分析、相應決策3個環(huán)節(jié)，提高了該算法的魯棒性。

2)該算法可對高斯噪聲影響、分散式局部移動、和大面積滑坡3種現(xiàn)象進行有效甄別，省略了計算三維空間坐標的繁瑣過程，提高了運行速度，節(jié)省生產(chǎn)成本和保證人身安全的前提下，對突發(fā)性滑坡事故應急工作能起到一定的幫助。

3)由于受試驗條件及時間限制，目前僅利用模擬試驗進行測試，且對于閾值設(shè)定的合理性和普適性還需要更多的實地場區(qū)進行驗證，這也是后期的優(yōu)化工作。