基于UVA与植被指数的排土场DEM模型构建关键技术

2019-02-11宫传刚卞正富卞和方雷少刚张周爱郭海桥

煤炭学报 2019年12期

宫传刚，卞正富，卞和方，雷少刚,黄赳，张周爱，郭海桥，张浩

(1.中国矿业大学国土资源研究所，江苏徐州 221116; 2.矿山生态修复教育部工程研究中心，江苏徐州 221116; 3.神华宝日希勒能源有限公司，内蒙古呼伦贝尔 021500)

随着无人机(Unmanned Aerial Vehicle,UAV)与传感器的不断发展，基于无人机平台的低空(100～1 000 m)/超低空(<100 m)摄影测量技术已显示出独特优越性。无人机摄影测量技术结合了无人机驾驶飞行器技术、遥感传感器技术、遥测遥控技术、通信技术、POS定位定姿态技术、GPS差分定位技术和遥感应用技术，具有系统灵活、起降方便、成本低、质量高、自动化和智能化等特点[6-8]，在快速获取小区域、飞行困难地区和复杂地形区域高分辨率影像方面具有明显优势[9]。目前，无人机摄影测量技术已被广泛用于地形测量[10-13]，环境监测[14-15]，农作物信息提取及应用[16-17]，森林植被信息提取[18-19]，湿地、水源地信息提取[20-21]，灾害应急救援[22-24]，电力巡视[25]，文物保护[26]，建筑物风险评估[27]等方面，其空间信息密度与精度可达机载雷达水平[28-29]。

然而，基于可见光立体像对匹配获取的地形数据大多未考虑植被与其他地物对高程的影响，因此在高精度地形提取方面难以保证垂直精度。针对这一问题，在植被类型单一且无建筑物干扰区，可以利用可见光植被指数间接反映植被高度，利用植被指数-植被高度回归模型剔除植被对地形的影响，从而获得真实地形。

以宝日希勒露天煤矿外排土场北坡为研究对象，采用旋翼无人机平台与地面GPS RTK系统，构建外排土场地形测量作业平台。建立可见光植被指数-植被高度修正模型对无人机摄影测量地形进行修正。最后用检查点对测量结果进行精度检验，同时与GPS RTK地形测量数据进行比较分析，对该测量方法的实用性和不足进行思考和研究，可为以后露天煤矿或其他复杂地形的测量提供参考。

1 研究区概况

宝日希勒露天矿位于内蒙古自治区呼伦贝尔市陈巴尔虎旗宝日希勒镇境内，地理坐标为东经:119.689°～119.761°，北纬:49.364°～49.412°。研究区位于宝日希勒露天煤矿外排土场北坡(图1)，总面积约2.2 km2，海拔为620～720 m，坡面角度为0°～30°，区内土地类型以草地与裸土地为主，植被以多裂叶荆芥(Schizonepeta multifida)、贝加尔针茅(Stipa baicalensis)和扁蓿豆(Pocokia ruthenica)为主。由于坡面起伏、土质松散，在雨水冲刷下形成多条侵蚀沟，使得坡面形态复杂。

2 数据采集及预处理

数据采集空中飞行单元采用旋翼无人机飞行平台;测量作业单元由机载数码相机、机载数据采集存储子系统组成;地面控制单元由地面站和数传电台组成;像控点测量单元、植被实测量单元由GPS RTK和钢尺组成。

图1 研究区地理位置示意Fig.1 Location map of study area

2.1 像控点、检查点和植被实测点量测

无人机上搭载的GPS精度较低，因此需对航拍添加地面像控点以及POS差分数据以保证精度[30]。此外，在研究区布置更多的控制点(检查点)，可以对解算的数据进行精度检验。像控点按照GB/T 7931—2008《1∶500，1∶1 000，1∶2 000 地形图航空摄影测量外业规范》采用内业人员设计、外业人员布设的方式进行。

干预前两组心理状态评分比较，差异无统计学意义（P＞0.05）。干预后，研究组SAS，SDS评分明显低于对照组，差异有统计学意义（P＜0.05）。见表1。

研究区共布设像控点19个，检查点8个，植被高度实测点43个(图1)。像控点均布设在无植被覆盖区且保证航向及旁向6片重叠范围内可见;检查点与植被高度实测点布置在坡顶、坡脚、平台及坡面等区域，其中1个检查点位于裸土地上、7个检查点位于植被区;植被高度实测点采用5点取样法在10 cm×10 cm样方内测得。所有点位置信息均采用GPS RTK进行测量，测量模式为“地面控制点模式”，单点测量180次取平均值。

2.2 影像获取

研究区东西长约2 400 m，南北宽约900 m，在设计航线时需根据航向重叠、旁向重叠、地面分辨率的要求来设定航高以及拍摄间隔，再由航高和相机参数来确定航线地面最低点分辨率[31]。无人机及相机等参数见表1，飞行中由GPS飞行控制系统控制相机快门进行定点曝光，相机设置为固定无穷远对焦、固定光圈以保证统一物镜畸变参数。除此之外，由于研究区地势起伏大，为减小影像畸变，重点区域采用倾斜拍摄与垂直拍摄相结合的影像获取方式，充分利用侧视影像来获取地物的侧面纹理信息，保证整体精度。实验数据采集于2017-08-30，天气晴朗无云，微风，共采集航片1 282张。

表1 无人机及飞行参数
Table 1 UAV and ground survey details

飞行参数型号及数值传感器FC300X_3.6_4000x3000(RGB)航线相对高度/m110航向重叠度/%80旁向重叠度/%75采集照片数量1 282地面像控点数量19地面检查点数量8影像分辨率/(cm·pix-1)5

2.3 数据处理

将外业采集影像与像控点导入Pix4D Mapper，在软件中对相机参数进行检校，并将每个像控点转刺到5～8张影像上。利用软件进行空中三角测量与密集匹配[32]、SfM摄影测量处理、多视图三维立体重建和几何模型建立等处理，进而生成高精度三维点云数据和数字地表模型(Digital Surface Model，DSM)，再对影像进行数字微分纠正、拼接、镶嵌等处理，可获得整个研究区数字正射影像(Digital Orthophoto Map，DOM)。其中，DSM数据为包含地面植被高度信息的地面高程模型，DOM数据为含有RGB波段的正射影像(图2)，其空间分辨率均为5 cm。

3 可见光植被指数与地形修正模型

由于Pix4D Mapper解算获得的DSM数据中高程值包含了地表植被高度，为获取研究区真实地形(DEM数据)，需建立可见光植被指数-植被高度模型修正DSM数据。

3.1 可见光植被指数-植被高度相关性分析

植被指数是地物反射波段间的不同组合方式，目前已公开发表的植被指数模型超过150种[33]，其中基于红光、近红外波段反射率差异构建的植被指数NDVI得到广泛应用。同样，绿色植被在可见光波段的光谱差异性也有大量研究，其部分常见可见光植被指数及其变形式见表2。

图2 研究区数字正射影像与数字地表模型示意Fig.2 Orthophotomap and DSM data of study area

表2 主要可见光植被指数
Table 2 Major vegetation indices of visible bands

植被指数名称计算公式绿红比值指数[34](Green-Red Ratio Index,GRRI)ρgreenρred绿蓝比值指数[35](Green-Blue Ratio Index,GBRI)ρgreenρblue归一化绿红差异指数[36](Normalized Green-Red Difference Index,NGRDI)ρgreen-ρredρgreen+ρred归一化绿蓝差异指数[37](Normalized Green-Blue Difference Index,NGBDI)ρgreen-ρblueρgreen+ρblue超绿指数[38](Excess Green,EXG)2ρgreen-ρred-ρblue红绿蓝植被指数[39](Red Green Blue Vegetation Index,RGBVI)ρ2green-(ρredρblue)ρ2green+(ρredρblue)

注:ρgreen,ρred,ρblue分别为无人机影像在绿、红、蓝波段的反射率。

表2中，为了使所有植被指数与植被高度呈正相关关系，取原文献中GRRI与GBRI植被指数的倒数。利用表2所示公式分别计算出研究区地面实测样本的6种植被指数，利用SPSS建立植被指数-植被高度相关性模型，其相关性见表3。由表3可知，研究区植被高度与植被指数均呈正相关关系，其中植被高度与RGBVI的相关性系数最高(R2=0.942)，其次为EXG(R2=0.908)，最低的为GRRI(R2=0.531)。在所列植被指数之间，GRRI与NGRDI,GBRI 与NGBDI相关性较高，绝对值均大于0.99。对比发现，与植被高度相关性较高的植被指数均由红、绿、蓝3个波段通过不同组合方式计算所得，如RGBVI和EXG。

表3 植被指数与植被高度的相关系数
Table 3 Correlation coefficient between vegetation index and vegetation height

注:**表示在置信度(双测)为0.01时显著相关。

利用植被高度实测样本(随机布置在草地和裸土地上，共43个)分为两部分，随机选择70%(30个)用于建立植被指数-植被高度回归模型，30%(13个)用于检验。利用SPSS回归分析模块拟合得到植被高度与各种植被指数之间最佳拟合曲线为线性拟合。统计回归模型的相关性系数R2和检验样本的平均绝对误差(Mean Absolute Deviation，MAD)及均方根误差(Root Mean Squared Error，RMSE)，公式如下，统计结果见表4。

(1)

(2)

式中，E为模型预测值与真实值的差值;n为验证样本数。

表4 植被指数-植被高度回归模型与精度检验
Table 4 Vegetation index-vegetation height regression model and precision test

植被指数VI模型R2MAD/mRMSE/mGRRIy=0.757x-0.6650.2740.075 0.097GBRIy=0.511x-0.4620.7670.046 0.056NGRDIy=1.796x+0.0730.2170.080 0.102NGBDIy=1.424x+0.0260.7890.040 0.053EXGy=0.007x+0.0230.9520.020 0.025RGBVIy=1.280x+0.0220.9500.022 0.026

由表4可知，不同植被指数—植被高度回归模型的R2,MAD,RMSE表现出明显差异。其中,R2最高的为EXG模型(R2=0.952)，其次为RGBVI模型(R2=0.950)，最低的为NGRDI模型(R2=0.217)。验证样本中MAD和RMSE最小值为EXG模型(MAD=0.02 m，RMSE=0.025 m)，其次为RGBVI模型(MAD=0.02 m，RMSE=0.025 m)，最大值为NGRDI模型(MAD=0.08 m，RMSE=0.102 m)。

综合分析3种评价指标及3.2节中植被指数与植被高度相关性分析结果可知，EXG模型与RGBVI模型能有效的反映研究区植被高度。下面以RGBVI模型为例估算研究区植被高度，进而得到研究区修正后数字高程模型(Digital Elevation Model，DEM)。

3.2 地形修正模型

利用RGBVI模型指数估算得到研究区内植被高度分部如图3所示。

由图3可知，研究区内裸土地(植被高度为0)面积占总面积32.7%，植被高度低于15 cm区域面积占总面积60.9%，植被高度高于15 cm区域面积占总面积6.4%。整体来看研究区东部比西部植被高度高，坡顶和平台比坡面高度高。图3中A,B，C,D为研究区4块典型区域，其植被高度估算结果细节如图4所示。

由图4可知，在研究区不同区域，植被高度对地形的影响程度与实际情况基本一致。区域A位于排土场坡顶，植被高度差异性明显，部分区域植被高度较高，对地形测量影响明显。区域B接近坡脚，内有冲蚀沟发育，地形起伏较大。区域C位于坡顶与边坡过渡区，内有防护林与少量人工建筑，植被高度差异性大。区域D位于坡脚，区内有公路、花坛等地物，其植被高度为0且轮廓非常清晰。

图3 基于RGBVI的研究区植被高度估算结果Fig.3 Estimation results of vegetation height in the study area based on RGBVI model

根据以上分析，可建立地形修正模型以除去研究区内植被高度对地形的影响。假设无人机摄影测量所获DSM数据任意点O(x,y)高程值为H0，利用植被指数-植被高度修正模型估算该点植被高度为h0，则修正后该点真实坐标为

A′(x,y,z)=A(x,y,H0-h0)

4 结果与分析

4.1 影像精度评估

为确定修正后地形数据的准确性，采用像控点误差统计和检查点验证2种方式进行精度评估。像控点误差统计见表6，检查点高程值与影像解算对应高程比较见表7。

由表6，7可知，所有像控点和检查点均有效，未出现异常值。其中像控点投影误差均控制在0.5个像素之内，精度满足设计要求;修正后影像高程大多小于修正前，点CheckP-5处为裸土地，植被指数为0，修正前后影像高程值不变。

表6 像控点误差统计
Table 6 Photo-control-points error statistics

点号高程误差/m投影误差/像素校准影像数GCP-01-0.0230.4817GCP-020.0210.4275GCP-030.0390.4435GCP-040.0190.3695GCP-05-0.0240.1455GCP-060.0260.2635GCP-070.0180.3067GCP-080.0220.2655GCP-090.0140.1945GCP-10-0.0060.4745GCP-110.0190.4385GCP-12-0.0180.4575GCP-130.0190.3155GCP-14-0.0570.3625GCP-15-0.0610.2166GCP-16-0.0770.4795GCP-17-0.0250.4778GCP-180.0210.4708GCP-190.0210.3498

表7 检查点与解算点高程比较
Table 7 Comparison between checkpointsand solution points elevation

检查点号修正前影像高程/m修正后影像高程/m检查点高程/mCheckP-1651.110651.014651.038CheckP-2625.712625.614625.591CheckP-3687.992687.771687.861CheckP-4631.266631.254631.221CheckP-5716.697716.697716.667CheckP-6697.706697.614697.631CheckP-7720.503720.411720.352CheckP-8698.755698.547698.579

使用均方根误差式(式(2))对检查点的可靠性进行评估[40]，检查点均方根误差由修正前的0.111 m变为0.045 m，精度显著提高，根据《数字航空摄影测量空中三角测量规范》(GB/T 23236—2009)的要求，测量结果符合1∶500地形数据制作要求。

4.2 基于GPS RTK的DEM数据精度验证

为进一步评估修正后无人机摄影测量地形数据的可靠性，将修正后地形数据与同期使用GPS RTK所测地形数据进行比较和分析(图5)。由于无人机摄影测量数据由点云生成，而GPS RTK地形测量数据由多点插值所得，因此无人机摄影测量(图5(a))解算数据比已有GPS RTK地形测量数据(图5(b))更加精细，在测量复杂地形时无人机摄影测量体现出明显优势。同时，将两组DEM数据相重叠(图5(c))可知，两组数据在不同区域呈现出不同高程值，下面进行具体论述。

图5 研究区DEM数据比较Fig.5 DEM data comparison in the study area

利用ArcGIS在研究区内取1 000个随机点(图5(c))，随机点处地面植被高度见表8。分别将修正后无人机解算高程数据和GPS RTK地形测量高程数据提取至随机点，建立同名点高程属性对照表，得出两组数据同名点高程差。除去异常值后，统计随机点高程差正态分布直方图如图6所示。

表8 随机点地面植被高度统计
Table 8 Vegetation height statistics of groundrandom points

植被高度/cm0(裸地)0～55～1010～20>20所占比例/%35.9024.7019.9015.603.90

由图6可知，无人机摄影测量解算高程数据与GPS RTK地形测量高程数据高程差主要集中在0.1 m内，但也有部分点的差值超过0.1 m。

图6 随机点高程差正态分布曲线Fig.6 Difference of elevation normal curve of random points

为进一步研究误差来源，在研究区中心区随机取3条剖面线(图7)进行地形对比。根据实地调查，研究区排土场以平台加斜坡的形式堆放，其中斜坡呈20°～30°，在从平台过度到斜坡时地形变化剧烈、棱角分明。由图7可知，修正后无人机解算地形与现实相符，而GPS RTK测量地形数据只能在整体上体现地形起伏。在地形平坦处两组数据高度重合，在地形突变处差异较大。

图7 无人机摄影测量解算数据与已有GPS RTK测量数据剖面对比示意Fig.7 Cross-sections through modelled surface with perpendicular distance between the surface obtained from UVA photogrammetric data and the reference surface