本发明涉及采矿安全领域，尤其涉及一种基于地面沉降监测数据的两层采空区地面沉降的预测方法。
背景技术
在过去的十年里，中国煤炭的年产量占全球年产量的40％以上。由于煤炭开采所产生的地面沉降不仅造成财产和生命的损失，而且还造成一系列地质环境问题。因此，对采空区地面沉降进行预测具有经济效益、社会效益、环境效益。现有的单层采空区地面沉降预测方法主要包括理论分析方法、物理模拟方法、数值模拟方法。
为了提高煤炭开采量以及经济发展的需要，煤炭的多层开采是很有必要的。然而，由于多层开采的剖面、形状和大小与单层开采不一样，单层开采时地面沉降的预测方法不适用于多层开采。
目前，一些学者对多层采空区地面沉降进行了研究，他们采用的方法主要是理论计算(主要是影响函数法)、物理模拟方法、数值模拟方法。影响函数法是一种理论方法，通用性较好，但是该方法受到地质条件、开采条件以及开采边界等的影响，其适应性受到限制。而物理模拟则费时费钱，其适应性也受到影响。
相比而言，数值模拟方法在适应性和经济性都具有优势，然而目前有至少2条不利因素制约了数值模拟的应用。第一，数值模拟的计算结果受岩体力学参数的影响。一般计算时参数选取工程勘察阶段测试的岩体力学参数，然而这些参数都是通过实验室内对岩块进行测试得到的参数，不是现场岩体的实际参数，所以数值模拟的计算结果往往与实际监测的地面沉降值相差较大。第二，在开采完上层煤后，部分上覆岩体会屈服，则这部分屈服的岩体力学参数要折减。如果未对已屈服岩体的力学参数进行折减，则计算结果偏小。
因此，目前对多层采空区地面沉降进行的研究尚有诸多不足，难以相对准确和客观地预测沉降值，有待对研究方法做进一步的改进。
技术实现要素：
本发明的目的在于解决以上的问题，提供一种可靠的两层采空区地面沉降的预测方法，包括以下步骤：
步骤100，在上层煤层即将开采前，在待开采范围内划定一定面积的区域，在该区域内设置地面沉降观测站和观测线，对观测站的高程进行长期观测，得到相应观测点因上层煤层开采导致的地面沉降值；
步骤200，对煤层上覆岩土体钻孔取样，通过室内实验确定上覆岩土体每一地层的物理和力学参数；
步骤300，根据钻孔资料和地形数据建立整个煤矿的三维地质模型，并对三维地质模型划分网格；
步骤400，选择煤层上覆岩土体的屈服准则，确定边界条件及初始应力；
步骤500，对步骤200得到的力学参数进行折减，结合步骤200得到的物理参数，定为初始参数；并且，所述步骤500相对于前述的步骤300，以及相对于前述的步骤400，无先后顺序限制；
步骤600，由步骤500得到的初始参数，利用步骤400确定的屈服准则、边界条件、初始应力，利用数值模拟技术，基于步骤100观测得到的地面沉降值，通过反演试算的方法确定上覆岩土体最终的力学参数；
步骤700，根据步骤600得到的力学参数，利用数值模拟计算整个煤矿开采上层煤导致的地面沉降值s1；
步骤800，开采上层煤后，对两层煤之间地层赋物理力学参数，对已经屈服的单元的力学参数进行折减，利用数值模拟计算和预测整个煤矿开采下层煤导致的地面沉降值s2，则s2和步骤700得到的s1之和即为开采两层煤的地面沉降最终值s。
优选地，步骤100中，地面沉降观测手段选择gps、北斗导航、全站仪、水准仪中的一种或多种，根据采深采厚比初步确定的地面沉降强度进行选择。
优选地，步骤100中，地面沉降观测时间为从上层煤开采前直到开采结束并且地面沉降稳定；观测频率5-10天观测一次；地面沉降稳定的判断标准是单个变形观测点6个月的沉降量总和小于30mm。
优选地，步骤200中，物理参数主要包括密度，力学参数主要包括单轴抗压强度、抗拉强度、弹性模量、泊松比、内摩擦角和粘聚力。
优选地，步骤300中，三维地质模型的边界根据煤层开采范围确定，网格划分的数量由计算机的计算能力确定。
优选地，步骤400中，屈服准则选择摩尔-库伦准则。
优选地，步骤500中，所述折减具体为弹性模量折减50％，泊松比增加15％，作为煤层上覆岩土体力学参数初始值。
优选地，步骤600中，所述反演试算可以通过正交试验来完成，具体为：
步骤601，步骤400屈服准则选择摩尔-库伦准则，试验因素选择弹性模量、泊松比、内摩擦角、粘聚力4个参数；
步骤602，定义步骤100观测得到的地面沉降平均值与数值模拟得到的地面沉降平均值的比值为地面沉降平均比，并选为试验因子；
步骤603，根据正交试验设计理论，生成多因素多水平正交试验表，多因素是指弹性模量、泊松比、内摩擦角、粘聚力4个试验因素，多水平是指每一个试验因素均对应多种梯度的数值；
步骤604，把正交试验表中的试验方案确定的多组参数分别给三维地质模型赋值，并进行数值模拟，计算得到每一组试验方案的地面沉降平均比；
步骤605，比较分析地面沉降平均比的数值大小，该数值最接近于1时的力学参数与实际岩土体的力学参数最接近，则该组力学参数取为上覆岩土体最终的力学参数。
优选地，步骤800中，对两层煤之间地层所赋的物理力学参数与上层煤顶板地层的物理力学参数相同。
优选地，步骤800中，所述折减具体为弹性模量折减99.5％。
本发明与现有技术相比，其显著特点是：
(1)本发明数值计算所需要的上覆岩土体力学参数是基于现场地面沉降监测数据通过反演方法得出的，参数的可靠性高。所以，通过数值计算得到的地面沉降值可靠性高。
(2)本发明并非建立地面沉降监测区域的三维地质模型，而是建立整个煤矿的三维地质模型，并在此基础上依据观测点的地面沉降数据进行力学参数反演，故所反演的力学参数可靠性高。
(3)将实验室内测得的岩块的力学参数进行折减后作为反演的初始力学参数，这大大提高了计算效率、节省计算时间，也防止反演得出的最终力学参数不是煤层上覆岩土体的实际参数。
(4)本发明提出用正交试验的方法来反演煤层上覆岩土体力学参数，该方法逻辑简单、可操作性强。
(5)开采上层煤后，将发生屈服的上覆岩土体单元的参数进行折减，然后模拟和计算开采下层煤的沉降数值，这与事实相符，所以计算得到的最终沉降值可靠性高。
附图说明
图1是根据本发明一个实施例的预测两层采空区地面沉降的方法的流程示意图；
图2为地层剖面图；
图3上层煤(2-2煤)工作面布置图，图中的abcd即位于222201工作面的gps地面沉降观测区域，m1～m6为数值模拟的地面沉降观测线；
图4下层煤(4-1煤)工作面布置图，m1～m6为数值模拟的地面沉降观测线；
图5为222201工作面的gps地面沉降观测区域的基准点、观测点、观测线布置图；
图6为划分网格后的三维地质模型图；
图7为gps观测得的地面沉降曲线与flac3d数值模拟得到的地面沉降曲线对比图；
图8(a)为岩石材料屈服后的应变软化现象示意图，(b)为破坏前、后弹性模量的相对变化关系示意图；
图9-图14为数值模拟的地面沉降观测线m1-m6中的开采两层煤后最终地面沉降预测值s和开采上层煤时地面沉降预测值s1对比示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
在本发明的描述中，需要理解的是，术语“包括/包含”、“由……组成”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的产品、设备、过程或方法不仅包括那些要素，而且需要时还可以包括没有明确列出的其他要素，或者是还包括为这种产品、设备、过程或方法所固有的要素。在没有更多限制的情况下，由语句“包括/包含……”、“由……组成”限定的要素，并不排除在包括所述要素的产品、设备、过程或方法中还存在另外的相同要素。
如图1所示，本发明的一个实施方式的两层采空区地面沉降的预测方法，包括如下步骤：
步骤100，在上层煤层即将开采前，在待开采范围内划定一定面积的区域，在该区域内设置地面沉降观测站和观测线，对观测站的高程进行长期观测，得到相应观测点因煤层开采导致的地面沉降值。
本步骤中，观测站可设置硬化混凝土墩，观测点主要包括基准点和变形观测点，其中基准点为远离变形区域且稳定可靠的已知点，变形观测点位于变形区域内观测站上，观测点的连线构成地面沉降观测线，如图4所示，地面沉降观测线一般有两条，一条沿煤层走向方向，另一条沿煤层倾向方向；相邻两个观测站的距离由开采深度确定，一般5～25米；观测线的长度要大于地面沉降的范围。
本步骤中，地面沉降观测方式可采用gps、北斗导航、全站仪或静力水准仪，具体采用哪种方式，以及采用其中的一种或几种方式组合，根据采深采厚比初步确定的地面沉降强度进行选择。
本步骤中，地面沉降观测时间为从上层煤开采前直到开采结束并且地面沉降稳定；观测频率一般5-10天观测一次；地面沉降稳定的判断标准是单个变形观测点6个月的沉降量总和小于30mm。
步骤200，对煤层上覆岩土体钻孔取样，通过室内实验确定上覆岩土体每一地层的物理和力学参数。
本步骤中的物理参数包括密度，力学参数包括单轴抗压强度、抗拉强度、弹性模量、泊松比、内摩擦角和粘聚力。
步骤300，根据钻孔资料和地形数据建立整个煤矿的三维地质模型，并对三维地质模型划分网格。
本步骤中的三维地质模型的边界根据煤层开采范围确定，网格划分的数量由计算机的计算能力确定，网格划分类型可以选择四面体和六面体混合网格。
步骤400，选择煤层上覆岩土体的屈服准则，确定边界条件及初始应力。
本步骤中的屈服准则可以选择摩尔-库伦准则，也可以选择岩土力学中适用的其他准则。
步骤500，对步骤200得到的力学参数进行折减，结合步骤200得到的物理参数，定为初始参数。
本步骤中，综合考虑步骤200实验得到的力学参数，并参考相关的研究，对步骤200得到的岩石的力学参数中的弹性模量折减50％，泊松比增加15％，作为煤层上覆岩土体力学参数初始值，这种折减方案更接近岩石的真实力学参数。
需要说明的是，步骤500相对于步骤200和步骤300并没有严格的先后顺序限制，步骤500也可置于步骤200和步骤300之间进行，即步骤200，对煤层上覆岩土体钻孔取样，通过室内实验确定上覆岩土体每一地层的物理和力学参数；步骤300，对步骤200得到的力学参数进行折减，结合步骤200得到的物理参数，定为初始参数；步骤400，根据钻孔资料和地形数据建立整个煤矿的三维地质模型，并对三维地质模型划分网格；步骤500，选择煤层上覆岩土体的屈服准则，确定边界条件及初始应力；步骤600，由步骤300得到的初始参数，利用步骤500确定的屈服准则、边界条件、初始应力，利用数值模拟技术，基于步骤100观测得到的地面沉降值，通过反演试算的方法确定上覆岩土体最终的力学参数。
步骤500还可置于步骤300和步骤400之间进行，步骤200，对煤层上覆岩土体钻孔取样，通过室内实验确定上覆岩土体每一地层的物理和力学参数；步骤300，根据钻孔资料和地形数据建立整个煤矿的三维地质模型，并对三维地质模型划分网格；步骤400，对步骤200得到的力学参数进行折减，结合步骤200得到的物理参数，定为初始参数；步骤500，选择煤层上覆岩土体的屈服准则，确定边界条件及初始应力；步骤600，由步骤400得到的初始参数，利用步骤500确定的屈服准则、边界条件、初始应力，利用数值模拟技术，基于步骤100观测得到的地面沉降值，通过反演试算的方法确定上覆岩土体最终的力学参数。
步骤600，由步骤500得到的初始参数，利用步骤400确定的屈服准则、边界条件、初始应力，利用数值模拟技术，基于步骤100观测得到的地面沉降值，通过反演试算的方法确定上覆岩土体最终的力学参数；
本步骤中，数值模拟技术可以通过flac3d软件来实现。
本步骤中，反演试算可以通过正交试验来完成，具体为：
步骤601，步骤400屈服准则选择摩尔-库伦准则，试验因素选择弹性模量、泊松比、内摩擦角、粘聚力4个参数；
步骤602，定义步骤100观测得到的地面沉降平均值与数值模拟得到的地面沉降平均值的比值为地面沉降平均比，并选为试验因子；
步骤603，根据正交试验设计理论，生成多因素多水平正交试验表，多因素是指弹性模量、泊松比、内摩擦角、粘聚力4个试验因素，多水平是指每一个试验因素均对应多种梯度的数值；
步骤604，把正交试验表中的试验方案确定的多组参数分别给三维地质模型赋值，并进行数值模拟，计算得到每一组试验方案的地面沉降平均比；
步骤605，比较分析地面沉降平均比的数值大小，该数值最接近于1时的力学参数与实际岩土体的力学参数最接近，则该组力学参数取为上覆岩土体最终的力学参数。
步骤700，根据步骤600得到的上覆岩土体最终的力学参数，利用数值模拟计算整个煤矿开采上层煤导致的地面沉降值s1。
步骤800，开采上层煤后，对两层煤之间的地层赋物理力学参数，对已经屈服的单元的力学参数进行折减，利用数值模拟计算和预测整个煤矿开采下层煤导致的地面沉降值s2，则s2和步骤700得到的s1之和即为开采两层煤的地面沉降最终值。
本步骤中，两层煤之间地层岩性与上层煤顶板地层岩性相同，所以两层煤之间地层所赋的物理力学参数与上层煤顶板地层的物理力学参数相同。
本步骤中，关于对已经屈服的单元的力学参数进行折减，由于岩石材料在屈服后一般呈现应变软化现象(图3a)，即受应力作用超过屈服点后，其应力-应变曲线斜率(或仅仅是应力本身)减小的现象。对于洞室而言，屈服的岩体可以自由地从顶板跌落，所以一般假定应力超过屈服点后便降为零或接近零(图3a)。将屈服准则选择摩尔-库伦准则，屈服状态后，岩体破坏，其参数发生改变，所以需要按照破坏后的状态对岩体参数进行更新(即折减)。
根据胡克准则，弹性的岩石材料有：e＝σ/ε，其中e为弹性模量，σ为应力，ε为应变。同样，根据胡克准则，破坏后的岩石材料可以认定为弹性模量减小的弹性材料。具体来说，对于应力σ达到极限应力σu时，岩石材料的弹性模量减小到一个很小的值，即ed<<ee，其中ee和ed分别为岩石材料破坏前、后的弹性模量。于是，把屈服后发生应变软化阶段的岩石材料看作为弹性模量很小的弹性材料，即把图3a所示的洞室破坏假定为图3b所示的斜率为很小值ed。根据图3b，可知破坏前、后应力的比值同弹性模量的变化具有对应关系，即σd/σu＝ed/ee＝d，其中σd和σu分别为刚刚破坏时的应力和极限应力σu。只需把比值d定为一个很小的值(比如0.5％)，即可保证应力超过屈服点后便降为零或接近零，并且屈服后的应力值不随着应变增加而增加。
基于上述分析，本发明认为对多层煤层开采进行数值模拟时，当上覆岩土体发生屈服时，需要对相应的弹性模量进行折减99.5％，即屈服后的弹性模量仅为初始弹性模量的0.5％，其余的参数则不变，对于多层煤层的开采，这种折减方案能最大程度地接近上覆岩土体发生屈服时岩石的真实力学参数。本步骤中，通过flac3d内置fish语言，将开采上层煤而发生屈服的单元找到，将这些单元的弹性模量替换成屈服后的弹性模量。相关代码如下：
由于在flac3d模拟过程中，当模型中一部分单元的弹性模量发生变化，则那部分单元将发生应力重分布，进而产生沉降，导致对下层煤开采进行数值模拟计算时，地面沉降值包括了由于单元应力重分布而产生的沉降。在本步骤中，在对下层煤开采进行模拟之前，对模型进行以下处理。
(1)将弹性模量折减之后的模型进行solve，让参数折减后的单元产生沉降位移；
(2)将上层煤开采所产生的地面沉降与由于参数折减后网格产生的沉降位移清零。
对模型进行上述处理后，下层煤开采之后，便可得到由下层煤开采导致的地面沉降s2，最终的地面沉降值即为上述s2与由上层煤开采导致的地面沉降s1之和。
以下给出本发明的一个具体应用的实施例。
某煤矿主采2-2煤与4-1煤。其中2-2煤煤层平均厚度6.0m，煤层平均埋深165m。4-1煤平均厚度5m，煤层平均埋深218m。该煤矿具体地层剖面图见图2。该煤矿于2008年11月开工建设，2011年5月矿井建成进行煤炭开采。采煤方法均采用走向长壁后退式综采的一次采全高开采方法，全部垮落法管理顶板。该煤矿的地形和工作面布置方式见图3和图4。
在2011年5月开采上层煤层之前，在待开采范围内划定一定面积的区域(即图3中222201工作面)，在该区域内设置地面沉降观测站和观测线，如图5所示。地面沉降观测线一般有两条，一条沿煤层走向方向，另一条沿煤层倾向方向；走向线一共布置110个变形观测点(r1-r110)，其中r1～r80观测点间距为20m，r80～r110观测点间距为10m。倾向线一共布置52个变形观测点(d1-d52)，观测点间距10m。
采用gps和全站仪对测点的高程进行测量。观测时间一般在一年以上，每隔7天对变形观测点进行一次高程测量。对沉降量较大的观测点，应当增加高程测量频率。当每一个变形观测点6个月的沉降量总和小于30mm，认为地面沉降基本完成。
通过对比煤层开采前、后观测点的高程数据，可得到由于煤层开采所产生的地面沉降值。
对煤层上覆岩土体钻孔取样，通过室内实验确定上覆岩土体每一地层的物理和力学参数。根据钻孔资料和地形数据建立整个煤矿的三维地质模型，并对三维地质模型划分网格。选择煤层上覆岩土体的屈服准则为摩尔-库伦准则，确定边界条件及初始应力。图6为划分网格后的三维地质模型。
对上述力学参数中的弹性模量折减50％，泊松比增加15％，作为煤层上覆岩土体力学参数初始值。
将物理参数和力学参数的初始参数，利用flac3d数值模拟技术，结合正交试验通过反演试算的方法确定上覆岩土体最终的力学参数，表1给出了通过反演试算得到的上层煤上覆岩土体最终的力学参数，图7给出了该煤矿地面沉降观测区域内的沉降曲线与flac3d数值模拟的沉降曲线对比。
表1通过反演试算得到的上层煤上覆岩土体最终的力学参数
将上覆岩土体最终的力学参数，利用数值模拟计算整个煤矿开采上层煤导致的地面沉降值s1。
开采上层煤后，对两层煤之间地层(即图2所示的第7层)赋于与第5层相同的物理力学参数，对已经屈服的单元的力学参数进行折减，利用数值模拟计算和预测整个煤矿开采下层煤导致的地面沉降值s2，则s2和s1之和即为开采两层煤的地面沉降最终值s。
图9～图14给出了数值模拟地面沉降6条观测线m1～m6(具体位置见图3和图4)中的开采两层煤后最终地面沉降预测值s(图中实线所示)和开采上层煤时各测线的地面沉降预测值s1(图中虚线所示)。
沉降预测结果表明，在开采下层煤后，图3和图4中的m3监测线所处的工作面的最大沉降可达11.6m，这比仅开采上层煤后的最大沉降增加了6.4m，较上层煤层开采后地面沉降值增加超过120％，如图11所示。这表明该处工作面顶板会发生整体垮落，这会殃及地面的设施和人员安全。为此，煤矿在开采下层煤之前对该处采空区做好了防范措施，并制定了相应的应急预案。这有力杜绝了发生物毁人亡的可能，极大地提高了煤矿生产的安全性。
通过上述预测，上层煤开采后，导致地面沉降，观测点很可能会被破坏，此时已无法继续使用，无法获得比较准确的沉降观测值，通过前期的观测并借助数值模拟，计算和预测整个煤矿开采下层煤将导致的地面沉降值，由此能够有效地获取开采下层煤将导致的地面总的沉降值，为后期在此区域内进行其他的生产生活作业提供可靠且必要的参考。
以上实例，仅用为本发明较佳的具体实施方式，但本发明的保护范围并不仅限于此，任何熟悉本技术领域的技术人员，在本发明披露的技术范围内，可想到的变化或等同替换，都应涵盖在本发明的保护范围之中。因此，本发明的保护范围应该以权利要求书的保护范围为准。