0 引 言
矿区地表移动变形是复杂的时间和空间过程[1-2]，常规变形监测技术手段主要有导线测量、水准测量等 [3]。受外界因素影响，矿区地表埋设的变形监测点易受损害[4]。同时，以点的测绘来表达地表形态时，存在片面性和偶然性，难以全面地反映地表变形特征。近年来，近景测量技术[5]、合成孔径雷达差分干涉测量(D-InSAR)技术[6]、测量机器人[7]、无人机摄影测量[8]等技术在变形监测领域应用广泛，这以面状测绘代替了传统的点状测绘。目前，衡量矿区地表破坏程度仍使用倾斜变形、水平变形、曲率等参数组成的指标体系，但这些指标参数在新技术体系下的求算、提取等是亟待解决的问题。
三维激光扫描技术在工业制造、文物保护及医学可视化等领域有广泛应用[9]，在矿区开采沉陷监测方面具有独特优势，科研工作者已开展了一系列有价值的研究工作。李秋等[10]结合矿区地表沉陷监测，就其技术优势与存在问题进行讨论，总结了其应用于开采沉陷监测中的技术优势。戴华阳等[11]应用地面三维激光扫描仪开展了沉陷区房屋移动变形监测研究，提出了基于房屋特征点的采动区房屋移动变形提取方法。吴侃团队[12-14]利用三维激光扫描技术在开采沉陷监测方面展开了针对性研究，在沉陷盆地数据获取、监测精度分析、沉陷规律预测等方面取得了卓有成效的研究成果，研究表明其监测精度可达毫米级，且可通过建立的建筑物点云模型，对建筑物进行整体形变分析。在基于点的信息表示中，点所处区域的法向量是主要属性指标 [15]，除基于点的高精度、高质量绘制方法[16-17]外，许多表面重建算法[18-19]也需借助法向量来获取精确的重建结果。国内学者多利用点云数据的海量性来整体研究沉陷问题，但忽视了点云数据“点”的重要信息。如能得到正确的法向量信息，在具有噪声和外点出现的情况下，点云的几何外形仍可被有效感知[20]。基于此，基于三维激光扫描点云数据，利用法向量对大柳塔矿区地表倾斜值进行研究，建立了基于法向量的开采沉陷倾斜变形提取方法，为判定矿区地表倾斜变形程度提供了有效技术途径。
1 三维激光扫描数据滤波处理
1.1 地表沉陷数据获取
利用Leica Scanstation2激光扫描仪采集陕西大柳塔矿某工作面地表沉陷数字高程模型(DEM)数据，共2次。扫描角度为360°×270°，采用360°方式扫描，并将数据的相对坐标转换为WGS84大地坐标。
1.2 Savitzky-Golay滤波
SAVITZKY和GOLAY[21]于1964年提出Savitzky-Golay滤波(S-G滤波)，该滤波算法被广泛应用于数据流平滑除噪，是在时域内基于局域多项式最小二乘法拟合的滤波方法[22]，其原理[23]为：将原始数据中xi设为原点，取其左边nl个样本点及右边ni个样本点，构建包含nl+ni+1个采样点的窗口数组，利用P阶多项式q(n)拟合该窗口数组，见式(1)：
(1)
式中，-nl<n≤ni；P≤nl+ni+1；ak为多项式q(n)中第k次方项的系数，k为n的次数。
定义函数S表示最小二乘多项式拟合曲线的残差：
(2)
其中,q(n)表示第n点的拟合值，x(n)为第n点测量值。当式(2)取最小值时，拟合效果最好，通过移动该窗口数组可得原始数据的所有拟合点。在拟合过程中，过于偏离正常趋势曲线的噪声部分会被剔除。
1.3 大柳塔矿点云数据S-G滤波处理
三维激光扫描技术采集的数据精度高、密度大，但矿区地表存在的大量石块、树木、线杆等地物造成点云数据中存在大量杂乱点和噪声点，影响地表建模精度和质量。因此，需对点云数据进行去噪处理[24]。第1期三维激光扫描数据包含32 231个点数据，第2期包含65 537个点数据。基于MATLAB采用S-G滤波方法对点云数据进行滤波处理，滤波窗口取nl=5、ni=5，滤波前后结果对比如图1所示。由图1可知，滤波后的数据异物点已有效去除，近地表面噪声点也得到了有效抑制。滤波后的数据更加契合地表起伏变化及走势变化，在此基础上进行建模，可有效提高DEM精度。
采用Roberts梯度算子来衡量S-G滤波边缘保持能力，Roberts梯度即交叉查分的方法：
|
g(x,y)|={[g(x,y)-g(x-1,y-1)]2+
(3)
其中，g(x,y)表示在(x,y)处的变化率，|g(x,y)|表示在(x,y)处的梯度。计算滤波数据中突出的边缘区梯度，Roberts梯度由47.82降低为32.60，数据边缘梯度降低31.8%。
2 数字高程模型(DEM)表面的法向量提取
2.1 下沉盆地DEM的构建
经滤波后的点云数据可有效表达地面高程坐标，采用不规则三角网TIN(Triangulated Irregular Network)构建DEM。图2为以点云数据构建地表沉降前后的DEM模型。为更加直观地表达地表变形信息，将生成的DEM模型导入ArcScene中，建立地表三维数字模型，如图3所示。
图1 S-G滤波数据对比
Fig.1 Comparison of S-G filtering processing
图2 沉降前后数字高档模型对比
Fig.2 Comparison of digital elevation mode
由图2a知，变形前地表起伏较小，地势较为平坦；由图2b知，随地下开采地持续，地表发生了明显的沉降变形。盆地右侧下沉量大于左侧下沉量，最大下沉处位于最先开采处，随着工作面不断推进，地表下沉盆地逐渐向前发育。由图3知，下沉盆地内存在大量带状褶皱，表明地表变形在时间和空间上具有不连续性，下沉盆地分布有裂缝和凸起。为分析地表变化过程，对这两期DEM数据叠加分析，得到工作面在开采过程中的地表沉降量，如图4所示，最大沉降量可达3.703 m。
图3 数据滤波后地表三维模拟
Fig.3 3D simulation of ground surface after filtering
图4 地表沉降量
Fig.4 Amount of surface subidence
2.2 局部表面拟合法提取法向量
局部表面拟合法是利用点云数据某邻域范围内的点拟合平面，通过计算平面法向量来估算点法向量的方法，也被称为“微切平面法”[25]，其原理为：给定点集 Pi{1,2,3,…,n}(n>3)，求平面S，使得各点Pi与平面S距离的平方和最小，平面S的方程可表示：
ax+by+cz+d=0
(4)
式中：a、b、c和d为平面方程系数。
则平面S的法向量为(a,b,c)，目标函数：
(5)
式中，pi为第i个点坐标；i为点集序号。
即使得
最小，应满足：
=0 (v=a,b,c,d)
(6)
于是有
∑2(axi+byi+czi+d)xi=0
∑2(axi+byi+czi+d)yi=0
∑2(axi+byi+czi+d)zi=0
∑2(axi+byi+czi+d)=0
(7)
式中，xi、yi和zi分别为pi点坐标分量，m。
解得：
(8)
可求得平面S的法向量F=(a,b,c)，即为通过局部平面拟合法求得的待求点Pi的法向量。
3 基于法向量变化的倾斜变形值计算
在传统的开采沉陷地表变形评价指标体系中，倾斜是重要指标之一，由相邻2点的下沉差及水平距离计算所得，但在面状测绘技术中较难求算。为定量计算下沉盆地各区域倾斜值，选取8个不同位置样本点，如图5所示，点A、B位于下沉盆地的边界区域，点C、G、H位于盆地底部区域，D、E、F三点为于盆地裂缝或凸起褶皱区域。
图5 样本点选取
Fig.5 Selection sample point
设各点法向量的模为1，利用地表观测的2期数据，求得8个样本点的单位法向量ei=(ai,bi,ci)及其变化量，见表1。
为验证上述方法的有效性，在2期DEM上提取8个点的坐标(Xi,Yi,Zi)，并利用各点坐标变化量求得其倾斜角，见表2。
表1 样本点法向量变化量
Table 1 Angle of sample point normal vector
样本点第1期单位法向量第2期单位法向量abcabc法向量变化量/(°)A0.113 9-0.183 80.999 9-0.183 10.102 40.987 023.42B0.195 8-0.242 50.899 9-0.202 20.232 20.879 238.47C-0.001 4-0.002 20.999 90.220 1-0.004 80.999 718.21D0.200 9-0.236 30.999 9-0.131 30.243 80.998 531.14E0.104 50.104 70.999 9-0.104 0-0.231 30.994 022.38F0.100 7-0.182 30.999 9-0.157 40.117 10.998 222.36G-0.181 60.205 30.989 90.162 5-0.276 70.995 036.45H0.100 40.171 60.999 9-0.115 0-0.138 80.999 121.41
表2 样本点下沉倾角
Table 2 Dip angle of sample point
样本点样本点坐标/m相邻点坐标/mXYZXYZ倾斜角/(°)A46 074.31251 673.2311 084.40746 076.25851 675.0041 083.25723.60B46 066.42151 538.1291 082.85546 067.46151 539.1841 081.69838.01C46 025.24551 523.3491 082.08146 026.89151 526.9341 080.82717.64D46 015.51251 473.8741 082.90046 018.82451 473.5821 080.79632.39E45 885.23151 528.5461 082.98445 886.42251 529.4881 082.33123.26F45 925.47251 618.3891 083.35545 924.89551 621.6411 082.03821.74G45 995.34251 628.1241 083.83545 997.48051 630.2661 081.50837.55H45 960.34251 703.2131 083.70945 958.70751 706.0271 082.47520.76
对法向量变化量与倾斜角做回归分析。相关性系数(Pearson Correlation Coefficient，r)是反映模型相关性的重要统计量。r值在-1到1之间，其数值大小反映了回归贡献的相对程度。
(9)
其中，Cov(X,Y)为X与Y的协方差，Var[X]为X的方差，Var[Y]为Y的方差。计算得到法向量夹角与真实下沉倾角相关性系数为0.97，验证了利用法向量夹角研究下沉倾角的科学性。为评价法向量提取下沉倾角方法的可靠性，计算法向量夹角与真实下沉倾角偏差如图6所示，其中误差为0.78°(表3)能满足实际工作需求。
图6 法向量夹角与下沉倾角对比
Fig.6 Comparison between angle of normal vector and dip angle
表3 精度分析
Table 3 Precision analysis
样本点ABCDEFGH误差/(°)-0.180.460.57-1.25-0.880.62-1.10.65中误差/(°)0.78
4 结 论
1)利用三维激光扫描等面状测绘技术，基于局部表面拟合法提取地面模型法向量，提出了矿区地表倾斜值的提取技术方法，解决了在面状测绘技术背景下地表倾斜值的有效、精确提取问题，为判定矿区地表倾斜变形程度提供了有效技术方法，对InSAR、无人机测绘等技术应用于矿区变形监测有一定的理论价值。
2)利用Savitzky-Golay滤波数据处理技术，有效剔除非地面噪声和异常地物噪声，数据边缘梯度降低31.8%，实现了基于三维激光扫描技术的地表点云模型的精确建立，为点云模型法向量的精确提取建立了保障。
3)研究并应用了基于三维激光扫描技术的倾斜值提取方法，基于地表沉降前后2期点云数据，通过计算位于下沉盆地内不同区域特征点的法向量变化量，求取了地表倾斜值，并基于实例验证了其精度，中误差达0.78°，能够满足矿区地表变形监测的实际工作需求。
[1] 崔希民,邓喀中.煤矿开采沉陷预计理论与方法研究评述[J].煤炭科学技术, 2017, 45(1):160-169.
CUI Ximin, DENG Kazhong. Research review of predicting theory and method for coal mining subsidence [J]. Coal Science and Technology, 2017, 45(1):160-169.
[2] 余学义,郭文彬,赵兵朝,等. 厚黄土层煤层开采沉陷规律研究[J].煤炭科学技术, 2015, 43(7):6-10.
YU Xueyi, GUO Wenbin, ZHAO Bingchao, et al. Study on mining subsidence law of coal seam with thick overlying loess stratum[J]. Coal Science and Technology, 2015, 43(7):6-10.
[3] 施展宇.地面三维激光扫描技术在开采沉陷应用研究[D]. 西安科技大学, 2014.
[4] 胡大贺,吴 侃,陈冉丽.三维激光扫描用于开采沉陷监测研究[J].煤矿开采,2013(1):20-22,35.
HU Dahe, WU Kan, CHEN Ranli. Data disposal of mining subsidence monitoring by 3-d laser scanning[J]. Coal Mining Technology,2013(1):20-22,35.
[5] 杨化超,邓喀中,张书毕,等.数字近景摄影测量技术在矿山地表沉陷监测中的应用研究[J].中国图象图形学报,2008, 13(3):519-524.
YANG Huachao, DENG Kazhong, ZHANG Shubi, et al.Research on subsidence monitoring of mine surface using digital close-range photogrammetric technique[J]. Journal of Image and Graphics, 2008, 13(3):519-524.
[6] 陈炳乾,邓喀中,范洪冬.基于D-InSAR技术和SVR算法的开采沉陷监测与预计[J].中国矿业大学学报,2014, 43(5):880-886.
CHEN Bingqian,DENG Kazhong, FAN Hongdong. Combining D-InSAR and SVR for monitoring and prediction of mining subsidence[J]. Journal of China University of Mining and Technology, 2014, 43(5):880-886.
[7] 杨晓玉.测量机器人在矿区地表动态变形监测中的应用[J].煤矿安全,2016,47(5):137-140.
YANG Xiaoyu. Application of measurement robot in surface dynamic deformation monitoring of mining area[J]. Safety in Coal Mines,2016,47(5):137-140.
[8] 田帅帅,赵艳玲,李亚龙,等.高潜水位矿区采煤沉陷地DEM的无人机构建方法[J].测绘通报,2018(3):98-101.
TIAN Shuaishuai, ZHAO Yanling, LI Yalong, et al.DEM establishing method of mining subsidence in high underground water mining area with UAV[J]. Bulletin of Surveying and Mapping,2018(3):98-101.
[9] 党晓斌. 三维激光扫描技术在建筑物形变监测中的应用研究[D]. 西安：长安大学, 2011.
[10] 李 秋,秦永智,李宏英.激光三维扫描技术在矿区地表沉陷监测中的应用研究[J].煤炭工程,2006(4):97-99.
LI Qiu, QIN Yongzhi, LI Hongying. Research and application of laser 3D scan technology to monitoring and measurement of mine ground subsidence[J]. Coal Engineering, 2006(4):97-99.
[11] 戴华阳,廉旭刚,陈 炎,等.三维激光扫描技术在采动区房屋 变形监测中的应用[J].测绘通报,2011(11):44-46.
DAI Huayang, KANG Xugang, CHEN Yan, et al. Study of the deformation of houses induce by mining based on 3D laser scanning[J]. Bulletin of Surveying and Mapping,2011(11):44-46.
[12] 蔡来良, 吴 侃, 张 舒.点云平面拟合在三维激光扫描仪变形监测中的应用[J]. 测绘科学, 2010, 35(5):231-232.
CAI Lailiang, WU Kan, ZHANG Shu. Application of point cloud plan fitting to deformation monitoring using 3D laser scanner[J]. Science of Surveying and Mapping, 2010, 35(5):231-232.
[13] 张 舒,吴 侃,王响雷,等.三维激光扫描技术在沉陷监测中应用问题探讨[J].煤炭科学技术,2008,36(11):92-95.
ZAHNG Shu, WU Kan, WANG Xianglei, et al. Discussion on application of 3D laser scanning technology to ground subsidence monitoring[J]. Coal Science and Technology, 2008,36(11):92-95.
[14] 陈冉丽,吴 侃.三维激光扫描用于获取开采沉陷盆地研究[J].测绘工程,2012,21(3):67-70.
CEHN Ranli, WU Kan. Research and application of 3D laser scanning to subsidence basin[J].Engineering of Surveying and Mapping，2012,21(3):67-70.
[15] 李 泉,程效军.自定位手持式三维激光扫描仪精度测试与分析[J].测绘通报, 2016(10):65-68.
LI Quan,CHENG Xiaojun. Investigating hand-held laser scanner accuracy [J]. Bulletin of Surveying and Mapping, 2016(10):65-68.
[16] Rusinkiewicz S,Levoy M. QSplat:a multiresolution point rendering system for large meshes[C]// In Proceedings of SIGGRAPH 2000,2000,343-352.
[17] Zwicker M, Pfister H, Baar J V, et al.Surface splatting[C]// Conference on Computer Graphics and Interactive Techniques. ACM, 2001:371-378.
[18] OHTAKE Y, BELYAEV A, ALEXA M, et al. Multi-level partition of unity implicits[J]. Acm Transactions on Graphics, 2003, 22(3):463-470.
[19] SHEN C, O′BRIEN J F, SHEWCHUK J R. Interpolating and approximating implicit surfaces from polygon soup[C]// ACM SIGGRAPH. ACM, 2005:204.
[20] 李 宝,程志全,党 岗,等.三维点云法向量估计综述[J]. 计算机工程与应用, 2010, 46(23):1-7.
LI Bao, CHENG Zhiquan, DANG Gang, et al.Survey on normal estimation for 3D point clouds[J], Computer Engineering and Applications, 2010, 46(23):1-7.
[21] SAVITZKY A, GOLAY M J E. Smoothing and differentiation of data by simplified least squares procedures[J]. Analytical Chemistry, 1964, 36(8):1627-1639.
[22] 赵安新,汤晓君,张钟华,等. 优化Savitzky-Golay滤波器的参数及其在傅里叶变换红外气体光谱数据平滑预处理中的应用[J].光谱学与光谱分析, 2016, 36(5):1340-1344.
ZHAO Anxin, TANG Xiaojun, ZHANG Zhonghua, et al.Optimizing savitzky-golay parameters and its smoothing pretreatment for ftir gas spectra[J]. Spectroscopy and Spectral Analysis, 2016, 36(5):1340-1344.
[23] 朱红运,王长龙,王建斌,等.基于奇异值分解和Savitzky-Golay滤波器的信号降噪方法[J].计算机应用, 2015, 35(10):3004-3007.
ZHU Hongyun, WANG Changlong, WANG Jianbin, et al.Signal denoising method based on singular value decomposition and Savitzky-Golay filter[J]. Journal of Computer Applications, 2015, 35(10):3004-3007.
[24] 敖建锋.动态沉陷区地面激光扫描数据处理关键问题研究[D].徐州：中国矿业大学, 2013.
[25] 王兆丰,闫 镔,童 莉,等. 自适应邻域尺寸选择的点云法向量估计算法[J].红外与激光工程, 2014, 43(4):1322-1326.
WANG Zhaofeng, YAN Bin, TONG Li, et al.Normal estimate method of point clouds based on adaptive neighbor size[J]. Infrared and Laser Engineering, 2014, 43(4):1322-1326.
