李昱昊，安士凱，周大偉，詹少奇，高銀貴

（1.中國礦業(yè)大學(xué)江蘇省資源環(huán)境信息工程重點實驗室，江蘇 徐州 221116；2.平安煤炭開采工程技術(shù)研究院有限責任公司煤礦生態(tài)環(huán)境保護國家工程實驗室，安徽 淮南 232001；3.鄂爾多斯市華興能源有限責任公司，內(nèi)蒙古鄂爾多斯 010399）

我國煤炭資源賦存現(xiàn)狀是“西部多、中部富、東部區(qū)域枯竭”。以晉陜蒙為代表的西部地區(qū)，煤炭產(chǎn)量已占到全國總產(chǎn)量的70%以上。隨著大量煤炭的產(chǎn)出，所引起的地表沉陷及災(zāi)害問題也日益突出，如地表塌陷和塌陷區(qū)積水[1-2]。因此，對開采沉陷進行快速、準確、全面的監(jiān)測，合理地建立下沉盆地模型和反演參數(shù)，可以定量研究受開采影響的巖層的時空變化規(guī)律?，F(xiàn)有的GNSS、水準等大地測量方法野外工作量大，且不足以獲取整個沉陷區(qū)的變形特征，因此難以支撐大規(guī)模西部礦區(qū)煤炭開采沉陷監(jiān)測。InSAR[3-8]等遙感檢測方法在變形劇烈的區(qū)域容易出現(xiàn)是失相干的現(xiàn)象，許多高分辨率影像成本過高。無人機攝影測量技術(shù)[9-12]以其自動化、智能化、快速化的特點，通過對地面周期性快速航拍獲取高分辨影像，構(gòu)建密集點云，從而在礦區(qū)沉陷監(jiān)測中發(fā)揮巨大優(yōu)勢；但在產(chǎn)品生產(chǎn)過程中，點云常常出現(xiàn)噪聲，野點的情況，同時因其精度僅達到分米級，對礦區(qū)盆地參數(shù)反演產(chǎn)生極大的影響[13]。在沉陷參數(shù)預(yù)計方面，目前有很多求取參數(shù)的研究方法[14]，但存在或多或少的問題；正交參數(shù)[15]試驗當參數(shù)較多時，工作量較大，試驗次數(shù)較多；模矢法[16]容易把局部最優(yōu)解當成全局最優(yōu)解，對參數(shù)初值有一定的要求；遺傳算法[17-18]需要確定種群數(shù)量等參數(shù)，否則會影響收斂速度與參數(shù)精度，且存在局部搜索能力差的缺陷。模擬退火算法[19-21]（Simulated Annealing，SA）是由Kirkpatrick Z 提出的一種求解組合優(yōu)化問題的算法，其優(yōu)點是：根據(jù)Metropolis 接受準則，采用一定的概率接受比當前解差的解，從而跳出局部最優(yōu)解來搜索全局最優(yōu)解。

針對西部礦區(qū)地形起伏且植被較少的地理環(huán)境，本研究利用索尼a7r2 數(shù)碼相機搭載垂直起降固定翼無人機開展礦區(qū)沉陷監(jiān)測建模研究。以內(nèi)蒙古唐家會煤礦工作面地表為實驗區(qū)，利用無人機獲取2 期影像數(shù)據(jù)，制作DEM 并獲取下沉盆地，針對存在噪點問題，對全盆地點云去噪方法進行對比擇優(yōu)。同時，為了使求參結(jié)果更準確，融合概率積分與模擬退火算法，選取全盆地下沉值進行盆地參數(shù)反演，最后對沉陷參數(shù)進行抗差分析[22-23]。

1 方法原理

1.1 無人機攝影測量獲取下沉盆地原理

無人機攝影測量技術(shù)流程圖如圖1。

在無人機攝影測量技術(shù)應(yīng)用過程中，為了確保礦區(qū)監(jiān)測工作的順利進行，首要任務(wù)是相機標定，設(shè)置相控點布設(shè)方案以及無人機飛行路線的規(guī)劃處理。在具體過程中，需要明確監(jiān)測區(qū)域邊界，設(shè)定航向重疊率、旁向重疊率、飛行航高，以此為基礎(chǔ)開展礦區(qū)監(jiān)測工作。

在獲取區(qū)域影像數(shù)據(jù)后，結(jié)合無人機飛控系統(tǒng)導(dǎo)出的POS 數(shù)據(jù)和校正過后得到的相機參數(shù)數(shù)據(jù)，進行影像自動內(nèi)定向，然后依次進行輸出坐標系、快速檢測、導(dǎo)入像控點、刺相控點、高精度處理，生成高密度點云，最終生成測區(qū)數(shù)字高程模型。通過多期對該區(qū)域的重復(fù)監(jiān)測，獲取不同時期工作面推進到不同位置的DEM，將DEM 相減得到監(jiān)測區(qū)域的地表下沉盆地[24]。

1.2 基于模擬退火算法的概率積分求參原理

根據(jù)概率積分法原理[25-26]，若煤層是水平的，單元B（s，t）的開采所引起的地表任意一點A 的下沉We（x，y）為：為煤層走向長為，L 為煤層傾向水平長。

模擬退火算法反演概率積分參數(shù)步驟如下：

1）給出概率積分的初始參數(shù)，確定初始溫度T，最低溫度T′，降溫系數(shù)δ 等退火參數(shù)，計算初始解。

2）判斷T＜T′，若是，則輸出結(jié)果X，若不是，執(zhí)行步驟3）～步驟5）。

3）若達到迭代長度，則降溫T=Tδ，δ 為介于0～1的實數(shù)。

4）生成相鄰方案。假設(shè)在某一狀態(tài)下參數(shù)值為X，則它的相鄰狀態(tài)X′為X 在一定鄰域隨機變化產(chǎn)生的值，并且該值必須滿足上下限約束。

式中：W0為煤層的頂板下沉量，W0=mqcosα；m為實際開采厚度；q 為下沉系數(shù)；α 為煤層傾角；D

2 內(nèi)蒙古唐家會煤礦概率積分法求參

2.1 獲取唐家會下沉盆地

唐家會區(qū)域影像圖圖2。

圖2 唐家會區(qū)域影像圖Fig.2 Image map of Tangjiahui area

唐家會煤礦位于內(nèi)蒙古自治區(qū)準格爾煤田東孔兌普查區(qū)的西南部，其地理坐標為：東經(jīng)：111°10′27″～111°14′34″，北緯：39°52′45″～39°57′22″。選擇61202工作面作為觀測地點，該工作面尺寸為240 m×1 960 m，平均采深538 m，平均采厚19 m，平均傾角2°，計劃開采時間2020 年7 月，計劃結(jié)束時間2021 年7月，在測區(qū)東南部有矸石山堆積。

用垂直起降固定翼無人機于2020 年8 月采集第1 期無人機影像數(shù)據(jù)，2021 年3 月采集第2 期影像數(shù)據(jù)，監(jiān)測區(qū)域?qū)捈s為1 000 m，長約為2 900 m，面積約為3 km2。其中航向重疊率、旁向重疊率均設(shè)為80%，航高設(shè)置為140 m。在工作面四角、航飛邊緣四角布設(shè)通視較好的像控點，其他像控點均勻布設(shè)在旁向重疊范圍內(nèi)，共計29 個。第1 次數(shù)據(jù)采集共獲得2 079 張有效影像，第2 次數(shù)據(jù)采集共獲得2 531 張有效影像，通過內(nèi)業(yè)處理獲取DEM，無人機影像處理生成DEM 如圖3。用檢查點驗證DEM 精度，經(jīng)計算，第1 期DEM 中誤差為0.066 m，第2 期DEM 中誤差為0.059 m，用2 期DEM 相減獲得唐家會地區(qū)地表下沉盆地，地面下沉盆地如圖4。

圖3 無人機影像處理生成DEMFig.3 The DEM generated by UAV image processing

圖4 地面下沉盆地Fig.4 Surface subsidence basin

2.2 點云去噪方法對比

經(jīng)過濾波后的DEM 仍帶有噪點，將DEM 復(fù)原為下沉點云，下沉點云圖如圖5，其中噪點區(qū)域用圓圈圈出。

圖5 下沉點云圖Fig.5 Subsidence point cloud

產(chǎn)生噪點的主要原因有：

1）點云濾波算法的局限性。植被點以及低矮建筑物不能完全去除。在研究區(qū)域里，由于2 期DEM監(jiān)測時間在夏冬兩季，植被區(qū)域覆蓋不同，疊加后的DEM 偏差較大。

2）插值引起的誤差。無人機穿透能力較低。在研究區(qū)域西部方位有落差較大的溝壑，其覆蓋植被較多，無人機攝影測量難以穿透此區(qū)域。不含點云數(shù)據(jù)的區(qū)域依賴于DEM 插值。在點云密度較低的區(qū)域，其差值DEM 存在較多偏差。

3）無人機精度原因。無人機精度為分米級，在保證整體精度較高的情況下，有少數(shù)點高程誤差較高，且水平偏移引起的誤差較為明顯。

這些噪點對后續(xù)處理或者預(yù)計參數(shù)效果帶來了很大的不穩(wěn)定性，因此需要進行去噪處理。BP 神經(jīng)網(wǎng)絡(luò)是一種按誤差反向傳播（簡稱誤差反傳）訓(xùn)練的多層前饋網(wǎng)絡(luò)，其算法稱為BP 算法，它的基本思想是梯度下降法，利用梯度搜索技術(shù)，以期使網(wǎng)絡(luò)的實際輸出值和期望輸出值的誤差均方差為最小。針對沉陷點云去噪，選用MATLAB 中BP 神經(jīng)網(wǎng)絡(luò)工具箱對原始下沉點云進行去噪處理，設(shè)置訓(xùn)練次數(shù)為1 000，學(xué)習(xí)速率為0.05，設(shè)置顯示頻率為500 次顯示1 次，訓(xùn)練目標最小誤差為10-5。同時選擇中值去噪、最小二乘去噪、最近鄰去噪與BP 神經(jīng)網(wǎng)絡(luò)去噪對比，其中中值濾波窗口設(shè)置為20，最小二乘濾波窗口設(shè)置為3，最近鄰點數(shù)設(shè)置為10，離群點閾值設(shè)置為1。去噪對比圖如圖6。

圖6 去噪對比圖Fig.6 Comparison of denoise

經(jīng)過4 種去噪方法的測試對比發(fā)現(xiàn)：最小二乘去噪效果最差，去噪窗口增大后容易出現(xiàn)點云紊亂的現(xiàn)象，其次為中值去噪及最近鄰去噪算法，BP 神經(jīng)網(wǎng)絡(luò)去噪效果最好，同時保留了下沉盆地數(shù)據(jù)的完整性。

2.3 全盆地擬合求參

61202 工作面走向達到充分采動，傾向方向為非充分開采，因此計算時走向需要乘上傾向采動系數(shù)。為使結(jié)果更可靠，均勻選取全盆地觀測值作擬合初值，避免選點帶來的參數(shù)不穩(wěn)定性，其中由于東南出矸石山在工作面開采范圍內(nèi)，因此選點時避開矸石山區(qū)域，共計取到368 個點。根據(jù)概率積分法，在Matlab 中多次使用模擬退火算法預(yù)計，選取目標函數(shù)值最小的參數(shù)作為最終參數(shù)，最終求得：下沉系數(shù)q 為0.6，走向主要影響角正切tanβ 為0.6，傾向下山方向主要影響角正切tanβ1為1.337，傾向上山方向主要影響角正切tanβ2為3.43，開采影響傳播角為89°。根據(jù)預(yù)計參數(shù)畫出下沉等值線，并模擬下沉盆地，模擬下沉盆地圖如圖7。

圖7 模擬下沉盆地圖Fig.7 Simulated subsidence basin graph

用配對T 檢驗判斷2 組數(shù)據(jù)相似程度，在執(zhí)行T 檢驗之前需要驗證實測值與預(yù)計值的誤差是否符合或近似符合正態(tài)分布，因此隨機選取40 組數(shù)據(jù)執(zhí)行Shapiro-Wilk（W 檢驗），利用本次實驗樣本計算得出檢驗統(tǒng)計量Wα=0.928。取顯著性水平0.05，根據(jù)樣本數(shù)量查表得到界限值為0.980 131，可知Wα=0.928＜0.981 031，從而在該顯著性水平下，接受原假設(shè)，即樣本服從正態(tài)分布，實測數(shù)據(jù)與預(yù)計數(shù)據(jù)間的誤差為偶然誤差，對實測值與預(yù)計值做配對T 檢驗，得到結(jié)果為0.704＞0.05，配對數(shù)據(jù)無顯著性差異，擬合結(jié)果較好。

3 UAV 求取參數(shù)的可靠性分析

實測數(shù)據(jù)由于存在誤差，與真實模擬盆地存在一定的差距，而誤差也會對求參結(jié)果帶來很大的影響，因此需要對參數(shù)進行抗差分析。各參數(shù)對誤差的敏感度不同，有些參數(shù)敏感度很高，而有些參數(shù)在誤差很大時也不會產(chǎn)生很大的影響。選取唐家會模擬盆地，分別加入不同中誤差來對比參數(shù)的影響程度，對每一次加入的誤差多次求參，優(yōu)選出擬合程度最好的參數(shù)來作為在該誤差時的參數(shù)，全盆地加入隨機誤差后的求參結(jié)果見表1。

由表1 可得，當對實測下沉值加入約（1%～10%）W0，即100～700 mm 的測量中誤差時，求取的下沉系數(shù)q、主要影響角正切tanβ、均與參數(shù)真值吻合度較高，比較穩(wěn)定。當對實測下沉值加入超過10%W0，即800 mm 中誤差時，下沉系數(shù)與主要影響角正切誤差較大，故認為此時的求參結(jié)果不可靠，而原始擬合下沉中誤差為589 mm，占最大下沉值8.1%，因此求參結(jié)果可靠。

表1 全盆地加入隨機誤差后的求參結(jié)果Table 1 Parameters acquired after adding random error to the whole basin

4 結(jié) 語

以內(nèi)蒙古唐家會煤礦為例，用BP 神經(jīng)網(wǎng)絡(luò)算法對下沉點云做去噪處理，并與中值去噪、最小二乘去噪、最近鄰去噪方法做對比，提高了下沉模型精度。針對傳統(tǒng)依靠走向、傾向主斷面求參精度不高的情況，融合模擬退火算法與概率積分模型，選用全盆地特征點進行參數(shù)反演，并對參數(shù)結(jié)果作抗差分析。結(jié)果表明，擬合模型測量中誤差占最大下沉值的8.1%，預(yù)計結(jié)果滿足工程要求，克服了無人機精度不高對求參帶來的影響，為后續(xù)采煤工作、災(zāi)害治理提供技術(shù)支持。