蔡音飛，VERDEL Thierry，OLIVIER Deck，李曉靜

(1.太原理工大學(xué) 礦業(yè)工程學(xué)院，山西 太原　030024；2.GeoRessources,Ecole des Mines de Nancy,Université de Lorraine,Nancy　54000；3.山西農(nóng)業(yè)大學(xué) 資源環(huán)境學(xué)院，山西 太谷　030800)

?

地形影響下的開采沉陷影響函數(shù)法優(yōu)化

蔡音飛1,2，VERDEL Thierry2，OLIVIER Deck2，李曉靜3

(1.太原理工大學(xué) 礦業(yè)工程學(xué)院，山西 太原030024；2.GeoRessources,Ecole des Mines de Nancy,Université de Lorraine,Nancy54000；3.山西農(nóng)業(yè)大學(xué) 資源環(huán)境學(xué)院，山西 太谷030800)

摘要:影響函數(shù)法廣泛應(yīng)用于地下層狀礦體的開采沉陷計算。理論上，該方法適用于水平地表條件下的沉陷預(yù)計。將地形變化納入影響函數(shù)法，使此方法可以應(yīng)用于非水平地表條件下。為了剔除其他開采和地質(zhì)因素的影響，除地形變化外可能影響沉陷形態(tài)的因素都被固定了的簡化數(shù)值模型被應(yīng)用于研究中。根據(jù)這些數(shù)值模擬實驗的結(jié)果，若干影響函數(shù)法的參數(shù)，包括影響半徑、影響角和下沉率，通過4個地形相關(guān)的修正系數(shù)被重定義為地表相對礦體高程的函數(shù)。優(yōu)化后的影響函數(shù)法可以更好的擬合非水平地表條件下的數(shù)值模擬結(jié)果和現(xiàn)場觀測數(shù)據(jù)。相較于經(jīng)典影響函數(shù)法，該方法需要將地表點相對礦體的高程作為輸入數(shù)據(jù)。在其他開采、地質(zhì)條件下，構(gòu)建簡單數(shù)值模型的方法可以被用于估算相應(yīng)條件下的地形影響修正系數(shù)。

關(guān)鍵詞:影響函數(shù)法；數(shù)值模擬；影響半徑；影響角；下沉率

地下采礦引起的地表沉陷對采區(qū)的土地、建筑、道路等的安全和功能造成了不同程度的影響。隨著可持續(xù)發(fā)展的這一科學(xué)發(fā)展觀在國民生產(chǎn)生活中日益得到重視，開采沉陷的研究在礦山生產(chǎn)活動中也顯得愈加重要。影響函數(shù)法[1-3](含概率積分法)廣泛應(yīng)用于地下層狀礦體的開采沉陷預(yù)計。國內(nèi)外有很多專家學(xué)者對非水平地表條件下的開采沉陷進行了研究，分析了沉陷規(guī)律[4-8]，并對影響函數(shù)法在受地形影響的開采沉陷計算中的應(yīng)用提出了一些優(yōu)化方法[8-11]。這些方法多考慮地質(zhì)條件與地形的綜合影響，導(dǎo)出的公式和參數(shù)的選取也相對比較復(fù)雜。本文利用數(shù)值模擬，考慮近水平礦體開采，并將地形對沉陷的影響單獨分離出來，對影響函數(shù)法進行優(yōu)化，改進沉陷計算的結(jié)果。相對于經(jīng)典影響函數(shù)法，本方法需將地形作為一個已知輸入數(shù)據(jù)。不同地質(zhì)條件下的地形影響修正參數(shù)的取值可以依據(jù)該條件下的數(shù)值模擬結(jié)果(沒有已知參照的情況下)，一定程度上避免了人工判斷的不確定性。

1經(jīng)典影響函數(shù)法

經(jīng)典的影響函數(shù)法是基于水平地表的，用于計算該條件下地表點的下沉(豎向移動)值。其他的地表移動變形值，包括傾斜、曲率、水平移動、水平變形，都可以通過數(shù)學(xué)方法由下沉計算。

很多數(shù)學(xué)函數(shù)都可以被用作影響函數(shù)，大部分是指數(shù)函數(shù)[2-3]，其中最常用的2個函數(shù)的二維表達式列于表1[10]：Knothe等所用的是標(biāo)準(zhǔn)正態(tài)分布的概率密度函數(shù)，Litwiniszyn等所用的是該函數(shù)的變形形式。本質(zhì)上，一個影響函數(shù)描述的是單元開采造成的任一地表點的下沉值，并可由此函數(shù)確定單元下沉盆地。某一地表點處由全采區(qū)開采造成的下沉是該點處所有單元下沉(由影響函數(shù)決定)的累積，見式(1)；由于其最大值永遠為1(由函數(shù)性質(zhì)決定)，需乘上一個系數(shù)mq(水平礦體)來調(diào)節(jié)大小。按此方法可求得所有地表點由全采區(qū)開采造成的下沉值。

表1　沉陷計算中常用的影響函數(shù)

注：x為地表點橫坐標(biāo)；s為開采單元橫坐標(biāo)；x-s為地表點到開采單元的水平距離；R為單元開采的影響半徑，當(dāng)|x-s| >R時，f(x)極小可忽略；t∈[1,2]，由地質(zhì)條件決定。

由表1中所列的影響函數(shù)的數(shù)學(xué)性質(zhì)可知，理論上，單元開采的下沉曲線是關(guān)于通過該單元的豎軸對稱的，水平礦體全采區(qū)開采的下沉曲線也是關(guān)于通過開采中心的豎軸對稱的。

(1)

式中，W為地表點由全采區(qū)(包括n個開采單元E1～En)開采引起的下沉值；m為礦體開采厚度；q為下沉率；f(x)為表1所列的任一影響函數(shù)。

2地形變化對沉陷影響的數(shù)值模擬分析

實際礦山生產(chǎn)活動中，地表沉陷受到很多因素的影響，主要包括開采條件、覆巖條件、地形條件等，難以將地形變化對沉陷的影響單獨分離出來，進行研究。所以本文擬使用FLAC3D數(shù)值模擬軟件，應(yīng)用一系列簡化的數(shù)值模型來分析僅地形變化對開采沉陷造成的影響。

FLAC3D是由Itasca公司為巖土工程應(yīng)用開發(fā)的連續(xù)介質(zhì)顯式有限差分軟件。該軟件能較好的模擬巖土體材料的力學(xué)行為及材料達到強度或屈服極限時產(chǎn)生的破壞和塑性流動，可應(yīng)用于大變形情況的模擬，如本文涉及的開挖工程。軟件計算所采用的數(shù)學(xué)模型是根據(jù)彈塑性理論的基本原理(應(yīng)變定義、運動定律、能量守衡定律、平衡方程及理想材料的連續(xù)性方程等)建立的[12-13]。本文將采用莫爾-庫侖彈塑性材料本構(gòu)模型進行計算。FLAC3D(包括FLAC)已被廣泛應(yīng)用于各種條件下的沉陷研究，如“三下”開采[13]、急傾斜礦體開采[14]、條帶開采[15]、非充分開采[16]等。

(1)模型分層。

研究中，所用模型均由3個地層組成，自下而上分別為底板層、礦體層和覆巖層，其地表傾向主斷面如圖1所示。

圖1　模型分層和尺寸(地表傾向主斷面)Fig.1　Strata and dimensions of the models in the middle vertical section along the dip direction of the surface

(2)模型尺寸。

標(biāo)準(zhǔn)模型的平面尺寸為1 000 m×1 000 m，當(dāng)傾角增大時，地表傾向上山方向的長度會適當(dāng)延長以保證沉陷邊界在模型范圍內(nèi)。模型底板厚度為100 m；礦體厚度為5 m，呈水平層狀，開采范圍300 m×300 m并位于標(biāo)準(zhǔn)模型正中；覆巖平均厚度D(即礦體平均采深，可由中心軸處地表到礦體豎直距離表示)可變，為100～350 m，上表面傾角α(即地表傾角)可變，為0°～20°。

(3)模型物理力學(xué)性質(zhì)。

模型使用摩爾-庫倫模型，各分層的物理、力學(xué)性質(zhì)均固定不變，見表2。覆巖層被簡化為一層連續(xù)、巖性相同的巖層，其物理、力學(xué)性質(zhì)可用原地層的層厚為權(quán)重，采用加權(quán)平均法求取。模型的左右邊界的橫向移動和底部的豎向移動被限制，上表面為自由面，所有內(nèi)部移動均由自重引起，無附加外力。

表2　模型的物理和力學(xué)性質(zhì)

簡言之，模型中除了地形因素(覆巖的平均厚度和上表面傾角)，其他可影響沉陷形態(tài)的因素都是固定不變的。地表傾向上山方向的長度增加并不影響沉陷形態(tài)。

本文運行了2組數(shù)值模擬實驗，每組為20個模型，2組間大部分模型可公用。第1組模型，固定平均采深，改變地表傾角；第2組模型，固定地表傾角，改變平均采深。根據(jù)結(jié)果，地表傾角增加使得下沉和水平移動的數(shù)值輕微變大，其形態(tài)由原來的對稱變?yōu)槊黠@的非對稱，表現(xiàn)為地表上、下山方向的影響范圍不對等和極值位置偏移；平均采深變大使得得下沉和水平移動的數(shù)值變小。

3基于數(shù)值模擬優(yōu)化影響函數(shù)法

為了將影響函數(shù)法用于非水平地表條件下的開采沉陷計算，本節(jié)依據(jù)上述數(shù)值模擬在地表傾向主斷面上的結(jié)果，對影響函數(shù)進行優(yōu)化。

(1)影響半徑R。

表1所列的經(jīng)典影響函數(shù)中，單元開采的影響半徑R是常數(shù)(圖2中R0)。對于傾斜地表，各地表點的影響半徑卻是不同的[11]：例如圖2中，為了計算點Pa由單元Ex開采引起的下沉，其影響半徑應(yīng)取Ra；對于點Pb，其影響半徑應(yīng)取Rb。假設(shè)任意地表點相對開采單元的高程已知，為H(x)，則任意點的影響半徑可由式(2)計算。

(2)

式中，φ為影響角，定義為開采單元和地表影響邊界的連線到水平線的夾角。

圖2　變化的影響半徑Fig.2　Varying influence radius for one mining element

(2)影響角φ。

由式(2)可知，非水平地表的影響半徑計算時需要考慮影響角。根據(jù)對數(shù)值模擬結(jié)果的分析，發(fā)現(xiàn)影響角隨開采深度變化不大，其隨地表傾角的變化如圖3所示。隨著地表傾角的增大，地表上山方向的影響角增大，而下山方向的影響角減小，地表上、下山方向的影響角不再相等，其增減幅度與地表傾角成線性關(guān)系，見式(3)。由數(shù)值模擬結(jié)果求得k，c分別為0.27，47.5°。

(3)

式中，φu，φd分別為地表上、下山方向的影響角；k，c為地形變化對影響角的修正系數(shù)。

圖3　影響角隨地表傾角變化規(guī)律Fig.3　Variation of the influence angles with the surface dip angle

對于起伏地表，任意地表點的傾角為H(x)導(dǎo)數(shù)的反正切函數(shù)，見式(4)，將其代入式(3)，計算起伏地表任意點沉陷時所用的地表上下山影響角應(yīng)如式(5)所示。

(4)

(5)

將式(5)代入式(2)，起伏地表上任意點的影響半徑計算公式為分段函數(shù)，見式(6)。以圖2為例，計算開采單元Ex引起的地表沉陷時，如果地表點在Ex左側(cè)(xs)，應(yīng)使用φu計算影響半徑；規(guī)定Ex正上方的地表點，影響半徑為由φd，φu計算得到的影響半徑的均值。

(6)

特別地，水平地表條件下，任意點的α恒為0，H為常數(shù)，所以上下山φ相等，R亦相等。

(3)下沉率q。

上文已述，最終下沉的數(shù)值需要用下沉率來調(diào)整，見式(1)。根據(jù)對數(shù)值模擬結(jié)果的分析，發(fā)現(xiàn)下沉率隨地表傾角變化不大，其隨平均開采深度的變化可見圖4，成反比例關(guān)系。將該關(guān)系應(yīng)用于地表上任意點，得到下沉率與地表點高程(相對于礦體)的關(guān)系，見式(7)。由數(shù)值模擬結(jié)果求得a，b分別為0.21，44.8。特別地，水平地表條件下，H為常數(shù)，q也為定值。

(7)

式中，a，b為地形變化對下沉率的修正系數(shù)。

圖4　下沉率隨平均開采深度變化規(guī)律Fig.4　Variation of the vertical subsidence coefficient with the mean mining depth

(4)影響函數(shù)法的優(yōu)化和討論。

綜上所述，將式(6)引入表1所列的影響函數(shù)(用R(x)代替R)，即可得到針對地形變化優(yōu)化的影響函數(shù)。新函數(shù)為分段函數(shù)，考慮了地表的變化，包括每個地表點相對礦體的高程和傾角(傾角可由高程函數(shù)求導(dǎo)后獲得)。將優(yōu)化的影響函數(shù)和式(7)代入式(1)，即可得改進的地表任意點受全采區(qū)開采影響的下沉值計算公式。其他地表移動變形值的計算沿用經(jīng)典方法。

上述影響函數(shù)法優(yōu)化算法已在MathematicaTM中編程實現(xiàn)。相較于經(jīng)典影響函數(shù)法，優(yōu)化方法需要將地表點相對礦體的高程函數(shù)作為一個已知數(shù)據(jù)輸入計算。由于地形的影響是本文的研究對象，所以地形函數(shù)已知也是研究的前提。

本文引入的4個地形相關(guān)的修正系數(shù)(a，b，c，k)，在不同實例中應(yīng)用時，需要按照文中方法用數(shù)值模擬實驗進行驗證。由于本文分離了地形的影響，也就是說其他開采、地質(zhì)條件均為固定不變的，文中提供的這些系數(shù)的數(shù)值并不適用于所有條件。但應(yīng)用簡化的數(shù)值模型求這幾個修正系數(shù)的方法是通用的。實際中，礦體的上覆巖層往往多于一層，可用他們的物理、力學(xué)性質(zhì)的加權(quán)平均(層厚為權(quán)重)作為簡化的數(shù)值模型的覆巖性質(zhì)。

4結(jié)果對比

對于非水平地表，使用變化的影響半徑R是必須的，否則無法體現(xiàn)地表高程的變化，計算結(jié)果也必然成對稱形態(tài)。所以此處將優(yōu)化的影響函數(shù)法計算結(jié)果(僅調(diào)整R或同時調(diào)整R，q，φ)與數(shù)值模擬結(jié)果及現(xiàn)場觀測數(shù)據(jù)進行對比。

(1)與數(shù)值模擬結(jié)果比較。

本例中數(shù)值模型的基本形態(tài)如圖1所示，地表傾角為15°，平均采深250 m，礦體開采長度300 m。數(shù)值模擬結(jié)果與優(yōu)化的影響函數(shù)法計算結(jié)果(僅調(diào)整R或同時調(diào)整R，q，φ)對比如圖5所示，并統(tǒng)計幾者差值的平方和列于表3。

圖5　數(shù)值模擬結(jié)果與優(yōu)化的影響函數(shù)法結(jié)果對比Fig.5　Comparison between the subsidence curves obtained from numerical simulation and improved influence function method

比較對象沉陷差值的平方和與優(yōu)化的影響函數(shù)法(僅R)比較與優(yōu)化的影響函數(shù)法(R,q,φ)比較數(shù)值模擬下沉0.5060.255水平移動0.5110.253觀測數(shù)據(jù)下沉2.1080.955水平移動2.0210.796

(2)與現(xiàn)場觀測數(shù)據(jù)[17]比較。

本例中現(xiàn)場觀測數(shù)據(jù)為晉城煤業(yè)集團成莊煤礦2307綜放工作面傾向觀測線的開采沉陷數(shù)據(jù)。該工作面采用走向長壁綜采放頂煤采煤法，走向1 845 m，寬度147.5 m。開采煤層為山西組3號近水平煤層，平均厚度6.7 m，平均開采深度327.8 m。工作面上方地面標(biāo)高+875～+950 m?，F(xiàn)場觀測數(shù)據(jù)與優(yōu)化的影響函數(shù)法計算結(jié)果(僅調(diào)整R或同時調(diào)整R,q,φ)對比如圖6所示，并統(tǒng)計幾者差值的平方和列于表3。

圖6　現(xiàn)場觀測數(shù)據(jù)與優(yōu)化的影響函數(shù)法結(jié)果對比Fig.6　Comparison between the subsidence curves obtained from field measurement and improved influence function method

由圖5，6和表3可見，通過對R，q，φ的調(diào)整，影響函數(shù)法可以更好的擬合數(shù)值模擬結(jié)果和現(xiàn)場觀測數(shù)據(jù)。

5結(jié)論

(1)提出了一種優(yōu)化經(jīng)典影響函數(shù)法的方法，使其可以將變化的地表納入計算。相對于經(jīng)典方法，改進的方法需考慮地表點相對礦體的高程及其4個地形相關(guān)的修正系數(shù)。優(yōu)化后的影響函數(shù)為分段函數(shù)，可以更好的擬合非水平地形影響下的沉陷曲線。

(2)對數(shù)值模型進行了簡化，使其僅反映地形的變化。該簡化方法具有通用性，可應(yīng)用于其他開采、地質(zhì)條件，估算相應(yīng)條件下的地形影響修正系數(shù)。

(3)應(yīng)用數(shù)值模擬結(jié)果，修正了若干影響函數(shù)法的參數(shù)，包括影響半徑、影響角和下沉率，他們均可以表示為地表相對高程的函數(shù)。

參考文獻:

[1]National Coal Board.Subsidence engineers’ handbook[M].London:National Coal Board,1975.

[2]Whittaker N Barry,Reddish J David.Subsidence occurrence,prediction and control[M].Amsterdam,Oxford,New York,Tokyo:Elsevier Science Publishers,1989.

[3]何國清，楊倫，凌賡娣,等.礦山開采沉陷學(xué)[M].徐州：中國礦業(yè)大學(xué)出版社,1991.

[4]Holla L.Ground movement due to longwall mining in high relief areas in New South Wales,Australia[J].International Journal of Rock Mechanics and Mining Sciences,1997,34(4):775-787.

[5]戴華陽，翟厥成，胡友健.山區(qū)地表移動的相似模擬實驗研究[J].巖石力學(xué)與工程學(xué),2000,19(4):1027-1033.

Dai Huayang,Zhai Juecheng,Hu Youjian.Testing study on surface displacement of mountainous region with similar material[J].Chinese Journal of Rock Mechanics and Engineering,2000,19(4):1027-1033.

[6]黃珍珍，吳侃，徐靖.山區(qū)地表移動變形規(guī)律分析[J].煤炭科技,2007(1):26-29.

Huang Zhenzhen,Wu Kan,Xu Jing.Rules of epidermal deformation in mountainous areas[J].Coal Science and Technology Magazine,2007(1):26-29.

[7]Dai Huayang,Lian Xugang,Liu Jiyan,et al.Model study of deformation induced by fully mechanized caving below a thick loess layer[J].International Journal of Rock Mechanics and Mining Sciences,2010,47:1027-1033.

[8]何萬龍.山區(qū)開采沉陷與采動損害[M].北京：中國科技出版社，2003.

[9]戴華陽.負(fù)指數(shù)函數(shù)法預(yù)計山區(qū)地表移動[J].礦山測量,1990(3):48-51.

Dai Huayang.Mountains mining subsidence prediction by negative exponent method[J].Mining Surveying,1990(3):48-51.

[10]Cai Yinfei,Verdel Thierry,Deck Olivier.On the topography influence on subsidence due to horizontal underground mining using the influence function method[J].Computers and Geotechnics,2014,61:328-340.

[11]Saeidi Ali.The vulnerability of structures subjected to hazards caused by ground movements:Development of a simulator for damages[D].Nancy:Ecole des Mines,2010.

[12]謝和平，周宏偉，王金安,等.FLAC在煤礦開采沉陷預(yù)測中的應(yīng)用及對比分析[J].巖石力學(xué)與工程學(xué)報,1999,18(4):397-401.

Xie Heping,Zhou Hongwei,Wang Jin’an,et al.Application on FLAC to predict ground surface displacements due to coal extraction and its comparative analysis[J].Chinese Journal of Rock Mechanics and Engineering,1999,18(4):397-401.

[13]王生俊，賈學(xué)民，韓文峰,等.高速公路下伏采空區(qū)剩余沉降量 FLAC3D計算方法[J].巖石力學(xué)與工程學(xué)報,2005,24(19):3545-3550.

Wang Shengjun,Jia Xuemin,Han Wenfeng,et al.A FLAC3Dmethod for calculation of residual subsidence in a mining area beneath an expressway[J].Chinese Journal of Rock Mechanics and Engineering,2005,24(19):3545-3550.

[14]高明中，余忠林.急傾斜煤層開采對地表沉陷影響的數(shù)值模擬[J].煤炭學(xué)報,2004,28(6):578-582.

Gao Mingzhong,Yu Zhonglin.Numerical analysis of surface subsidence behavior relative to steep excavation[J].Journal of China Coal Society,2004,28(6):578-582.

[15]袁堅，馮濤，劉金海,等.條帶采寬及留寬對地表沉陷影響的研究[J].采礦與安全工程學(xué)報,2007,24(1):88-91.

Yuan Jian,Feng Tao,Liu Jinhai,et al.Effect of mining width and retaining width on ground subsidence[J].Journal of Mining & Safety Engineering,2007,24(1):88-91.

[16]郭俊廷，閻躍觀，戴華陽,等.基于FLAC3D的非充分采動地表移動求參方法[J].礦業(yè)工程研究,2010,25(4):43-46.

Guo Junting,Yan Yueguan,Dai Huayang,et al.Determination of ground movement permanents under insufficient mining condition based on simulation results of FLAC3D[J].Mineral Engineering Research,2010,25(4):43-46.

[17]宋新華，鄧喀中，閻躍觀,等.地形對地表移動變形的影響分析[J].礦山測量,2007(3):70-72.

Song Xinhua,Deng Kazhong,Yan Yueguan,et al.Analysis of topography influence on ground displacement and deformation[J].Mining Surveying,2007(3):70-72.

Improving the influence function method to take topography into the calculation of mining subsidence

CAI Yin-fei1,2，VERDEL Thierry2，OLIVIER Deck2，LI Xiao-jing3

(1.CollegeofMiningEngineering,TaiyuanUniversityofTechnology,Taiyuan030024,China;2.GeoRessources,EcoledesMinesdeNancy,UniversitedeLorraine,Nancy54000,France；3.CollegeofResourcesandEnvironment,ShanxiAgriculturalUniversity,Taigu030800,China)

Abstract:The classic influence function method is often used in the calculation of mining subsidence caused by stratiform underground excavations.Theoretically,its use is limited to the subsidence predictions under the condition of horizontal ground surface.In order to improve the original influence function method to take topographic variations into account.Due to real-world mining conditions that are usually complicated,it is difficult to separate topography influences from influences caused by other factors.The present work uses simplified numerical simulation models to analyze the topography influence only,avoiding the effects of other factors.Then several parameters of the influence function method,including the influence radius,influence angle,and vertical subsidence coefficient,are adjusted by considering four correction factors with respect to the topography.They all can be defined as functions of the relative surface elevations with ore body.After this,the improved influence function method can well fit the numerical simulation results and field measurement data.Comparing to the original method,the improved method only takes the topography as additional input data.The method of constructing the simplified numerical simulation models can be used to solve the topographic correction factors under other mining and geological conditions.

Key words:influence function method;numerical simulation;influence radius;influence angle;vertical subsidence coefficient

中圖分類號:TD325.2

文獻標(biāo)志碼:A

文章編號:0253-9993(2016)01-0271-06

作者簡介:蔡音飛(1983—),男,浙江嘉興人,講師，博士。E-mail:caiyinfeimines@163.com,yinfei.cai@univ-lorraine.fr

收稿日期:2014-12-15修回日期:2015-01-19責(zé)任編輯:常琛

蔡音飛，VERDEL Thierry，OLIVIER Deck，等.地形影響下的開采沉陷影響函數(shù)法優(yōu)化[J].煤炭學(xué)報,2016,41(1):271-276.doi:10.13225/j.cnki.jccs.2014.1736

Cai Yinfei，Verdel Thierry，Olivier Deck，et al.Improving the influence function method to take topography into the calculation of mining subsidence[J].Journal of China Coal Society,2016,41(1):271-276.doi:10.13225/j.cnki.jccs.2014.1736