基于Google Earth Engine的采煤沉陷水体方向变化自动识别

赵艳玲，丁宝亮，何厅厅，肖 武，任 河

(1.中国矿业大学(北京) 地球科学与测绘工程学院，北京 100083；2.浙江大学 公共管理学院，浙江 杭州 310058)

能源是一个国家可持续发展的基础和重要支柱，对于我国而言，煤炭资源的地位尤为突出，占我国化石能源基础储量的94%左右。但随着煤炭的持续开采，势必会造成严重的环境问题，如地表沉陷、土壤污染等。据统计，平均每开采1万t煤，就有0.2～0.3 ha土地受到开采沉陷的影响。由于煤矿开采引起的地表扰动一直是一个全球性的问题，在美国、澳大利亚等国家由于地表下沉产生的沉陷水体对当地生态造成了难以修复的破坏。在我国东部高潜水位地区，如安徽，地下水极易由于地表沉降而显露出来，从而形成地面沉陷积水区，大面积的沉陷水体会损毁大量优质耕地，影响当地的生态系统健康。据安徽省淮南市政府统计，至2017年淮南市采煤塌陷区面积已达3万ha，近一半塌陷坑被水覆没，且由于地形平坦，土地复垦和回填的成本较高，大部分塌陷坑长期处于积水状态。因此，为了制定合理有效的复垦方案，对塌陷区积水的时空变化进行持续监测是有必要的。

遥感技术在监测地表自然资源变化方面已经非常成熟，基于遥感图像的土地利用变化监测已广泛用于采矿方向，如基于SAR技术监测矿区地表沉降，其通过跟踪 SAR图像像素偏移来监测地表的位移。对沉降边界进行预测的方法主要有时间函数模型、基于实测数据的预测模型、力学和数值模拟等方法。但是这些方法往往需要采煤相关数据，且处理速度较慢，自动化程度低。对于矿区沉陷积水的变化研究一般是基于3S技术，从水体的“大小”、“形状”和“变化率”角度出发，在一定的时间和空间尺度进行监测，大量研究采用水体面积年均变化幅度、年均变化速度以及水体向其他土地利用类型转移的列联矩阵等方法来描述沉陷积水的时空变化。随着GIS技术的进一步发展，出现了Digital Shoreline Analysis System(DSAS)方法，其通过绘制基线来对研究区边界进行监测。STEELE和HEFFERNAN运用海岸线发展系数Shoreline Development Factor(SDF)来描述水体的大小分布、连通性和形状；BHAGAWAT以城市中心为中点，分别向东、东南、南、西南、西、西北、北、东北共8个方向发射射线，通过这8条射线来监测该城市边界的空间位置变化，并取得了较高的精度。但在采煤沉陷水体的研究上，对沉陷水体各方向的空间位置变化进行监测的研究还较少。有研究发现，随着煤炭的持续开采，地面形成的沉陷区通常会随着开采工作面而呈现规律性的变化，而沉陷积水区也会出现类似的变化。通常，煤炭开采分为横向单煤层多工作面开采和纵向多煤层工作面开采。在横向单煤层多工作面开采过程中，随着工作面推进，地表沿该方向依次沉降，沉陷积水多呈现条带状；在纵向多煤层开采过程中，随着开采程度加深，地表沿垂直工作面方向向四周沉降，沉陷积水多呈现圆状。在煤炭开采中，2种开采方式往往同时存在。以上的方法与结论为本次研究提供了一条新思路。谷歌大数据云计算平台GEE能很快地调出大量遥感数据，并且其内部封装了大量函数可供使用，也逐渐被应用到采矿领域。在进行矿区沉陷水体的研究上，YI等构建了具有空间特征的土壤水分时间序列轨迹，并通过LandTrendr算法，对采煤沉陷区进行时空动态分析；YANG等采用改进归一化水体指数(mNDWI)和最大类间方差(OTSU)图像分割算法提取了1988—2018年中国东部高地下水位矿区的塌陷湿地。由此可见，GEE平台在进行矿区沉陷水体的研究上具有处理速度快、处理数据量大等的特点，大有发展前景。但平台上还尚未出现类似DSAS等的相关函数，需要使用者通过编写代码绘制基线(射线)从而实现监测功能。

笔者利用遥感大数据平台GEE，基于1989—2016年潘谢矿区沉陷水体数据，综合几何学和统计学知识，对每个独立的沉陷水体分别绘制射线，通过多组射线来分析潘谢矿区沉陷水体各方向的年际形状变化，并对其进行拟合，预测其变化情况。此方法可在缺少采煤相关信息的条件下仅依据遥感技术实现，且能达到较高的自动化程度，为矿区内沉陷水体的监测和生态修复方案的制定提供一定的依据。

1 研究区概况及数据源

1.1 研究区概况

潘谢矿区位于安徽省淮南市(图1)，地理坐标：116.33°E—116.90°E，32.72°N—32.93°N，南北跨度25 km，东西跨度60 km。矿区地形平坦，标高30～40 m，无承压地下水埋深为1.5 m左右，当矿区地面沉降达1.5 m时可能出现积水。气候属温带季风气候，年均气温15 ℃，年均降水量970 mm。矿区内共有谢桥矿、张集矿、顾北矿、顾桥矿、丁集矿、潘三矿、朱集矿、潘北矿、潘二矿和潘一矿等10个矿山企业。潘谢矿区的开采历史多达110余年，煤炭资源丰富，有9～18层可采煤层，平均厚度为20～30 m，属于近水平厚松散层煤层群。由于长期开采，地面经历多次沉陷，长期处于不稳定沉陷状态，地面累计沉陷深度最大可到20～30 m，加之当地较高的地下水位，形成了大面积的沉陷积水区，对当地房屋、耕地造成了严重影响。

图1 研究区位置示意Fig.1 Schematic diagram of the location of the study area

1.2 数据源及数据处理

Google Earth Engine为一个存储卫星数据以及针对卫星数据进行批量处理的云计算平台，它存储了全球40多年的遥感数据，包括Landsat系列、Sentinel,MODIS,NOAA AVHRR等卫星产品。平台提供Python和JavaScript两种客户端库，用户可以自由在GEE代码编辑器进行代码编写，之后传至云端进行大量数据的并行计算，提高了遥感数据的处理效率。相较于传统的本地计算模式，GEE的高存储以及云计算能力，使大量遥感数据的快速处理和保存得以实现。

本次研究数据源为GEE平台提供的1986-01-01—2019-12-30的Landsar SR数据，这些数据已经进行了几何校正，之后去除50%以上高云像素的数据，最后根据研究区位置共选取了1 460个可用数据。首先使用LandTrendr方法提取变化水体(人工挖掘水体和沉陷水体)，然后采用形态学方法提取出沉陷水体，详细过程请参考文献[40-41]，最终得到1989—2017年共29年的沉陷水体数据。选择1989—2016年数据用于沉陷水体空间位置变化研究，并将2017年沉陷水体数据作为真实值进行精度评价。

由于初始的沉陷水体数据存在细小面积斑块，且本次研究注重沉陷水体最外层边界的变化，不考虑区块内部的其他土地类型，因此选取半径为2个像素的圆形核对初始数据进行先腐蚀后膨胀的形态学开操作(图2)，用以消除细小面积斑块，并且将研究区内由于道路用地等其他土地类型产生的狭长区域一起融入积水斑块，能减少数据量，处理后的结果如图3所示。

图2 形态学开操作前后对比Fig.2 Morphological comparison before and after operation

图3 潘谢矿区1989—2016年沉陷水体分布Fig.3 Distribution of subsided water bodies in Panxie mining area from 1989 to 2016

2 研究方法

本研究以射线法为核心，提出了采煤沉陷水体方向变化自动识别方法，具体技术流程如图4所示。

图4 技术流程Fig.4 Technical flow chart

图5 射线法示意Fig.5 Schematic diagram of ray method

2.1 射线法基本原理

对于单个沉陷水体区块而言，射线法实现的基本步骤为：① 选取一个原点，并以正东方向为0°线；② 选择一定的间隔角度，沿逆时针方向依次发射射线；③ 记录原点到射线相交于该区块各年的沉陷水体边界的距离数据。如图5所示，以丁集矿某沉陷积水区3 a的变化为例，选取角度为的射线，其中内部的红色边界代表沉陷水体某区块最初产生年(2013年)的边界，次外层的绿色代表该区块2014年的边界，最外层蓝色代表该区块2015年的边界。图5中角度为的射线分别交3 a的沉陷水体边界于,，三点，记录发射原点到各交点的距离列表，即,,的距离。按照这种规则，将研究区内所有沉陷水体区块均进行射线法操作，最后得出每一区块各方向上的距离列表(列表长度=2016-+1，其中为最高产生的年份)。

2.2 射线原点的选取

在采用射线法进行沉陷水体各区块边界距离变化提取时，选取的射线原点对结果有很大的影响，因此选取适当的射线原点十分关键。由于研究区开采历史很长，缺少采煤的相关信息，无法以各矿最初的开采位置为原点进行分析，因此，对于每个沉陷水体区块而言，射线原点的选取约定2个原则：一是在进行射线法拟合时，拟合结果应最接近该区块的演化过程；二是以区块的最初产生年为起始年开始研究，能最大程度地接近矿区最初开采位置，原因是在观察研究区沉陷水体的演化过程中发现，各沉陷水体区块是由最初产生的多个小面积区块逐渐演化合并而成的。

以某一区块为例，首先以2016年最终的区块为基础，观测其是由哪些最初形成的小面积区块演化合并而成的，并以这些最初的小面积区块的中心点为射线原点，分别采用射线法构建距离列表，并对每个射线原点采用决定系数来评价，最后选取最大值，即拟合程度最好的点作为该区块的射线原点。其中决定系数的计算公式为

=

(1)

其中，为预测数据与原始数据均值之差的平方和；为原始数据与其均值之差的平方和。取值为0～1，取值越接近1，表明方程的自变量对因变量的解释能力越强，这个模型对数据拟合程度越高。和可通过以下公式计算：

(2)

(3)

2.3 射线间隔角度的选取

采用射线法对沉陷水体进行监测时，射线间隔角度与两条射线间的沉陷水体边界形状复杂性成正比，即沉陷水体边界形状越复杂，选取的射线间隔角度越小。因此需要选取合适的射线间隔角度。

在采用射线法进行沉陷水体边界拟合时，由于相邻2条射线选取的间隔角度较小，可以将射线间曲线型的沉陷水体边界用图6中的直线段来进行拟合。其中，距离差越小，说明这2条射线夹角范围内的沉陷水体边界变化程度越小，则用一条直线段进行拟合的可靠性越强。研究中发现将1°设定为射线间隔角度时，两射线间的沉陷水体边界形状如图6(a)所示，即由下面的射线至上面的射线监测过程中，沉陷水体边界一直保持向外(右)扩张，趋势不变，此时可用直线拟合沉陷水体边界；图6(b)中沉陷水体边界是先向外(右)扩张至点，再向内(左)扩张的形状，即趋势存在变化。当趋势存在多次变化时，用直线拟合沉陷水体边界精度较低，就需要在这2条射线间加密一条甚至多条射线。

图6 相邻射线间沉陷水体边界变化情况Fig.6 Boundary changes of subsided water bodies between adjacent rays

为了确定合适的射线间隔角度，首先选择面积较大的5个沉陷水体，从产生年到2016年，以1°为间隔构建射线，统计了相邻射线间的距离差，共计18 360个数据，并按从小到大的顺序排列编号。如图7所示，当距离差大于103.2 m时发生突变，故认为距离差小于103.2 m的部分是正常的，距离差大于103.2 m的部分视为异常值，即需要进一步加密射线的部分。之后共选取0.25°,0.5°,1°,2°以及4°为备选角度，分别统计其各自总体上相邻射线距离差之和以及距离差大于103.2 m(异常值)的数量，随后将0.25°和0.5°组成第1组，0.5°和1°组成第2组，以此类推共4组。以第1组为例，假设间隔角度为0.5°时，异常值的数量为，间隔角度为0.25°时，异常值的数量为，则-代表射线在由间隔角度为0.5°加密到0.25°过程中对异常值的消减程度，最后选择一个既能最大程度地减小相邻射线距离差之和又能大量消减异常值的角度。

图7 相邻射线距离差统计Fig.7 Statistics of distance difference between adjacent rays

2.4 距离变化数据与年份的相关性检验

为了确定沉陷积水区块是否持续变大，对距离变化与年份之间进行了相关性分析。可直接采用GEE平台提供的皮尔逊相关性分析(ee.Reducer.pearsonsCorrelation())，对每条射线进行相关性检验。其中，在每条射线上，将射线原点到边界的距离()和相应的年份()作为2个变量。

(4)

其中，为和的皮尔逊相关性;Cov(,)为和的协方差;Var(),Var()分别为和的方差。取值为-1～1，其中正负号分别代表正相关和负相关，的绝对值越接近1，说明和的相关性越强。其中Cov(,)可以通过以下公式计算：

(5)

Var()=Cov(,)

(6)

Var()=Cov(,)

(7)

2.5 预测沉陷水体提取

为了获得每个沉陷水体各方向的变化规律，利用射线法计算出的各沉陷水体区块各方向的距离列表数据，采用一元线性最小二乘回归方法进行沉陷水体的边界拟合及趋势分析。GEE提供了线性拟合的ee.Reducer.linearFit()功能。

=+

(8)

其中，为自变量，表示时间(区块产生年份)；为因变量，表示沉陷水体边界到原点距离的回归值；为斜率，表示沉陷水体边界到原点距离的变化趋势，>0，越大说明边界扩张速度越快；为截距，表示水体最初产生年份边界到原点距离。和可以通过以下公式计算：

(9)

(10)

在每个距离列表中，根据其对应的变化规律来预测下一年即2017年的距离，之后结合该射线的角度，计算出该角度下2017年沉陷水体边界点，将每个区块按照其所有射线方向计算出的边界点依次连接起来，作为预测的2017年沉陷水体边界结果。

2.6 沉陷水体扩张程度分析

基于得到的2017年预计沉陷水体数据和2016年沉陷水体数据，通过扩张系数来衡量各沉陷水域的扩张速度。

(11)

其中，和分别为每条射线上2016年沉陷水体边界到射线原点的距离以及预测的2017年沉陷水体边界到射线原点的距离。其中>0,且的取值越小说明沉陷水体在该方向的扩张速度越慢，沉陷水体越接近稳定。统计研究区内所有沉陷水体区块各方向射线的扩张系数，并进行分级，分级标准见表1。

表1 扩张系数分级标准Table 1 Expansion coefficient grading standard

3 结果分析

研究以1989—2016年沉陷水体为基础数据，在对其进行射线法监测时，先将沉陷水体所有区块从其产生年至2016年间所有沉陷水体边界数据按年份整合。

3.1 射线原点选取结果

至2016年，潘谢矿区各个小面积沉陷水体已经演化成面积可观的共41块沉陷水体，各区块选取的射线原点结果及相应的决定系数见表2。

表2 各沉陷水体区块选取的射线原点情况及决定系数Table 2 Ray origin situation and determination coefficient selected in each mining area

由表2可看出，41块沉陷水体每个区块的决定系数均大于0.6。统计每个区块的面积并将其作为权重，将决定系数进行加权平均操作，计算公式为

(12)

计算出加权平均决定系数为84.56%，回归方程拟合程度较好。

3.2 射线间隔角度选取结果

由上所述，选取0.25°，0.5°，1°，2°，4°为备选角度，分别统计其各自总体上相邻射线距离差之和以及距离差大于103.2 m(异常值)的数量，见表3。

表3 射线间隔角度选取Table 3 Ray interval angle selection

可见，从0.5°划分到0.25°时，相邻射线距离差之和减少了4 087.8 m，消除了7个异常值；从1°划分到0.5°时，相邻射线距离差之和减少了3 299.8 m，消除了16个异常值；从2°划分到1°时，相邻射线距离差之和减少了7 522.5 m，消除了76个异常值；从4°划分到2°时，相邻射线距离差之和减少了3 019.3 m，消除了20个异常值。因此，当射线间隔划分到1°时，相邻射线的距离差之和减少得最多，且能最大程度地消减异常值，故最终选取1°为射线间隔角度。

3.3 皮尔逊相关性检验结果

经过对沉陷水体边界至原点的距离与其相应的年份进行皮尔逊相关性分析，结果如图8所示，其中86.33%以上的射线相关性系数达到0.8以上，67.45%的相关系数达0.9以上，说明每条射线上，原点距该沉陷水体边界距离与相对应年份具有较强的相关性，之后进行的回归分析可靠性强。

3.4 预测沉陷水体提取及其对比分析

基于前面的线性拟合提取出2017年预测沉陷水体，并将其与通过遥感影像提取出的2017年沉陷水体进行叠加(图9)，统计出2者相交的面积和各自的面积，并计算精度。精度计算公式为

=

(13)

式中，为预测沉陷水体的精度；为真实2017年沉陷水体中预测正确部分的面积；为真实2017年沉陷水体的面积。

精度评价的结果见表4。其中面积较大的矿区，如谢桥矿、张集矿、顾桥矿、潘三矿和潘一矿，其预测精度均在80%以上，潘谢矿区总体预测与遥感影像提取出的相交面积为8 836.68 hm，遥感影像提取面积为10 466.56 hm，总体预测精度为84.43%。

图8 沉陷积水区边界变化与相应年份的相关性统计Fig.8 Correlation statistics between the boundary change of the subsidence water area and the corresponding year

图9 潘谢矿区2017年沉陷水体边界Fig.9 Boundary of subsided water body in Panxie mining area in 2017

表4 潘谢矿区沉陷水体预测与遥感提取面积比较Table 4 Comparison of subsidence water body prediction and real area in panxie mining area

3.5 沉陷水体扩张性分析

沉陷水体各方向的扩张性分析结果如图10所示。图10(a)中A区域整体呈现蓝色，沉陷水体扩张速度慢，沉陷水体发育较完全，稳定性强；B区域大面积呈现红色和黄色，扩张速度中等偏慢，处于较慢的发育阶段，整体较稳定；C区域扩张速度较慢，稳定性强；D区域大部分呈现棕色，处于急剧扩张状态，活跃程度高；E区域大面积呈现绿色，扩张速度较快，活跃性较强。图10(b)中A区域大面积呈现紫色和棕色，扩张速度快，活跃性强；B,C，E三个区域扩张速度均较慢，整体较稳定；D区域大部分为紫色，扩张速度快，活跃性强。图10(c)中A,C，D,E，I,J六个区域扩张速度较慢，其中C,D，E,I均处于较稳定的状态，A,J扩张速度缓慢，有较弱的稳定性；B区域从东北方向沿逆时针至西南方向沉陷水体扩张速度较慢，其余的部分扩张速度较快；F区域左边扩张速度快，右边扩张速度中等偏慢；G区域大面积呈绿色，扩张速度较快，活跃程度高；H区域处于急剧扩张状态，极不稳定，活跃性强；K区域从东北方向沿逆时针至西方向区间内，扩张速度中等偏快，剩余区间扩张速度慢。

谢桥矿、潘三矿、潘北矿、潘二矿、潘一矿建成投产时间分别为1997年、1992年、2004年、1989年、1983年，开采历史较长，且近几年开采量较小，沉陷水体扩张速度慢。张集矿、顾北矿、丁集矿、朱集矿建成投产时间分别为2001年、2007年、2007年、2007年，且其中张集矿开采时间较晚，这些矿区近几年仍处在煤炭大量开采阶段，沉陷水体扩张速度较快。顾桥矿于2007年建成投产，且处在煤炭大量开采阶段，其西北方向扩张速度较快、东南方向扩张速度慢，但图3中显示，该矿区沉陷水体是由多个平行“条带”状沉陷水体汇聚而成，射线法对连续扩张型沉陷水体的监测效果较好，而对前面这种沉陷水体的演化方式监测效果一般。

图10 潘谢矿区沉陷水体扩张速度Fig.10 Expansion speed of subsidence water body in Panxie mining area

4 结 论

(1)提出了射线法应用于采煤沉陷水体空间位置变化监测时原点和射线角度间隔确定方法。其中原点的选取是以2016年独立的沉陷水体为基础，选取出各独立沉陷水体最初形成的小面积区块，再从这些小面积区块中基于决定系数筛选出最能代表(拟合)这一独立沉陷水体演化过程的区块，将这一区块中心作为射线原点；射线间隔角度选择既能最大程度地减小相邻射线距离差之和又能大量消减异常值的角度。据此，本研究对潘谢矿区41个沉陷水体区块选取了各自的射线原点，并以1°为射线间隔，构建了共14 760条射线。

(2)通过对构建出的射线统计出的年际距离变化数据进行皮尔逊相关性分析，发现射线原点距该沉陷水体边界距离与其相对应年份具有较强的相关性，其中86.33%以上的射线相关性系数达到0.8以上，67.45%的射线相关系数达0.9以上。

(3)构建一元线性最小二乘法回归方程拟合了2017年沉陷水体边界，总体决定系数为84.56%，拟合程度良好。并将预测出的2017年沉陷水体数据与遥感影像提取的2017年沉陷水体数据进行对比，预测精度为84.43%。

(4)利用扩张系数评价了沉陷水体各方向的扩张速度。其中，谢桥矿、潘三矿、潘北矿、潘二矿、潘一矿沉陷水体扩张速度慢，张集矿、顾北矿、丁集矿、朱集矿沉陷水体扩张速度较快，顾桥矿西北方向扩张速度较快、东南方向扩张速度慢，与矿山企业的开采情况基本对应。

在缺少矿区采煤相关信息的条件下，本研究将射线法的思想引入采煤沉陷积水的研究中，提供了一个利用遥感影像数据，对采煤沉陷水体的空间位置进行拟合和预测的新思路。实现了在GEE平台上射线的自动绘制及应用，从沉陷水体的提取到射线法监测，整个研究基本实现自动化，只需使用者选取合适的遥感影像数据集、划定研究区即可。

由于影响采煤沉陷水体变化的因素不仅仅是地下采煤，比如自然降水、区域水资源调配等，对预测精度有一定影响，后续可考虑加入改进。