劉吉波 王志紅 任傳建

(貴州工程應(yīng)用技術(shù)學院 礦業(yè)工程學院, 貴州 貴陽 551700)

0 引言

沉陷預(yù)計可以在地下資源采出前掌握地表破壞情況,在礦山資源開發(fā)利用中具有非常重要的作用。概率積分法是目前較為成熟,應(yīng)用廣泛的預(yù)計方法,已有很多預(yù)計軟件,多是采用Visual C++、C#、Visual Basic.NET等軟件平臺開發(fā)[1-6],也有采用Matlab、AutoCAD、ArcGIS等進行二次開發(fā)實現(xiàn)的[7-10]。

Python是一種跨平臺的計算機程序設(shè)計語言,是結(jié)合了解釋性、編譯性、互動性和面向?qū)ο蟮哪_本語言,具有免費、可移植、功能強大、易于使用的特點。Python除了math等標準程序庫外,還提供了大量成熟的專業(yè)程序包,其中numpy和scipy主要用于科學計算,sympy程序庫具有符號計算功能[11-15]。

本研究利用scipy模塊的積分功能實現(xiàn)了開采沉陷預(yù)計計算,對算法進行了優(yōu)化設(shè)計,并采用matplotlib模塊實現(xiàn)了等值線繪制。

1 利用Python進行概率積分預(yù)計的基本算法

1.1 地表任意點移動變形概率積分計算公式[16-17]

概率積分法地表下沉值計算公式為:

(1)

式(1)中,Wmax為地表充分采動時的最大下沉值;r為主要影響半徑;D為開采區(qū)域;X、Y為預(yù)計點地面坐標。傾斜、水平移動、曲率等變形可依據(jù)下沉公式計算而得。

1.2 scipy積分計算功能及精度分析

scipy模塊提供了豐富的積分運算,其中integrate子模塊包含一重、二重及三重積分函數(shù),二重積分可以用dblquad函數(shù)計算。dblquad的一般調(diào)用形式是dblquad(func,a,b, gfun, hfun),其中,func是待積分函數(shù)名,b、a是x變量的上下限,hfun、gfun為定義y變量上下限的函數(shù)名。dblquad函數(shù)的返回值是一個tuple類型的變量(result, abserr),result是積分結(jié)果,abserr是積分誤差。

import scipy as sp

def func(y,x):

return sp.exp(-x-y)

gfun=lambda x:0.0

hfun=lambda x:x

(a,b)=(0.0,2.0)

value=sp.integrate.dblquad(func,a,b,gfun,hfun)

print(“%.3f”% value[0])

上述代碼輸出為0.374，計算誤差小于4.0×10-15，scipy的dblquad函數(shù)積分精度很高，滿足沉陷預(yù)計要求。

1.3 利用scipy實現(xiàn)概率積分預(yù)計

概率積分法適用疊加原理,故對開采區(qū)域進行相應(yīng)劃分后,可使每個子開采區(qū)域符合積分計算要求。因開采區(qū)域一般是由多條線段構(gòu)成的多邊形,故邊界可用直線方程表示。為簡化問題,以單一區(qū)域開采下沉預(yù)計為例進行研究。

設(shè)有一水平工作面,煤厚m=3.0 m,下沉系數(shù)q=0.7,開采深度H=200.0 m,主要影響角正切tanβ=2.0。X方向位于區(qū)間[300 m,1000 m],Y方向位于直線y1=0.1x+270和y2=0.1x+470之間,地面預(yù)計格網(wǎng)左下角坐標為(0,0),網(wǎng)格間距20 m,X方向65個網(wǎng)格,Y方向40個網(wǎng)格。

為確保程序通用性,上述工作面開采地表下沉計算的程序?qū)崿F(xiàn)為:

import scipy as sp

(m,q,tanb,H)=(3.0,0.7,2.0,200.0)

(W0,r)=(m*q*1000,H/tanb)

(x1,x2)=(300.0,1000.0)

y1=lambda x:a1*x+b1

y2=lambda x:a2*x+b2

(a1,b1,a2,b2)=(0.1,270,0.1,470)

(nx,ny,dx,dy)=(65,40,20,20)

(X0,Y0)=(0,0)

(Xn,Ym)=(X0+(nx+1)*dx,

Y0+(ny+1)*dy)

data=[]

for X in range(X0,Xn,dx):

for Y in range(Y0,Ym,dy):

def pintegral(y,x):

return sp.exp(-sp.pi*(((x-X)/r)

**2+((y-Y)/r)**2))/(r*r)

c=sp.integrate.dblquad(pintegral,

x1,x2,y1,y2)

data.append([X,Y,W0*c[0]])

根據(jù)預(yù)計結(jié)果,繪制下沉等值線圖。

1.4 預(yù)計點密度選擇

利用概率積分法進行地表移動變形預(yù)計時,預(yù)計點的密度對于預(yù)計效果和運行效率至關(guān)重要。密度過小,雖用時短,但精度會大幅降低,導致變形等值線呈現(xiàn)劇烈的鋸齒狀;密度過大,精度高,但預(yù)計耗時多,預(yù)計工作面多,預(yù)計區(qū)域范圍大時影響更大。根據(jù)試驗,一般情況下網(wǎng)格密度50 m×50 m時即滿足需要,當要求較高時可適當增加網(wǎng)格密度,反之網(wǎng)格密度可降低。

2 預(yù)計算法優(yōu)化

2.1 積分函數(shù)構(gòu)造

Python的符號計算模塊sympy、科學計算模塊scipy、數(shù)值計算模塊numpy和數(shù)學函數(shù)模塊math中都定義了指數(shù)函數(shù)exp。在其他條件不變情況下,選用不同模塊的exp函數(shù)構(gòu)造積分函數(shù),程序運行時間有很大差異,如表1所示。從表1可以看出采用math模塊的exp比scipy的計算效率提高1倍,numpy和scipy的效率相當,而sympy因要進行符號運算,效率極低,不建議采用。

表1 不同exp計算時間(循環(huán)內(nèi)部)

將預(yù)計點X、Y坐標定義為全局變量,積分函數(shù)定義于循環(huán)體之外,程序運行時間如表2所示。通過表2和表1可知,程序性能提高1倍左右。

表2 不同exp計算時間(循環(huán)外部)

2.2 利用通用函數(shù)提高效率

numpy模塊提供了通用函數(shù)功能,其具有與輸入數(shù)組形狀相同的輸出數(shù)組,可以一次性對所有數(shù)組數(shù)據(jù)進行計算[12-13],避免了循環(huán)操作,從而提高程序計算效率。通過vectorize函數(shù)可以快速創(chuàng)建通用函數(shù)。

#定義通用函數(shù)

def cal(n):

(row,col)=(int(n/(nx+1)),n-row*(nx+1))

(X,Y)=(X0+col*dx,Y0+row*dy)

def pintegral (y,x):

return mh.exp(-mh.pi*(((x-X)/r)**2

+((y-Y)/r)**2))/r/r

c=dblquad(pintegral,x1,x2,y1,y2)

return round(c[0]*W0,0)

#預(yù)計點數(shù)組

points_array=np.arange((ny+1)*(nx+1))

#生成通用函數(shù)

vgauss=np.vectorize(cal)

#預(yù)計計算

result=vgauss(points_array)

由表3和表1可知,通用函數(shù)操作可極大提高程序性能,減少運行時間。

表3 調(diào)用通用函數(shù)運行時間

2.3 利用影響圓確定計算范圍

在水平煤層開采時,對地表點A有影響的煤層開采范圍是一個圓,其以A在煤層的垂直投影點O為圓心,半徑R=Hctanδ0,δ0為邊界角,如圖1所示。只有位于圓內(nèi)的煤層開采對A點有影響,圓外的煤層開采對A點影響為0。當多工作面開采尤其是土地復(fù)墾等需要進行全井田預(yù)計時,預(yù)計范圍大,計算點數(shù)多,利用影響圓法可以大幅提高程序效率[14-16]。

圖1 水平煤層開采對地表點的影響圓

設(shè)計工作面如圖2所示,取邊界角δ0=55°,則R=140.0 m。將工作面邊界向外側(cè)偏移R,得到新邊界C和D,位于邊界C和D內(nèi)的點進行預(yù)計計算,位于C和D外的點移動變形直接賦值為0。

圖2 有效預(yù)計范圍確定

全部預(yù)計,共10 201個預(yù)計點,用時45 s,利用影響圓法,共2 652個預(yù)計點,用時6 s。

實際預(yù)計時,對于非水平煤層開采,可分別計算工作面走向方向、上山方向和下山方向的影響圓半徑,并確定有效預(yù)計范圍。

2.4 坐標系選擇

在實際開采中,工作面走向經(jīng)常是任意方向的,造成沉陷預(yù)計復(fù)雜性提高。將y方向預(yù)計網(wǎng)格數(shù)改為ny=75,選用math.exp函數(shù),當Y積分限斜率a=0.0,即工作面走向沿X方向時,預(yù)計用時16 s;當Y積分限斜率a=1.0時,預(yù)計用時20 s,說明工作面走向與Y軸的不垂直度增加,程序用時增加。同時,對于矩形工作面,工作面需進行分割才能進行積分運算。

選擇合適的工作面坐標系進行積分計算可以解決此問題。如圖3所示,當工作面為矩形時,可以工作面左下角點為原點,X軸沿工作面走向方向;當工作面為非矩形時,可以令X軸沿長對角線方向;當工作面為多邊形時,選最長的兩點連線為X軸。此時需要將計算點坐標轉(zhuǎn)化至工作面坐標系后進行積分計算。從地面坐標系到工作面坐標系轉(zhuǎn)換公式為[17]:

圖3 坐標系選擇

(2)

其中,φ為工作面坐標系x軸順時針與大地坐標系X軸的夾角,(X0,Y0)為工作面坐標系原點O的大地坐標系下的坐標。

3 等值線繪制

matplotlib是Python實用的圖形和表格繪制軟件包,可用于便捷地繪制等值線圖、等值線云圖和三維曲面圖[15]。

繪制二維等值線圖的部分代碼如下,生成等值線圖如圖4所示。

圖4 下沉等值線圖

#datax,datay,datav分別存放X坐標、Y坐標和變形值

#建立三角剖分

triang=tri.Triangulation(datax,datay)

#實現(xiàn)等值線繪制

con=plt.tricontour(datax,datay,triang.triangles,datav,levels=listlevels,cmap='rainbow')

#標注等值線

label=plt.clabel(con, inline=False, fmt='%.0f', fontsize=20)

繪制三維下沉曲面和下沉等值線云圖的部分代碼如下,生成下沉曲面和等值線云圖如圖5所示。

圖5 下沉曲面圖和下沉等值線云圖

#添加繪圖子窗口

ax=fig.add_subplot(111, rojection='3d')

# 設(shè)置圖像z軸的顯示范圍

ax.set_zlim(3000,0)

#繪制下沉曲面

ax.plot_surface(X3d,Y3d,Z3d,linewidths=1, rstride=1,cstride=1, cmap=plt.get_cmap('rainbow'))

#繪制下沉等值線云圖

cset=ax.contourf(X3d,Y3d,Z3d, zdir='z',

levels=vlevels,linewidths=1,offset=3000, cmap='rainbow')

4 結(jié)束語

對利用scipy的積分模塊進行開采沉陷預(yù)計計算和利用matplotlib繪制等值線圖及曲面圖進行了研究,得出以下結(jié)論:

(1)scipy科學計算模塊功能強大、高效,用少量代碼即可完成復(fù)雜的概率積分預(yù)計;利用math.exp比numpy.exp、sicpy.exp及sympy.exp計算效率高;將積分函數(shù)定義于循環(huán)體外,可以進一步提升程序性能;利用通用函數(shù)進行數(shù)組整體操作可提高計算效率。

(2)多工作面大范圍預(yù)計時,利用影響圓法可以大幅減少參與計算的點數(shù);合理選擇坐標系、設(shè)計預(yù)計點密度也能提高程序性能。

(3)利用matplotlib繪圖模塊,可以快捷繪制變形等值線圖、等值線云圖或三維曲面圖,圖形美觀實用。