0 引 言
煤矿资源开采后，开采区域周围的岩层和地表发生移动和变形，导致地表产生沉降，建构筑物遭到损毁，严重时还会引发地质灾害，因此有必要开展矿区地表沉降的监测和预计等工作[1]。矿区地表形变监测通常沿工作面走向或倾向布置剖面线状观测线，使用传统测量手段如水准仪和全站仪等定期或不定期观测，虽然该方法获取的数据精度较高，但工作量大、覆盖范围有限、且不易实现连续监测。合成孔径雷达差分干涉测量(Differential Interferometric Synthetic Aperture Radar,D-InSAR)是地表变形监测的新技术，可用于提取地表微小形变信息，因其全天时、全天候、监测范围广、精度高等优势，为矿区开采沉陷监测提供了新的手段。
20世纪90年代开始，国内外学者对InSAR在矿区的应用展开了大量研究，CARNEC等[2]最早利用D-InSAR技术监测煤炭开采导致的地表形变，得到了法国Gardanne地区的地表沉降分布；COLESANTI等[3]对57景ERS-1/2数据采用永久散射体InSAR(Persistent Scattered InSAR,PS-InSAR)技术监测法国Lorraine铁矿的地表稳定性，并对比水准数据验证了监测结果的可靠性；范洪冬等[4]对12景ERS1/ERS2数据采用小基线集InSAR(Small baseline subset InSAR，SBAS-InSAR)技术，获取了江苏省某矿区3年间的地表平均沉降速率分布，并发现在西部开采工作面附近沉降量最大；李达等[5]利用13景TerraSAR-X数据，通过SBAS-InSAR监测了陕西省某矿区工作面开采引起的地表沉降，结果表明该方法在矿区地表沉降监测方面的良好应用前景。
概率积分法(Probability Integral Method，PIM)是目前开采沉陷预计中应用最广泛的方法之一，预计时需确定多个参数的取值，其中与沉降有关的参数包括下沉系数q、主要影响角正切tan β、拐点偏移距Si(i=1,2,3,4)和开采影响传播角θ0。为得到上述参数，通常以线状观测站的实测数据为基础，利用线性近似法、模矢法和智能优化算法来确定最佳参数[6-7]。线性近似法和模矢法存在易发散、对初值要求高等不足，而智能优化算法模拟了生物觅食等自然现象或种群间合作与竞争过程，通过随机搜索达到全局最优的目的。近年来，学者们引入了许多智能优化算法来求取预计参数，如王正帅等[8]将粒子群算法(Particle Swarm Optimization，PSO)和文化算法(Cultural Algorithm,CA)结合，提出概率积分法文化-随机粒子优化算法，成功获取了山东省某矿区4326工作面的预计参数；王光磊[9]通过遗传算法(GeneticAlgorithm,GA)的运用反算出概率积分法的预计参数，并以此为依据建立了开采沉陷的预计模型；陈涛等[10]利用果蝇优化算法(Fruit Fly Optimization Algorithm,FOA)构造了拟合值与实测值均方差最小的适应度函数模型，获取了安徽省某矿区1013工作面的预计参数，同时比较遗传算法和粒子群算法的结果，得出该算法的精度最高。
以上工作在算法上的改进使得参数精度得到很大程度提升，但是仅依靠线状观测站的实测数据难以体现开采导致的全盆地沉降性质，如能获取开采影响范围内的全盆地沉降，以此为基础获取最佳预计参数，精度可得到进一步提升。鉴于此，笔者融合SBAS-InSAR技术和FOA，以我国内蒙古某矿区为例，选用27景Sentinel-1A数据，利用SBAS-InSAR技术获取2017年9月至2018年7月的矿区时序沉降，将其与实测水准数据对比验证后，提取221上06A工作面的全盆地沉降值，融合改进步长的FOA确定该地质采矿条件下的最佳预计参数。
1 试验方法与计算原理
1.1 SBAS-InSAR基本原理
时序InSAR技术基于D-InSAR技术发展而来，根据主影像的数量分为单主影像和多主影像2大类，其中多主影像以SBAS-InSAR技术为代表。SBAS-InSAR技术采用奇异值分解的方法，联合多个小基线集进行地表时序形变的反演[11-13]。假设获取某一地区的n+1景SAR影像，将影像按获取时间排列，选定其中某n幅景为公共主影像，完成剩余影像的配准与重采样，再设定时空干涉基线阈值，生成M对干涉图，其中M满足下式，即
(1)
设干涉图j由tA、tB(tA>tB)时间获取的影像生成，干涉图j中任意1个像素的干涉相位为
δφj=φ(tA)-φ(tB)=φdef,j+φtopo,j+φflat,j+φorbit,j+φatm,j+φnoise,j
(2)
式中：φ(tA)和φ(tB)分别为tA和tB时刻SAR影像的相位，φdef,j为卫星视线向的形变相位，φtopo,j、φflat,j、φorbit,j、φatm,j和φnoise,j分别为地形、平地效应、轨道误差、大气延迟和噪声导致的相位。
地形相位可通过引入外部DEM消除得差分干涉图；接着利用差分干涉图的幅度信息筛选失相干表现缓慢的候选像元点；再进行相位解缠，并通过理论公式和滤波等方式去除或减弱由地形、轨道和大气延迟等误差导致的相位，最后得到形变相位。通常设相邻时间间隔内的速率为线性，则形变相位可表示为
(3)
其中，矩阵A为系数阵，由于时空基线的设定使得矩阵A秩亏，故利用奇异值分解法联合求解上式的最小范数最小二乘解，最终得到观测时间内的平均沉降速率和时间序列沉降值。
1.2 融合FOA的概率积分参数优化原理
果蝇优化算法(Fruit Fly Optimization Algorithm，FOA)是基于果蝇群体在觅食过程中的合作与竞争行为的智能优化算法，果蝇群体利用嗅觉感知食物气味，根据气味判断食物大致位置，再利用视觉向气味浓度最大的地方飞行，从而准确判断食物位置[14-15]。FOA的大致流程为：设置初始参数，包括果蝇种群位置、种群规模、适应度函数及迭代次数；定义搜索步长，果蝇个体按步长更新位置；计算个体的气味浓度值即适应度函数值，标记浓度最高的个体为最佳个体；记录最佳个体的位置及浓度值，群体按视觉飞往该位置；按步长重复更新种群位置，迭代得到最终结果。融合FOA求取概率积分预计参数的原理及步骤如下：
设概率积分预计参数向量p=[q,tan β,S1,S2,S3,S4,θ0]，P为p的搜索空间，由概率积分原理得矩形工作面地表点(x,y)的下沉量W计算公式为
(4)
且
式中：w0为地表最大下沉量，当地表为充分采动时w0=qmcos α，(m)为煤层厚度，q为下沉系数，α为煤层倾角，当地表为非充分采动时w0需乘以地表采动程度系数[16]；D1和D3分别为工作面倾向和走向长；S1和S2分别为倾向下山和上山拐点偏移距；S3和S4分别为走向左和右边界拐点偏移距；l和L分别为工作面走向和倾向计算长度[17]。
以此预计参数按式(4)进行开采沉陷预计，并按该预计值与SBAS-InSAR所得沉降值差值的平方和最小为原则，将该问题转化为约束优化问题，在空间P中找到向量P，使P满足下式[10]，即
(5)
式中：f为适应度函数值；W为InSAR所得沉降值；w为预计沉降值；N为选取的特征点数量。
基于SBAS-InSAR和FOA的全盆地开采沉陷监测与参数计算的流程如图1所示。
图1 基于SBAS-InSAR和FOA的开采沉陷监测与参数计算流程
Fig.1 Flowchart of subsidence monitoring and parameters calculating by SBAS-InSAR and FOA
2 试验区域和数据简介
试验区域位于内蒙古自治区某矿区，井田内岩层包含三叠系上统、侏罗系中下统(含煤地层)、侏罗系中统、白垩系下统及第四系。其中第四系表土层平均厚度14 m，白垩系下统地层平均厚度338 m，岩性主要为中生代厚层砂岩，存在强度低，胶结性差，遇水泥化等特点，为弱胶结岩层覆岩条件下开采[18]，邻近矿区开采资料显示该覆岩条件下开采导致的地表沉降远小于我国中东部。矿区221上06A工作面平均采深690 m，采高5.67 m，采宽285 m，煤层倾角为2°，2017年9月开始回采，到2018年4月停采，走向累计推进890 m。为监测地表移动变形，在工作面上方布设了走向和倾向2条地表移动观测线，共计测点106个，其中走向自南向北，共设置50个，命名为F1～F50，倾向自西向东，共设置59个，命名为G1～G59，2017年9月15日首次观测，至2018年7月26日共观测14次。
试验数据选用Sentinel-1A卫星从2017年9月8日至2018年7月29日获取的共27景SAR影像，该卫星由欧洲太空局(ESA)于2014年发射，载有C波段合成孔径雷达，影像覆盖范围达20 km×25 km，方位向、距离向分辨率分别为5和20 m。笔者将原始影像裁剪后作为研究区域，同时选用美国宇航局(NASA)航天飞机测图任务(SRTM)90 m分辨率的DEM数据去除地形相位。
3 试验结果与分析
3.1 矿区时序沉降提取与精度分析
为处理27景Sentinel-1A数据，先以2018年2月11日的影像为公共主影像完成所有影像的配准与重采样；再设置短基线时空阈值，经干涉计算生成干涉图；然后去除地形相位、滤波和相位解缠；最后设定适当阈值选取高相干点，利用奇异值分解方法求解范围内的形变速率，并利用时间基线计算得点的时序沉降值[19]。结果表明，2017年9月至2018年7月间，由于221上06A工作面的开采，导致的地表最大沉降速率为154 mm/a，最大沉降量为201 mm，在工作面上方形成了地表移动盆地。绘制了视线向矿区地表沉降速率图和矿区时序沉降图，分别如图2和图3所示。根据矿区资料显示，221上06A工作面非该矿首采面，2016年9月至2017年12月间，存在221上17(首采面)和221上18两
图2 视线向矿区地表沉降速率
Fig.2 Subsidence rate of line of sight to surface
图3 矿区时序沉降
Fig.3 Time series subsidence of mining area
个工作面的开采，老采空区的影响导致地表沉陷并未与工作面的形状一致，移动盆地偏向老采空区方向，且本区表土层较薄，覆岩为中生代厚砂岩，使得地表沉降较小。
将221上06A工作面地表移动观测站的实测水准数据用于分析SBAS-InSAR所得沉降值的精度，选取未遭破坏且附近存在高相干性点的水准点(未包含最大沉降点)共计40个，分别绘制走向和倾向方向走的沉降对比折线图，如图4所示，对比计算的结果见表1，SBAS-InSAR平均沉降量误差为3.48 mm，标准差为51.16 mm，均方根误差为10.64 mm。结果表明，在该工作面的应用上SBAS-InSAR的精度较高，效果较好。
表1 观测线上水准与SBAS-InSAR沉降对比
Table 1 Comparison of settlement between leveling and SBAS-InSAR on observation lines
点号W/mm水准SBAS-InSAR绝对误差/mm相对误差/%F11219758F316271169F924351146F2280941417F2388981011F249610266F2510310522F261051161110F3013813800F31158145-138F32162150-127F33173164-95F34180170-106F35184170-148F36188173-158F46154148-64F47140133-75F48133129-4-3F4911111433F5010110544G431421135G148178-33G17949622G1887981113G221151281311G24128121-75G2514014000G26148143-53G29171151-2012G37171162-95G40153140-13-8G41140130-107G42127115-129G4410297-55G505440-1426G514930-19-39G523922-1744G542713-1452G552210-1255G57149-5-36
图4 观测线上水准与SBAS-InSAR沉降对比
Fig.4 Comparison of settlement between leveling and SBAS-InSAR of observation lines
3.2 融合FOA的概率积分参数确定步骤及结果
利用第3.1节的结果融合FOA求取概率积分预计参数。先设置FOA相关参数，初始概率积分预计参数设置为q=0.6，tan β=1.6，s=0.05 H(H为采深)，θ0=89°；最大迭代次数设为N=100；将预计参数按p=[q,tan β,S1,S2,S3,S4,θ0]组成向量，再组成预计参数矩阵，其中矩阵行数为种群规模i=50，列数为参数个数j=7；考虑近水平煤层开采，为简化处理，未改变开采影响传播角θ0，其余种群初始搜索步长设置为b0=[0.01,0.01,1,1,1,1]；再按沉降值10 mm确定221上06A工作面地表移动盆地边界点，并剔除受周围老采空区影响的点，共提取位于盆地内的特征点3 723个，如图5所示；由于该工作面为非充分采动，预计参数向量按式(4)计算非充分采动条件下的地表沉降值；接着利用SBAS-InSAR
图5 221上06A工作面地表移动盆地
Fig.5 Surface moving basin of 221 upper 06A working face
得到的沉降值和预计参数得到的拟合值按式(5)计算浓度函数值，选出最优浓度函数值的预计参数向量；为提高参数求取的效率和精度，改进搜索步长并更新预计参数矩阵，设第k(2≤k≤N)次迭代的搜索步长为bk，最优浓度函数值为fk，若kfk<fk-1，则减小步长，
反之若fk≥fk-1，则增大步长，
最后按最大迭代次数N重复以上步骤，调试分析得最佳预计参数向量pbest，其中q=0.33，tan β=1.83，s=0.045H，θ0=89°，由于该工作面为非充分采动，下沉系数q已根据文献[16]修正。参数求取过程中所用点数达3 723个，随机展示了其中30个拟合值和SBAS-InSAR监测值的对比计算结果见表2。此外，由于地表移动变形尚未稳定，所求结果应为动态沉陷过程的预计参数，因此选取绘制了SBAS-InSAR监测得到的最大下沉点时序沉降情况如图6所示。结果表明，2018年5月工作面停采后，地表最大下沉点的残余沉降并不明显，其中5月底至7月底的沉降量仅为6 mm，故可认为所求动态参数虽略小于终态预计参数，但符合该地质采矿条件，适用于该矿区待采工作面的预计工作，后续可进一步根据地表稳定后的最终移动变形，完善终态预计参数的求取。
表2 部分拟合结果与SBAS-InSAR结果对比计算
Table 2 Comparison of some SBAS-InSAR results and fitted results
点号W(SBAS)/mm拟合值/mm绝对误差/mm相对误差/%1181.67214.27-32.60182173.93205.03-31.10183163.31170.68-7.3654151.92141.2110.7175147.32138.289.0466143.09146.49-3.3927130.86124.686.1758128.67139.33-10.6689119.80100.0219.771710113.49129.35-15.861411101.7390.0311.69111295.8580.1715.67161390.9085.025.8861477.3371.156.1781568.2769.38-1.11161656.4562.04-5.59101748.9950.99-2.0041846.1251.00-4.88111938.4945.19-6.70172036.0339.43-3.3992132.5636.68-4.12132230.8928.392.5082326.9633.52-6.57242420.5119.011.5072518.1316.671.4782614.8214.350.4732711.3311.51-0.182288.026.861.1614296.815.311.5022305.014.390.6312
图6 221上06A工作面最大下沉点时序沉降
Fig.6 Time series settlement of maximum subsidence point of 221 upper 06A working face
4 结 论
1)利用27景Sentinel-1A数据和SBAS-InSAR技术监测2017年9月至2018年7月间内蒙古某矿区地表变形，获得221上06A工作面的全盆地沉降情况，其最大沉降量为201 mm。研究表明SBAS-InSAR技术能够较好地监测矿区工作面开采引起的地表全盆地沉陷。
2)针对全盆地开采沉陷监测结果，采用改进步长的果蝇优化算法，实现了全盆地3 723个点共同求取概率积分预计参数，全盆地求参效果良好。
3)采用改进步长的果蝇优化算法进行全盆地求参获得的本区概率积分预计参数为：q=0.33，tan β=1.83，s=0.045 H，θ0=89°，为后续开采沉陷预计提供了科学依据。
参考文献(References):
[1] 何国清,杨 伦,凌庚娣,等.矿山开采沉陷学[M].徐州:中国矿业大学出社,1991:1-2.
[2] CARNEC C,MASSONNET D,KING C.Two examples of the use of SAR interferometry on displacement fields of small spatial extent[J].Geophysical Research Letters,1996,23(24):3579-3582.
[3] COLESANTI C,MOUELI S L,BENNANI M,et al.Detection of mining related ground instabilities using the permanent scattered technique :a case study in the east of France[J].International Journal of Remote Sensing,2012,26(1):201.
[4] 范洪冬,邓喀中,薛继群,等.利用时序SAR影像集监测开采沉陷的试验研究[J].煤矿安全,2011,42(2):15-18.
FAN Hongdong,DENG Kazhong,XUE Jiqun,et al.An experimental research on using time series SAR images to monitor mining subsidence[J].Safety in Coal Mines,2011,42(2):15-18.
[5] 李 达,邓喀中,高晓雄,等.基于SBAS-InSAR的矿区地表沉降监测与分析[J].武汉大学学报(信息科学版),2018,43(10):1531-1537.
LI Da,DENG Kazhong,GAO Xiaoxiong,et al..Monitoring and analysis of surface subsidence in mining area based on SBAS-InSAR[J].Geomatics and Information Science of Wuhan University,2018,43(10):1531-1537.
[6] 朱晓峻,郭广礼,方 齐.概率积分法预计参数反演方法研究进展[J].金属矿山,2015(4):173-177.
ZHU Xiaojun,GUO Guangli,FANG Qi.Recent progress on parameter inversion of probability integral method[J].Metal Mine,2015(4):173-177.
[7] 冉 典,吴 捷.模矢法在厚松散层地表移动参数解算中的应用[J].地矿测绘,2014,30(2):25-27,35.
RAN Dian,WU Jie.Application of pattern search method in determining surface movement parameters of thick loose bed[J].Surveying and Mapping of Geology and Mineral Resources,,2014,30(2):25-27,35.
[8] 王正帅,邓喀中,康建荣.概率积分法参数反演的文化-随机粒子群优化算法[J].辽宁工程技术大学学报(自然科学版),2013,32(3):311-315.
WANG Zhengshuai,DENG Kazhong,KANG Jianrong.Random PSO embedded cultural framework for parameters inversion of probability integral method.Journal of Liaoning Technical University (Natural Science),2013,32(3):311-315.
[9] 王光磊.基于遗传算法的矿区开采沉陷预计参数研究[D].阜新：辽宁工程技术大学,2014.
[10] 陈 涛,郭广礼,朱晓峻等.利用果蝇算法反演概率积分法开采沉陷预计参数[J].金属矿山,2016(6):185-188.
CHEN Tao,GUO Guangli,ZHU Xiaojun et al.Mining subsidence prediction inversion of the probability integral method based on fruit flies algorithm[J].Metal Mine,2016(6):185-188.
[11] 朱建军,李志伟,胡 俊.InSAR变形监测方法与研究进展[J]．测绘学报,2017,46(10):1717-1733.
ZHU Jianjun,LI Zhiwei,HU Jun.Research progress and methods of InSAR for deformation monitoring[J].Acta Geodaetica et Cartographica Sinica,2017,46(10):1717-1733.
[12] LANARI R,Mora O,MANUNTA M,et al.A Small-Baseline approach for investigating deformations on full-resolution differential SAR interferograms[J].IEEE Transactions on Geoscience and Remote Sening,2004,42(7):1377-1386.
[13] ZHANG L,ZHONG L U,DING X,et al.Mapping ground surface deformation using temporarily coherent point SAR interferometry:application to Los Angeles Basin[J].Remote Sensing of Environment,2012,117(1):429-439.
[14] PAN W.A new fruit fly optimization algorithm:taking the financial distress model as an example[J].Knowledge-based System,2011,26(7):67-74.
[15] 霍慧慧.果蝇优化算法及其应用研究[D].太原:太原理工大学,2015.
[16] 胡炳南,张华兴,申宝宏，等.建筑物、水体、铁路及主要井巷煤柱留设与压煤开采指南[M].北京:煤炭工业出版社,2017:1-2
[17] 邓喀中,谭志祥,姜 岩，等.变形监测与沉陷工程学[M].徐州:中国矿业大学出版社,2014:217-220.
[18] 孙利辉,纪洪广,杨本生.西部典型矿区弱胶结地层岩石的物理力学性能特征[J].煤炭学报,2019,44( 3) :865-873.
SUN Lihui,JI Hongguang,YANG Bensheng.Physical and mechanical characteristic of rocks with weakly cemented strata in Western representative mining area[J].Journal of China Coal Society,2019,44( 3) :865-873.
[19] BERARDINO P,FORNARO G,LANARI R,et al.A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms[J].IEEE Transactions on Geoscience and Remote Sensing,2003,40(11):2375-2383.
[20] 司刚全,李水旺,石建全,等.采用改进果蝇优化算法的最小二乘支持向量机参数优化方法[J].西安交通大学学报,2017,51(6):14-19.
SI Gangquan,LI Shuiwang,SHI Jianquan et al.Least squares support vector machine parameters optimization based on improved fruit fly optimization algorithm with application[J].Journal of Xi’an Jiaotong University,2017,51(6):14-19.
[21] 朱志同,郭 星,李 炜.新型果蝇优化算法的研究[J].计算机工程与应用,2017,53(6):40-45.
ZHU Zhitong,GUO Xing,LI Wei.Research on new fruit fly optimization algorithm[J].Computer Engineering and Applications,2017,53(6):40-45.
