基于Matlab的概率積分法開采沉陷預計參數(shù)解算

沈 震 徐良驥 劉 哲 秦長才

(安徽理工大學測繪學院，安徽 淮南 232000)

基于Matlab的概率積分法開采沉陷預計參數(shù)解算

沈 震 徐良驥 劉 哲 秦長才

(安徽理工大學測繪學院，安徽 淮南 232000)

通過在開采工作面地表建立觀測站獲取地表變形數(shù)據(jù)解算概率積分法開采沉陷預計參數(shù)，從而預計工作面周邊或相似開采條件下工作面的開采沉陷并評估和指導開采作業(yè)，其前提是精確獲取概率積分法開采沉陷預計參數(shù)。為此，基于Matlab軟件，采用最小二乘法擬合觀測點變形數(shù)據(jù)解算概率積分法開采沉陷預計參數(shù)，并結合Matlab軟件繪圖工具開發(fā)了集數(shù)據(jù)載入、坐標轉(zhuǎn)換、參數(shù)解算、結果輸出、開采沉陷預計及反演為一體的可視化開采沉陷預計系統(tǒng)。以淮南謝橋煤礦11316工作面為例，根據(jù)地表移動觀測點位移數(shù)據(jù)解算開采沉陷預計參數(shù)并與實測值進行對比分析，結果表明，該系統(tǒng)解算出的開采沉陷預計參數(shù)符合兩淮礦區(qū)開采沉陷的基本規(guī)律，反演結果與實測值基本吻合，該系統(tǒng)對于實現(xiàn)礦區(qū)開采沉陷高精度預計和反演具有一定的參考價值。

開采沉陷 概率積分法 開采沉陷預計系統(tǒng) 最小二乘法 Matlab

開采沉陷會改變礦區(qū)原始地表狀態(tài)，對礦區(qū)環(huán)境及周邊居民生產(chǎn)生活造成影響，因此有必要對開采沉陷進行精確預計以便有計劃地進行開采作業(yè)[1-3]。開采沉陷的預計方法有概率積分法、負指數(shù)函數(shù)法，典型曲線法等，其中概率積分法應用較廣泛[4-6]。本研究借助Matlab軟件豐富的內(nèi)置函數(shù)與強大的數(shù)據(jù)處理功能，實現(xiàn)根據(jù)任意點變形量解算概率積分法開采沉陷預計參數(shù)，并結合該軟件的繪圖工具開發(fā)出可視化開采沉陷預計系統(tǒng)。

1 系統(tǒng)主要功能

1.1 概率積分法模型

以自定義函數(shù)的形式建立概率積分法模型，并保存在單獨的M文件中，該模型程序代碼如下：

Function Wz= probability_integration_method (data2，xy)

l=data1 (1)-data2 (3)-data2 (4);

L=data1 (2)-data2 (5)-data2 (6))*sin (data2 (7) +data1 (3))/sin (data2 (7));

CX=1/2*(erf(pi^(1/2)*xy(:,1)*data2(2)/data1(4))-erf(pi^(1/2)*(xy(:,1)-l)*data2(2)/data1(4)));

CY=1/2*(erf(pi^(1/2)*xy(:,2)*data2(2)/data1(5))-erf(pi^(1/2)*(xy(:,2)-L)*data2(2)/data1(6)))

Wmax=data1 (7)*data2 (1)*cos (data1 (3));

WZ=Wmax*CX.*CY.

(1)已知參數(shù)。包括L0、L1、a、H0、H1、H2和M。L0、L1分別為走向和傾向線長度，m；a為煤層傾角，(°)；H0為平均開采深度，m；H1、H2分別為下山及上山方向開采深度，m；M為煤層開采厚度，m。已知參數(shù)以向量的形式保存于“data1”中，調(diào)用時按順序記為data1(n)。

(2)預計參數(shù)。包括q、tanβ、S1、S2、S3、S4和θ0。q為下沉系數(shù)；tanβ為主要影響角正切值；S1、S2為走向線兩端拐點偏距，m；S3、S4為傾向線下山和上山方向拐點偏距，m；θ0為影響傳播角，(°)[7]。待求參數(shù)以向量形式保存于“data2”中。

(3)點坐標。地表觀測點在工作面坐標系中的坐標為(x,y)。

1.2 數(shù)據(jù)載入和坐標轉(zhuǎn)換

將所有觀測點的坐標值、下沉量和水平移動量按一定格式保存于Excel或.txt文檔中，通過“xlsread”、“textread”等命令讀取數(shù)據(jù)。Matlab軟件可將所有數(shù)據(jù)保存在一個矩陣中，按需要從該矩陣中提取坐標值、下沉量和平移量等信息，并保存在相應的數(shù)據(jù)矩陣中，便于調(diào)用[8]。

概率積分法通常采用獨立坐標系(工作面坐標系)預計開采沉陷。開采工作面多為矩形，以工作面的西南端點為原點，X軸沿工作面走向線方向，Y軸沿工作面傾向線方向并與X軸垂直，建立工作面坐標系，見圖1。圖1中，矩形ABCD為開采工作面，XOY為測量坐標系，X′O′Y′為工作面坐標系[9-10]。

圖1 工作面坐標系

坐標轉(zhuǎn)換時將點D在XOY坐標系中的坐標作為平移參數(shù)，在CAD軟件中量取工作面的傾斜角度作為旋轉(zhuǎn)參數(shù)(如無特殊要求可不考慮尺度參數(shù))。當工作面范圍較大時，也可根據(jù)2個或2個以上公共點在2套坐標系中的坐標，采用最小二乘法求解坐標轉(zhuǎn)換參數(shù)，坐標轉(zhuǎn)換模型為[11-12]

(1)

式中，(x′,y′)、(x,y)分別為點E在XOY測量坐標系和X′O′Y′工作面坐標系中的坐標；α為旋轉(zhuǎn)角度，(°)；x0、y0為平移參數(shù)，m。

將獲取的坐標轉(zhuǎn)換參數(shù)以向量形式保存于“p1”中，定義坐標轉(zhuǎn)換函數(shù)，程序代碼如下：

Function x1=coordinate_transformation (p1,X0)

for k=1:length(X0);

X1(k，:) = (X0(k，:)-[p1 (1)，p1 (2)]) ;

X1(k，:) = [cos (p1 (3))，sin (p1 (3));-sin (p1 (3))，cos (p1 (3))]*X1(k，:)′.

end

調(diào)用函數(shù)“coordinate_transformation”可完成由測量坐標系至工作面坐標系的轉(zhuǎn)換。

1.3 解算概率積分法預計參數(shù)

監(jiān)測點在工作面坐標系中的坐標為(x,y)，實際下沉量為z，采用最小二乘法進行曲面擬合，經(jīng)過多次迭代求取1組參數(shù)，使得擬合曲面與實測結果偏差的平方和最小，數(shù)學模型為[13-14]

(2)

式中，V1為各監(jiān)測點實測下沉量與最小二乘擬合值的偏差，m；V2為各監(jiān)測點實測水平變形量與最小二乘擬合值的偏差，m;k為監(jiān)測點數(shù)目；Wzk、Uk分別為各監(jiān)測點實測下沉和水平變形量，m；W(x,y),U(x,y)分別為各監(jiān)測點實測下沉和水平變形量的最小二乘擬合值，m，可調(diào)用Matlab軟件中的“l(fā)sqcurvefit”函數(shù)[15]進行擬合得到。程序代碼如下:

X1=coordinate_transformation (p1，X0);

ff=optimset; ff.MaxFunEvals=2^12; ff.TolFun=1e-10; ff.TolX=1e-10;

[data2,resnorm，residual]=lsqcurvefit (@probability_integration_method，data0，X1，Wz).

程序中，“X1”、“Wz”為監(jiān)測點在工作面坐標系下的坐標和對應的下沉量；“data0”為參數(shù)迭代初始值包含了7個變量，其格式與“data2”一致，初值的選取直接影響預計參數(shù)的解算精度，可參考開采工作面地質(zhì)及技術資料確定[16]。函數(shù)返回值包括預計參數(shù)的擬合最佳值(data2)、監(jiān)測點擬合殘差(residual)、監(jiān)測點擬合殘差的平方和(resnorm)，將解算所得參數(shù)作為已知值，根據(jù)水平移動量可解算監(jiān)測點水平移動參數(shù)。

1.4 開采沉陷反演及預計

解算出開采沉陷預計參數(shù)后，將已知參數(shù)和預計參數(shù)代入概率積分法模型計算開采沉陷區(qū)域內(nèi)的地表變形量，并采用Matlab軟件中的“ezmesh”、“ezcontour”函數(shù)繪制下沉曲面及下沉等值線圖，從而預計或反演地表變形結果，程序代碼如下[17-21]:

syms x; syms y;

xy=[x,y]

ezmesh (coordinate_transformation(data2,x,y),[min_x,max_x],[ min_y,max_y]);

ezcontour (coordinate_transformation(data2,x,y)，[min_x,max_x],[ min_y,max_y]);

程序中，“min_x”、“max_y”分別為X、Y軸上的繪圖邊界，繪圖結果見圖2、圖3。

圖2 預計下沉曲面

2 應用實例

2.1 開采沉陷預計參數(shù)解算

以淮南礦區(qū)謝橋礦11316工作面為例，工作面走向長1 980 m，傾向長240 m，煤層傾角15°，平均采深718.3 m，煤層開采厚度3.9 m。地表移動觀測站共設置2條觀測線，計87個測點，選取其中30個有較高精度并能反映工作面地表變形特點的監(jiān)測點實測數(shù)據(jù)作為解算數(shù)據(jù)，見表1。根據(jù)礦區(qū)預計參數(shù)經(jīng)驗值，取迭代初始值data0=[0.7;2;25;25;10;5; 1.31]進行最小二乘擬合，可得q=0.798，tanβ= 1.755，S1=35.97，S2=35.47，S3=9.996、S4=4.996，θ0=1.133。

圖3 地表下沉等值線

表1 11316工作面部分觀測點實測數(shù)據(jù)

Table1 Part of measured date at observation points of 11316 working face m

點 號XY下沉量ml01-14.41130.850.793ml0217.00139.570.885ml0347.85135.891.010ml0476.58136.241.131ml08225.37132.861.642ml10346.98135.541.670ml11377.39128.581.661ml17555.05126.511.457ml18585.66128.761.695ml21675.94120.781.717ml23826.03137.071.637ml24856.07134.051.661ml291215.86136.741.689ml301245.91139.921.702ml311275.38133.611.716ml331335.98137.941.638ml351396.63123.481.923ml371456.57138.791.668ml401577.31133.361.613ml411606.03126.611.648ml461876.96122.240.984ml471906.16131.590.817ml491965.44134.130.619ms03939.2621.011.519ms04944.1843.191.445ms07947.51102.591.614ms10950.25162.651.550ms11938.92183.511.583ms13954.14221.461.501ms14941.84242.971.411

2.2 預計參數(shù)驗證及地表移動變形規(guī)律反演

將開采沉陷預計參數(shù)代入概率積分法模型計算上述30個監(jiān)測點的地表變形量的最小二乘擬合值，并與實測數(shù)據(jù)進行對比，計算擬合殘差，結果見表2。

表2 觀測點擬合殘差

Table.2 Fitting residual of observation points m

點號殘差點號殘差點號殘差ml01-0.01ml230.04ml46-0.01ml020.02ml240.02ml470.03ml030.02ml29-0.01ml49-0.01ml040.01ml30-0.02ms03-0.13ml08-0.10ml31-0.03ms040.05ml10-0.02ml330.04ms070.05ml110.01ml35-0.24ms100.10ml170.23ml370.01ms110.02ml18-0.01ml400.03ms13-0.04ml21-0.03ml41-0.02ms14-0.05

根據(jù)擬合值和擬合殘差計算擬合決定系數(shù)，公式為

(3)

式中，SSR為擬合值平方和，m2；SSE為殘差平方和，m2。

將相關數(shù)值代入式(3)，計算得R2=0.997 5，擬合程度較好，擬合值與實測值(表1)較為接近。

采用概率積分法對11316工作面開采沉陷進行反演，并將反演結果與實測值進行對比，結果見表3。

表3 概率積分法反演結果與實測值比較

由表3可知，由反演得到的最大下沉值與實測值的誤差僅為0.242 m，其余變形參數(shù)反演值和實測值也較為接近。對原始觀測數(shù)據(jù)進行分析可知，達到充分下沉的下沉盆地內(nèi)監(jiān)測點實測下沉量為1.60～1.72 m，最大下沉值(1.923 m)對應的是ml35點，為異常點。據(jù)此可認為，本研究開發(fā)的開采沉陷預計系統(tǒng)能夠滿足礦區(qū)地表沉陷預計的精度要求。

3 結 語

基于Matlab軟件，開發(fā)了開采沉陷預計系統(tǒng)，該系統(tǒng)采用最小二乘法擬合觀測點變形數(shù)據(jù)解算概率積分法開采沉陷預計參數(shù)。淮南謝橋煤礦11316工作面試驗結果表明，該系統(tǒng)對于礦區(qū)地表沉陷的預計精度較高。

[1] 何國清，楊 倫，凌賡娣，等.礦山開采沉陷學[M].徐州：中國礦業(yè)大學出版社，1995. He Guoqing，Yang Lun，Ling Gengdi，et al.Mining Subsidence in Mine [M].Xuzhou:China University of Mining and Technology Press，1995.

[2] 方 齊，郭廣禮，朱曉峻.礦區(qū)地表觀測站輔助設計與數(shù)據(jù)處理程序?qū)崿F(xiàn)[J].金屬礦山，2015(5)：129-134. Fang Qi，Guo Guangli，Zhu Xiaojun.Realization of aided design of surface movement observation station and data processing in mining area[J].Metal Mine，2015(5)：129-134.

[3] 楊俊凱，范洪冬，趙偉穎，等.基于D-InSAR技術和灰色Verhulst模型的礦區(qū)沉降監(jiān)測與預計[J].金屬礦山，2015(3)：143-147. Yang Junkai,Fan Hongdong,Zhao Weiying,et al.Monitoring and prediction of mining subsidence based on D-InSAR and gray verhulst model[J].Metal Mine，2015(3)：143-147.

[4] 姚 琦，馮 濤，李石林，等.基于概率積分法的煤礦“三下”開采沉陷預計[J].煤礦安全，2012(7)：188-190. Yao Qi，F(xiàn)eng Tao，Li Shilin，et al.Subsidence prediction of coal mine “three under” mining based on probability integral method [J].Safety in Coal Mines，2012(7):188-190.

[5] 朱曉峻，郭廣禮，方 齊.概率積分法預計參數(shù)反演方法研究進展[J].金屬礦山，2015(4)：173-177. Zhu Xiaojun，Guo Guangli，F(xiàn)ang Qi.Recent progress on parameter inversion of probability integral method [J].Metal Mine，2015(4):173-177.

[6] 楊俊凱，陳炳乾，鄧喀中，等.基于D-InSAR與概率積分法的開采沉陷監(jiān)測與預計[J].金屬礦山，2015(4)：195-200. Yang Junkai，Chen Bingqian，Deng Kazhong，et al.Monitoring and prediction of mining subsidence based on D-InSAR and probability integral method [J].Metal Mine，2015(4):195-200.

[7] 王永輝，倪岳暉，周建偉，等.基于概率積分法的橫河煤礦巨厚松散層下開采沉陷預測分析[J].地質(zhì)科技情報，2014(4)：219-224. Wang Yonghui，Ni Yuehui，Zhou Jianwei.et al.Subsidence prediction under thick and loose overburden of Henghe coal mine based on probability integral method [J].Geological Science and Technology Information，2014(4):219-224.

[8] 劉保柱.MATLAB 7.0從入門到精通[M].北京：人民郵電出版社，2006. Liu Baozhu.MATLAB 7.0 From Beginning to Master [M].Beijing:Post & Telecom Press，2006.

[9] 韓買俠，陳傳錄，王夏莉.基于Matlab的坐標轉(zhuǎn)換方法研究[J].測繪與空間地理信息，2013(9)：89-91. Han Maixia，Chen Chuanlu，Wang Xiali.Research of coordinate transformation methods based on Matlab [J].Geomatics & Spatial Information Technology，2013(9):89-91.

[10] 汪桂生.礦區(qū)開采沉陷觀測數(shù)據(jù)處理研究[D].西安：西安科技大學，2011. Wang Guisheng.Observation Data Processing of Mining Subsidence[D].Xi′an:Xi′an University of Science and Technology，2011.

[11] 梁月吉，謝劭峰，龐光鋒.基于Matlab的坐標轉(zhuǎn)換程序設計[J].地理空間信息，2014(2)：124-125. Liang Yueji，Xie Shaofeng，Pang Guangfeng.Design of coordinate transform program design based on Matlab[J].Geospatial Information，2014(2):124-125.

[12] 劉學軍.工程測量中平面坐標轉(zhuǎn)換及其應用[J].北京測繪，2014(6)：142-145. Liu Xuejun.Planar coordinate conversion and application in engineering surveying [J].Beijing Surveying and Mapping，2014(6):142-145.

[13] 許 冬，王臨清，吳 侃.任意形狀工作面沉陷預計計算方法[J].金屬礦山，2014(5)：55-59. Xu Dong，Wang Linqing，Wu Kan.Mining subsidence prediction calculation method of random shape working face [J].Metal Mine，2014(5):55-59.

[14] 王 友，王雙亭.抗差嶺估計在概率積分法預計參數(shù)求取中的應用研究[J].煤礦開采，2014(3)：17-19. Wang You，Wang Shuangting.Application of anti-poor-estimation in prediction parameters solution by probability integral method [J].Coal Mining Technology，2014(3):17-19.

[15] 譚衍濤，張興福.MATLAB擬合工具箱在GNSS高程轉(zhuǎn)換中的應用[J].工程勘察，2014(9)：65-68. Tan Yantao，Zhang Xingfu.Application of MATLAB surface fitting tool for GNSS height transformation [J].Geotechnical Investigation & Surveying，2014(9):65-68.

[16] 鄒友峰.開采沉陷預計參數(shù)的確定方法[J].焦作工學院學報：自然科學版，2001，20(4)：253-257. Zou Youfeng.Determining method of predication parameters on mining subsidence [J].Journal of Jiaozuo Institute of Technology:Natural Science Edition，2001，20(4):253-257.

[17] 劉 蕓.淺析Matlab繪圖功能在高等數(shù)學中的應用[J].吉林省教育學院學報，2012(7)：57-58. Liu Yun.Application of Matlab drawing function in higher mathematics[J].Journal of Educational Institute of Jilin Province，2012(7):57-58.

[18] 英 英.基于MATLAB的圖形圖像處理系統(tǒng)的實現(xiàn)[D].呼和浩特：內(nèi)蒙古大學，2013. Ying Ying.Realization of the Graphics Image Processing System Based on Matlab[D].Hohhot:Inner Monggulia University，2013.

[19] 徐光華.基于二次曲線擬合的隧道激光點云濾波方法及其應用[J].測繪通報，2015(5)：42-45. Xu Guanghua.Laser point cloud filtering and application in tunnel deformation monitoring based on quadratic curves fitting [J].Bulletin of Surveying and Mapping，2015(5):42-45.

[20] 顧 偉，譚志祥，鄧喀中.基于Active X Automation的三維可視開采沉陷預計[J].金屬礦山，2013(3)：119-122. Gu Wei，Tan Zhixiang，Deng Kazhong.3D visualization prediction of mining subsidence system based on Active X Automation [J].Metal Mine，2013(3):119-122.

[21] 查劍鋒，郭廣禮，趙海濤，等.概率積分法修正體系現(xiàn)狀及發(fā)展展望[J].金屬礦山，2008(1)：15-18. Zha Jiangeng，Guo Guangli，Zhao Haitao，et al.Present situation and prospect of correction system for probability integral method [J].Metal Mine，2008(1):15-18.

(責任編輯 王小兵)

Calculating on the Prediction Parameters of Mining Subsidence with Probability Integral Method Based on Matlab

Shen Zhen Xu Liangji Liu Zhe Qin Changcai

(SchoolofGeodesyandGeomatics，AnhuiUniversityofScience&Technology，Huainan232000，China)

The prediction parameters of mining subsidence are calculated by the probability integral method of the surface deformation data obtained by the observation stations established on the surface of mining working face，and the mining subsidence of the surrounding area of mining working face or the mining working face under the similar mining conditions are predicted to conduct evaluation and guidance on the mining operation.Therefore，obtaining the prediction parameters of mining subsidence with high precision is the precondition of prediction and inversion of mining subsidence prediction.Based on Matlab software，the deformation data of observation points are fitted by the least square method to calculate the prediction parameters of mining subsidence with probability integral method.Combing with the drawing tools of Matlab software，the mining subsidence visualization system with the functions of data loading，coordinate transformation，calculation parameters，the output of calculation results，prediction and inversion of mining subsidence is developed.Taking the 11316 working face of Xieqiao coal mine，Huainan city as the research example，according to the deformation data of the observation points of the surface observation stations，the prediction parameters of mining subsidence are calculated，and comparison between predicted parameters and measured parameters is conducted.The results show that the predicted parameters calculated by the mining subsidence visualization system is consistent with the basic rules of the mining subsidence in Huainan & Huaibei mining area，and the inversion results of the mining subsidence visualization system is identical to the measured parameters.The mining subsidence visualization system can provide reference for realizing the prediction and inversion of mining subsidence with high precision.

Mining subsidence，Probability integral method，Mining subsidence prediction system，Least square method，Matlab

2015-06-04

安徽省對外科技合作計劃項目(編號：1503062020)，國家自然科學基金項目(編號：41472323)。

沈 震(1990—)，男，碩士研究生。通訊作者 徐良驥(1978—)，男，副教授，博士，碩士研究生導師。

TD325

A

1001-1250(2015)-09-170-05