摘要：本文利用MATLAB 工具箱函數，設計了求取地表移動變形預計參數的程序，結合工作面任意點的觀測數據和地質采礦條件，能夠很好的求取地表移動變形預計參數。本文通過實例，對所擬合的參數進行檢驗，證明了此方法的可行性。

　　關鍵字：MATLAB,非線性最小二乘法,概率積分法,預計參數

　　1 引言

　　MATLAB有美國MathWorks公司推出，其與Mathematica、Maple并稱為三大數學軟件，具有優秀的數值計算能力和卓越的數據可視化功能。 MATLAB自問世以來，已經發展稱為適合多學科、多種工作平臺的功能強大的大型軟件，廣泛應用于大學科研、工程計算等領域，尤其在工程界。

　　概率積分法通過積分而得到下沉曲線表達式，因其所用的移動和變形預計公式中含有概率積分而得名。這種方法是將礦山巖層移動作為一種服從統計規律的隨機現象來討論。因此，此種方法是以隨即介質理論為基礎的一種預計方法。概率積分法自20世紀60年代引入我國以來，在許多礦山得到廣泛使用，并不斷發展，目前是我國進行地表移動變形預計主要的方法之一。地表移動變形預計參數對開采沉陷預計的準確性具有重要意義，因而本文以概率積分法數學模型為基礎，借助MATLAB平臺，編寫求取地表移動變形參數的程序。

　　2 概率積分法求參的數學模型

　　利用概率積分法關于任意點的地表下沉預計公式作為求取預計參數的基本數學模型，即：

　　地表移動觀測站任何一個測點的下沉實測值W全能表達成自變量(x,y)(觀測點的水平坐標)和概率積分法預計參數q、tgβ、s3、s4、tgβ1、tgβ2、s1、s2、θ的函數，見下式：

　　假設公式(3)里的參數q、tgβ、s3、s4、tgβ1、tgβ2、s1、s2、ctgθ，已經求出，那么上面的數學公式中就只剩下水平移動系數b這一待求參數。

　　3 Matlab曲線擬合求參函數

　　Matlab強大的繪圖和用戶自定義函數等功能，在數據處理領域被廣泛應用，其中當然包括對實驗、實測數據的擬合工作。鑒于本文研究內容，著重介紹Matlab中非線性曲線擬合(非線性最小二乘擬合)函數。

　　Matlab中提供了兩個求非線性最小二乘擬合的函數：lsqcurvefit和lsqnonlin。lsqcurvefit的定義如下：

　　已知擬合的數據點為：xdata和ydata，lsqcurvefit用以求含參量variable(向量)的向量值函數F(x,xdata)中的參變量x(向量)，使滿足下式：

　　完整調用格式為：

　　式中：x是擬合的參數最優解，norm是最優解的誤差的平方和，res為誤差向量，ef是程序結束時的狀態指示：ef>0表示結果收斂，ef=0表示函數調用次數或迭代次數達到最大值(該值在options中指定)，ef<0表示不收斂，out則包含數據的一個結構變量，包含實際的迭代次數等，lam是上下界所對應的Lagrange乘子，jac是結果(x點)處的雅可比矩陣，x0為初始參數(缺省時程序自動取x0=0)fun是給出的目標函數，xdata和ydata為擬合的數據，v1,v2是擬合參數的上下界，可以用兩個[]代替，options是包含算法控制參數的結構設定(或顯示)控制參數的命令為Optimset函數。

　　Matlab可以很快速的直接讀寫txt格式和Excel格式的數據文件，所以給開采沉陷數據處理帶來了方便，解決了求參和預計過程中繁多的數據手工輸入和輸出的問題。在Matlab中讀寫Excel文件可以使用xlsread函數和xlswrite函數，xlsread函數的使用格式為：[num txt]=xlsread(‘filename’，sheet，’range’)，讀一個Excel文件名為filename的文件，sheet用來指定頁，range用來指定區域，最后把指定的數值讀到num中，所有字符串到字符串單元數組txt中。使用xlswrite函數寫Excel文件，經常用到的格式為：xlswrite (‘filename’，M，sheet，’range’)，寫矩陣或字符串單元數組M到filename中的指定頁和指定區域。

　　4 地表移動變形參數求取與檢驗

　　4.1地表移動變形預計參數求取的算法流程

　　4.2求取沉陷預計參數的程序設計

　　本文主要用到的兩個自定義函數budgets和horizonal_budgets，函數budgets用于求取沉陷預計參數q、tgβ、s3、s4、tgβ1、tgβ2、s1、s2、θ，部分源代碼如下：

　　options=optimset('Display','iter','MaxFunEvals',2^8,'MaxIter',2^14,'TolX',1e-3,'TolFun',1e-4,'LargeScale','off','LevenbergMarquardt','on','LineSearchType','cubicpoly');

　　f=inline('1./(a(1)*x(3,4)*cos(x(3,5))).*((a(1)*x(3,4)*cos(x(3,5)))/2.*(erf(sqrt(pi).*x(1,:).*a(2)./x(3,1))-erf(sqrt(pi).*(x(1,:)-(x(3,6)-a(5)-a(6))).*a(2)./x(3,1)))).*((a(1)*x(3,4)*cos(x(3,5)))/2.*(erf(sqrt(pi).*x(2,:).*a(3)./x(3,2))-erf(sqrt(pi).*(x(2,:)-(x(3,7)-a(7)-a(8)).*(sin(a(9)+x(3,5))/sin(a(9)))).*a(4)./x(3,3))))','a','x');

　　[a,norm,res,ef,out,lam,jac]=lsqcurvefit(f,a,[x;y;z],w,[],[],options);

　　函數horizonal_budgets用于求取沉陷預計參數b，部分源代碼如下：

　　options=optimset('Display','iter','MaxFunEvals',2^8,'MaxIter',2^14,'TolX',1e-8,'TolFun',1e-12,'LargeScale','on','LevenbergMarquardt','on')

　　f=inline('b.*(x(4,1).*x(3,4).*cos(x(3,5))).*(exp(-pi.*(x(1,:).*x(4,2)./x(3,1)).^2)-exp(-pi.*((x(1,:)-(x(3,6)-x(4,5)-x(4,6))).*x(4,2)./x(3,1)).^2)).*(x(4,1)*x(3,4)*cos(x(3,5))/2.*(erf(sqrt(pi).*x(2,:).*x(4,3)./x(3,2))-erf(sqrt(pi).*(x(2,:)-(x(3,7)-x(4,7)-x(4,8)).*(sin(x(4,9)+x(3,5))/sin(x(4,9)))).*x(4,4)./x(3,3))))./(x(4,1)*x(3,4)*cos(x(3,5))).*cos(x(3,8))+b.*(x(4,1)*x(3,4)*cos(x(3,5))).*(exp(-pi.*(x(2,:).*x(4,3)./x(3,2)).^2)-exp(-pi.*((x(2,:)-(x(3,7)-x(4,7)-x(4,8))).*x(4,3)./x(3,3)).^2)).*(x(4,1)*x(3,4)*cos(x(3,5))/2.*(erf(sqrt(pi).*x(1,:).*x(4,2)./x(3,1))-erf(sqrt(pi).*(x(1,:)-(x(3,6)-x(4,5)-x(4,6))).*x(4,2)./x(3,1))))./(x(4,1)*x(3,4)*cos(x(3,5))).*sin(x(3,8))','b','x');

　　[b,norm,res,ef,out,lam,jac]=lsqcurvefit(f,b,[x;y;z;w],u,[],[],options);

　　特別注意的是在這兩個函數中都調用了LevenbergMarquardt(L-M)極小化方法，以避免出現局部解和病態問題。

　　4.3 預計參數的實例驗證

　　本文以皖北某礦井1013工作面為例，結合地質采礦條件和觀測點實測數據，利用MATLAB提供的函數求取地表移動與變形預計參數。此工作面走向長度575m，傾向長度150m，下山采深為415m，上山采深為361m，平均采深為388m，煤層平均傾角10o，煤層平均厚度3.1m。在該工作面上方地表布置兩條觀測線(走向和斜向)，由于受實地地形的限制，走向設置大半條觀測線，斜向觀測線觀測點布設在水渠上，走向觀測線與斜向觀測線的夾角為44o。地表移動觀測站布設如圖2所示,工作面實測數據整理如圖3所示。

　　預計的走向觀測線下沉值和水平移動值的中誤差分別為60.01mm和21.4mm，分別占最大值的2.6%和2.5%，預計的斜向觀測線下沉值和水平移動值的中誤差分別為73.3mm和34.4mm，分別占最大值的3.1%和4.5%。

　　5 結論

　　通過預計的擬合曲線與實測數據的可視化比較，可知只要充分利用MATLAB工具箱函數，結合任意工作面監測點實測數據及實際的地質采礦條件，就能夠擬合出較為準確的地表移動變形預計參數，極為方便，但在此需特別指出，擬合時初始化的參數非常重要，因此要充分利用周圍礦區的資料。

　　參考文獻：

　　[1] 鄒友峰,鄧喀中,馬偉民.礦山開采沉陷工程[M].徐州:中國礦業大學出版社,2003

　　[2] 吳侃,周鳴.礦山開采沉陷預測預報系統[M].徐州:中國礦業大學出版社,2000

　　[3] 康建榮,王金莊,溫澤民.任意形多工作面多線段開采沉陷預計系統(MSPS)[J]礦山測量,2000,1(3):24~27

　　[4] 朱劉娟,陳俊杰,鄒友峰.任意形狀工作面開采地表移動變形預計的算法實現[J].遼寧工程技術大學學報,2005,24(3):337~339

　　[5] 陳良松,汪青松,楊永國. 巨厚松散層條件下地表移動變形特征.煤礦安全,2011 42(8)

　　[6] 王明柱,郭廣禮,王磊,王彬.基于嶺估計的概率積分法預計參數的求取.煤礦開采, 2012 17(2)

　　[7] 周萬茂,張華興,仲惟林.多工作面開采影響下求取地表移動參數方法研究.煤炭學報.2000 25(3)

　　[8] 梁明,王成緒.厚黃土覆蓋山區開采沉陷預計.煤田地質與勘探.2001 29(2)

　　[9] 朱廣軼,王瑋,劉曉群.開采與地表沉陷分析的可視化軟件系統.煤炭工程.2006 (8)