摘要：本文利用MATLAB 工具箱函數(shù)，設(shè)計(jì)了求取地表移動(dòng)變形預(yù)計(jì)參數(shù)的程序，結(jié)合工作面任意點(diǎn)的觀測(cè)數(shù)據(jù)和地質(zhì)采礦條件，能夠很好的求取地表移動(dòng)變形預(yù)計(jì)參數(shù)。本文通過(guò)實(shí)例，對(duì)所擬合的參數(shù)進(jìn)行檢驗(yàn)，證明了此方法的可行性。

　　關(guān)鍵字：MATLAB,非線性最小二乘法,概率積分法,預(yù)計(jì)參數(shù)

　　1 引言

　　MATLAB有美國(guó)MathWorks公司推出，其與Mathematica、Maple并稱為三大數(shù)學(xué)軟件，具有優(yōu)秀的數(shù)值計(jì)算能力和卓越的數(shù)據(jù)可視化功能。 MATLAB自問(wèn)世以來(lái)，已經(jīng)發(fā)展稱為適合多學(xué)科、多種工作平臺(tái)的功能強(qiáng)大的大型軟件，廣泛應(yīng)用于大學(xué)科研、工程計(jì)算等領(lǐng)域，尤其在工程界。

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

　　2 概率積分法求參的數(shù)學(xué)模型

　　利用概率積分法關(guān)于任意點(diǎn)的地表下沉預(yù)計(jì)公式作為求取預(yù)計(jì)參數(shù)的基本數(shù)學(xué)模型，即：

　　地表移動(dòng)觀測(cè)站任何一個(gè)測(cè)點(diǎn)的下沉實(shí)測(cè)值W全能表達(dá)成自變量(x,y)(觀測(cè)點(diǎn)的水平坐標(biāo))和概率積分法預(yù)計(jì)參數(shù)q、tgβ、s3、s4、tgβ1、tgβ2、s1、s2、θ的函數(shù)，見(jiàn)下式：

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

　　3 Matlab曲線擬合求參函數(shù)

　　Matlab強(qiáng)大的繪圖和用戶自定義函數(shù)等功能，在數(shù)據(jù)處理領(lǐng)域被廣泛應(yīng)用，其中當(dāng)然包括對(duì)實(shí)驗(yàn)、實(shí)測(cè)數(shù)據(jù)的擬合工作。鑒于本文研究?jī)?nèi)容，著重介紹Matlab中非線性曲線擬合(非線性最小二乘擬合)函數(shù)。

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

　　已知擬合的數(shù)據(jù)點(diǎn)為：xdata和ydata，lsqcurvefit用以求含參量variable(向量)的向量值函數(shù)F(x,xdata)中的參變量x(向量)，使?jié)M足下式：

　　完整調(diào)用格式為：

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

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

　　4 地表移動(dòng)變形參數(shù)求取與檢驗(yàn)

　　4.1地表移動(dòng)變形預(yù)計(jì)參數(shù)求取的算法流程

　　4.2求取沉陷預(yù)計(jì)參數(shù)的程序設(shè)計(jì)

　　本文主要用到的兩個(gè)自定義函數(shù)budgets和horizonal_budgets，函數(shù)budgets用于求取沉陷預(yù)計(jì)參數(shù)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);

　　函數(shù)horizonal_budgets用于求取沉陷預(yù)計(jì)參數(shù)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);

　　特別注意的是在這兩個(gè)函數(shù)中都調(diào)用了LevenbergMarquardt(L-M)極小化方法，以避免出現(xiàn)局部解和病態(tài)問(wèn)題。

　　4.3 預(yù)計(jì)參數(shù)的實(shí)例驗(yàn)證

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

　　預(yù)計(jì)的走向觀測(cè)線下沉值和水平移動(dòng)值的中誤差分別為60.01mm和21.4mm，分別占最大值的2.6%和2.5%，預(yù)計(jì)的斜向觀測(cè)線下沉值和水平移動(dòng)值的中誤差分別為73.3mm和34.4mm，分別占最大值的3.1%和4.5%。

　　5 結(jié)論

　　通過(guò)預(yù)計(jì)的擬合曲線與實(shí)測(cè)數(shù)據(jù)的可視化比較，可知只要充分利用MATLAB工具箱函數(shù)，結(jié)合任意工作面監(jiān)測(cè)點(diǎn)實(shí)測(cè)數(shù)據(jù)及實(shí)際的地質(zhì)采礦條件，就能夠擬合出較為準(zhǔn)確的地表移動(dòng)變形預(yù)計(jì)參數(shù)，極為方便，但在此需特別指出，擬合時(shí)初始化的參數(shù)非常重要，因此要充分利用周圍礦區(qū)的資料。

　　參考文獻(xiàn)：

　　[1] 鄒友峰,鄧喀中,馬偉民.礦山開(kāi)采沉陷工程[M].徐州:中國(guó)礦業(yè)大學(xué)出版社,2003

　　[2] 吳侃,周鳴.礦山開(kāi)采沉陷預(yù)測(cè)預(yù)報(bào)系統(tǒng)[M].徐州:中國(guó)礦業(yè)大學(xué)出版社,2000

　　[3] 康建榮,王金莊,溫澤民.任意形多工作面多線段開(kāi)采沉陷預(yù)計(jì)系統(tǒng)(MSPS)[J]礦山測(cè)量,2000,1(3):24~27

　　[4] 朱劉娟,陳俊杰,鄒友峰.任意形狀工作面開(kāi)采地表移動(dòng)變形預(yù)計(jì)的算法實(shí)現(xiàn)[J].遼寧工程技術(shù)大學(xué)學(xué)報(bào),2005,24(3):337~339

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

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

　　[7] 周萬(wàn)茂,張華興,仲惟林.多工作面開(kāi)采影響下求取地表移動(dòng)參數(shù)方法研究.煤炭學(xué)報(bào).2000 25(3)

　　[8] 梁明,王成緒.厚黃土覆蓋山區(qū)開(kāi)采沉陷預(yù)計(jì).煤田地質(zhì)與勘探.2001 29(2)

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