劉吉波 王志紅 任傳建

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

0 引言

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

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

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

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

1.1 地表任意點(diǎn)移動(dòng)變形概率積分計(jì)算公式[16-17]

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

(1)

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

1.2 scipy積分計(jì)算功能及精度分析

scipy模塊提供了豐富的積分運(yùn)算,其中integrate子模塊包含一重、二重及三重積分函數(shù),二重積分可以用dblquad函數(shù)計(jì)算。dblquad的一般調(diào)用形式是dblquad(func,a,b, gfun, hfun),其中,func是待積分函數(shù)名,b、a是x變量的上下限,hfun、gfun為定義y變量上下限的函數(shù)名。dblquad函數(shù)的返回值是一個(gè)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，計(jì)算誤差小于4.0×10-15，scipy的dblquad函數(shù)積分精度很高，滿足沉陷預(yù)計(jì)要求。

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

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

設(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ù)計(jì)格網(wǎng)左下角坐標(biāo)為(0,0),網(wǎng)格間距20 m,X方向65個(gè)網(wǎng)格,Y方向40個(gè)網(wǎng)格。

為確保程序通用性,上述工作面開采地表下沉計(jì)算的程序?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ù)計(jì)結(jié)果,繪制下沉等值線圖。

1.4 預(yù)計(jì)點(diǎn)密度選擇

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

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

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

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

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

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

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

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

numpy模塊提供了通用函數(shù)功能,其具有與輸入數(shù)組形狀相同的輸出數(shù)組,可以一次性對(duì)所有數(shù)組數(shù)據(jù)進(jìn)行計(jì)算[12-13],避免了循環(huán)操作,從而提高程序計(jì)算效率。通過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ù)計(jì)點(diǎn)數(shù)組

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

#生成通用函數(shù)

vgauss=np.vectorize(cal)

#預(yù)計(jì)計(jì)算

result=vgauss(points_array)

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

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

2.3 利用影響圓確定計(jì)算范圍

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

圖1 水平煤層開采對(duì)地表點(diǎn)的影響圓

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

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

全部預(yù)計(jì),共10 201個(gè)預(yù)計(jì)點(diǎn),用時(shí)45 s,利用影響圓法,共2 652個(gè)預(yù)計(jì)點(diǎn),用時(shí)6 s。

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

2.4 坐標(biāo)系選擇

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

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

圖3 坐標(biāo)系選擇

(2)

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

3 等值線繪制

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

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

圖4 下沉等值線圖

#datax,datay,datav分別存放X坐標(biāo)、Y坐標(biāo)和變形值

#建立三角剖分

triang=tri.Triangulation(datax,datay)

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

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

#標(biāo)注等值線

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é)束語

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

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

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

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