重力梯度数据相较于传统的重力场数据，含有更多高频信息，在油气矿产资源勘探、大地测量学、地球物理学及国防建设中具有较强的优越性和应用价值[1]．航空重力梯度测量技术有着高空间分辨率和不受地形条件限制等优势，在资源能源勘探等领域得到了广泛应用[2]．航空重力梯度测量作为一项高精尖技术，在我国的相关研究正处于起步阶段[3-4]．随着我国航空重力与重力梯度测量仪器自主研发的力度加大，建设我国自主的航空重力梯度标定场势在必行．根据公开资料，澳大利亚已在2009年完成了Kauring航空重力梯度标定场的建设[5]．美国在内华达州建立了航空重力标定场，但无法进行航空重力梯度标定[6]．我国首个民用航空重力梯度标定场的建设是国家重大科技基础设施建设项目“精密重力测量研究设施(PGMF)”的重要组成部分．PGMF的建设目标是建成国际一流、综合指标国际领先的精密重力测量研究设施，为我国地球科学基础研究及精密重力仪器研制、测量与应用研究提供必要的实验条件，满足我国地质调查、资源勘探等对重力数据和重力基准的战略需求．此次建设的航空重力梯度标定场，位于黑龙江省黑河市五大连池风景区，其范围在北纬48°30'~48°51'、东经126°00'~126°25'．标定场内有14座火山，均分布在北东向和北西向连线的交汇处，地球内部的岩浆沿着北东和北西方向两组断裂带的交汇处喷溢而出，形成排列整齐的火山锥．最近的火山喷发，则于公元1719-1721年，发生在老黑山和火烧山．此次喷发溢流的熔岩在4个地方阻塞了区内的石龙江，形成了五个火山堰塞湖，最终形成“五大连池”[7]．该区域曾经的火山活动造就了地表与地下质量的不均匀分布，从而形成明显的重力梯度效应．经前期调研和专家论证，符合航空重力梯度标定场对重力变化的要求．地形的起伏会引起重力梯度发生一定的变化，高精度地形改正对重力梯度解算具有重要的意义[8]，航空重力梯度标定场DEM数据不可或缺．重力数据处理的地形改正过程中，DEM每1 cm误差，会给重力改正带来约3 µGal的误差．参照国外重力梯度标定场建设20 µGal的要求，并适当加以提升，则须要获取标定场核心区域5 cm高程精度的数字高程模型(digital elevation model，DEM)．另外，高分辨率的正射影像可以辅助区域内重力测量点的布设与定位，因此也须获取核心区域的1∶2 000数字正射影像(digital orthophoto map，DOM)．随着低空摄影测量与遥感装备的发展，无人机及低空遥感传感器正朝着高精度、小型化的方向发展．无人机以其轻便、灵活、低空的特性，在运行成本、飞行审批、操作难度等方面具有较大的优势；同时，专为无人机设计的航空摄影测量相机与激光雷达(light detection and ranging，LiDAR)可以实现精细的对地观测[9]．在高精度定位定姿传感器的辅助下，无人机可进行厘米级分辨率与精度的航空摄影测量与机载激光扫描[10]．针对航空重力梯度标定场核心区域高精度、高分辨率DEM和DOM的获取需求，采用无人机分别搭载激光雷达与相机的方式进行地形测量．本研究从技术角度介绍航空重力梯度标定场无人机精细地形测量的实施、数据处理流程，并进行数据精度分析．1 标定场无人机地形数据获取1.1　测量区域概况测量区域位于五大连池航空重力梯度标定场核心区，分为“中心区域”测区与“近区加密”测区，测区位置如图1所示．中心区域测区为边长约5.3 km的矩形区域，面积约28.4 km²；近区加密测区为边长约2.3 km的矩形区域，面积约5.5 km²．10.13245/j.hust.220912.F001图1五大连池航空重力梯度标定场范围及测区示意图测区内多为平原，地势较为平坦，平均海拔约330 m，内有三座高约140 m的火山．其中，笔架山、卧虎山分别位于中心区域测区的东北、东南角；老黑山位于近区加密测区中部．因在秋冬季节观测，测区内大豆、玉米等农作物已经收割．测区内地物以黑土地为主，并有树林、苔藓与火山杨．在火山周围存在熔岩凝固而成的火山石海，地表较为崎岖．测区内的三座火山与典型地物如图2所示，秋冬植物的落叶与裸露的地表给机载激光扫描与地面滤波带来了便利，但黑土地的颜色与纹理重复给航空摄影测量带来了挑战，下节详述相机参数的设定．10.13245/j.hust.220912.F002图2五大连池航空重力梯度标定场景观与地物1.2　无人机地形测量任务规划航空重力梯度标定场的建设需要高精度、高分辨率的DOM和DEM作为基础空间数据．其中，高DOM须由航空摄影测量正射镶嵌而来．尽管DEM也是航空摄影测量的产品，但测区内自然地物与植被较多，航空摄影往往无法绕开植被的遮挡获取地面信息，导致地面滤波后的DEM存在较多孔洞，极大降低了数据质量．机载激光雷达对植被有穿透性，可以获取林下地面点，为地面滤波提供真实数据[11]，因此必须使用机载激光雷达获取DEM数据．本项目选用SONY RX1RII全画幅非量测型相机与Riegl VUX-1LR激光雷达，分别搭载于飞马V100固定翼无人机及飞马D20多旋翼无人机，进行测区DOM和DEM的数据获取工作．SONY RX1RII相机配备焦距35 mm的镜头，对焦固定至无穷远，像元大小4.5 μm，像幅大小为7 952×5 304像素；Riegl VUX-1LR激光雷达具有最大1 350 m测程，无限次回波能力，理论测距精度1.5 mm．图3为两个测区的影像与激光雷达航线图，表1列举了无人机正射影像与激光雷达的航线参数．10.13245/j.hust.220912.F003图3无人机观测航线图10.13245/j.hust.220912.T001表1无人机观测航线参数航线名称相对飞行高度/m重叠度航带间距/m航速/(m∙s-1)航向/%旁向/%正射影像350806512618激光雷达150—4017910因测区纬度较高，为避免太阳高度角与方位角变化对影像匹配的影响，并保证测区边缘的重叠度，摄影测量航线按照东西方向布设，正射影像航向外扩2条基线、旁向外扩1条航线．为保证正射影像精度，地面分辨率(ground sample distance，GSD)设定为4.5 cm，大于1:2 000正射影像的需求．因测区内黑土地、黑色火山岩的反射率较低，且固定翼无人机平台相对不稳定[12]，容易产生像移（运动模糊），因此航空摄影测量的相机参数需要做出调整．首先，相机的快门速度δt与像移大小正相关，且无人机平移、俯仰、滚转及偏航方向的运动产生的像移是不同的，无人机水平、俯仰与滚转、偏航方向产生的像移分别为Δpv=cvδt/(Hl)；(1)Δpw=ω'δt/[asin(l/c)]；(2)Δpk=κ'rδt/l,  r=w2+h2/2，(3)式中：c为相机主距；v为无人机飞行速度；δt为相机曝光时间；H为航高；l为像素大小；ω'为俯仰或滚转角速度；κ'为偏航角速度．设置相机曝光时间δt=0.001  s，根据实际情况，代入v=18 m/s，ω'=5 (°)/s,κ'=2 (°)/s，可得Δpv=0.40,Δpw=0.68，Δpk=0.47，均在1像素以内，符合影像质量的要求．相机的拍照模式设置为快门优先，考虑到测区地物较低的反射率，将相机ISO设置为自动模式，其上限为1 600．采用以上的航线参数，在28.4 km2的中心区域布设47条摄影测量航带，获取了6 217张影像；而5.5 km2的近区加密区域布设26条航带，获取了1 285张影像．对于激光雷达航线，根据设备的测距精度[13]与DEM精度5 cm的需求，确定航高为150 m．利用已有卫星DEM数据规划无人机仿地飞行航线，保证无人机离地高度恒定(变高航点为图3中蓝色点)，从而让测区内点密度、光斑大小与精度保持一致．因测区主要由自然地物构成，为确保有足够数量穿透植被的地面点，同时避免激光雷达的多周期回波 (multiple-time-around，MTA)混淆，激光雷达扫描频率设置为400 kHz，航速为10 m/s，旁向重叠度40%，航带间距179 m．考虑观测的准确性，截取视场角为70°的点云，则点云理论密度为31 pt/m2．由于机载惯性测量单元(inertial measurement unit，IMU)累积误差的影响[14]，直线飞行时间不应超过500 s．因此，对于边长超过5 km的中心区域，将区域平均划分为4个子区域，每个子区域内布设20条航带．为方便子区域间的拼接与航带平差，子区域之间留有300 m重叠部分．近区加密区域不分割，布设16条航带．1.3　控制点设置为了评估正射影像的绝对精度，采用防水油布制作像控点靶标，大小为1 m×1 m，其在影像上约占22像素，远大于靶标大于6倍地面分辨率的尺寸建议[15]．像控点设计与布设效果如图4所示．虽然有对辅助激光雷达航带平差的地面控制标志的研究[16]，但尚无机载激光雷达控制点靶标的成熟工程应用．对于格网DEM而言，只须验证其高程精度；因此，激光雷达检查点布设至道路、景观平台等人工平坦地物上并测量其高程，再与DEM相同位置的高程进行对比，计算而得激光雷达点云生成DEM的高程误差．这里使用标定场建设时布设的控制点，采用实时动态技术(real time kinematic，RTK)测量像控点与激光雷达检查点的坐标，参考框架为CGCS2000，使用大地高，其平面精度约为2 cm，高程精度约为3 cm．虽然测区地势较为平坦，但耕地、火山区域人员进入困难，无法绝对均匀地布设地面点，且激光雷达检查点只能布设至人工平坦地物上．在中心区域设立了13个影像像控点，26个激光雷达检查点；近区加密区域设立了12个影像像控点，16个激光雷达检查点，见图5．10.13245/j.hust.220912.F004图4像控点设计与布设效果图10.13245/j.hust.220912.F005图5影像与激光雷达控制点/检查点设置2 地形测量数据处理与成果分析2.1　正射影像数据处理使用的无人机平台带有网络RTK功能，因此首先使用厂商自带软件进行动态后处理技术(post-processing kinematic，PPK)解算摄站点外方位元素初始值．采用Agisoft Metashape软件进行包括空中三角测量(下文称“空三”)、像控点刺点、影像密集匹配、正射影像生成的摄影测量流程．首先对影像外方位元素及像控点中误差进行估计，软件据此确定空三的权．虽然影像外方位元素是由网络RTK，经PPK解算而来，为应对可能的粗差，将外方位线元素中误差定为5 m，角元素中误差定为3°；控制点由网络RTK直接测量而来，定权时放大其中误差，则定其水平中误差5 cm，垂直中误差10 cm．针对非量测相机畸变的不确定性，采用附加参数的Brown’s模型对相机内方位元素与畸变参数进行建模，提高空三的精度[17]．下一步，采用影像匹配辅助刺点方式标记各像控点靶标．对于中心区域，在像控点中均匀选取8个作为控制点、5个作为检查点；近区加密区域则为6个控制点，6个检查点，其分布如图5左列所示．正射影像处理的后续密集匹配、正射纠正与影像镶嵌拼接步骤按照Metashape软件的标准流程进行．2.2　激光雷达点云数字地形模型生成本项目激光雷达点云数据量较大，具体统计列入表 2中．因测区地形起伏及无人机不稳定的干扰，获取的点云密度与理论值有显著差异．中心区域数据平均点云密度为120 pt/m2，最大处为637 pt/m2，最小处为14 pt/m2；近区加密区域数据平均点云密度为80 pt/m2，最大处462 pt/m2，最小处为8 pt/m2．10.13245/j.hust.220912.T002表2激光雷达点云统计测区航带数航带点云数总点云数占用空间/GiB中心区域811.9×106~5.8×1072.2×10970近区加密187×106~3.5×1073.3×10815从激光雷达原始观测数据生成精细DEM须要经过点云解算、航带平差、去冗余、去噪、地面滤波等流程，如图6所示．10.13245/j.hust.220912.F006图6激光雷达点云DEM制作流程图激光雷达的原始观测记录每个点的角度与距离数据，结合机载IMU与PPK解算获取的高精度位姿信息，针对每条航带进行点云解算，得到以航带为单位的三维点云．由于相邻航带之间具有一定重叠，航带点云的同名特征之间可能会存在空间误差，即点云分层现象，这会对后续高精度DEM生成带来较大影响．航带平差是消除相邻航带间点云误差的常用手段．本项目点云数据量较大，且中心区域分成了4个子区块，因此采用多架次航带平差思路进行处理．首先，对重叠航带内的点云数据进行剖面检查，若两个重叠平面之间的距离大于0.1 m，则可认为出现点云数据分层现象．一般而言，航带间分层现象较为常见．然后，提取两条航带之间的同名特征，删除其中偏差大于5 cm的特征对，削弱粗差对航带平差的影响，进而对各区块内部航带进行平差．确保区块内部点云无分层现象后，对各区块的接边区域进行航带平差，保证区块内、区块间的点云具有较高的一致性．航带重叠区域的点云因重复扫描，会出现局部点云密度过高，即冗余现象，因此须要对点云进行去冗余处理．对于重叠区域不同航带的点云数据，若来自不同航带的临近点距离小于0.5 m，则认为存在冗余现象．从航带重叠区域的中线开始进行裁切，以达到去除冗余点的目的．去冗余不仅可以有效减少数据量，还可以去除航带边缘误差较大的点云数据．在激光雷达工作过程中，由于设备、环境和人为操作等因素的影响，往往会使得获取的三维数据出现大量的离散点，即噪声点．噪声点通常不规则的分布在测量空间内，明显高于或低于地面点，且难以通过一定的数学模型进行区分，给后续的地形、地物分离造成干扰．这里采用基于统计学的点云噪声滤除方法[18]，首先使用K-d树方法构建点云数据的空间拓扑结构，并据此查询每个点的k邻域，本文中k设定为20．遍历点云中的每个点，将点云和邻域内所有点组成点集pi，并计算点集内所有邻域点到该点的平均距离Di．若该点邻域内平均距离超出总体平均距离80倍以上，则认为该点为噪点，将其滤除．经过点云预处理步骤后，为生成DEM，还须对点云进行地面滤波，去除树木及人工地物．这里采用布料模拟滤波(cloth simulation filter，CSF)方法[19]对点云数据进行滤波处理．由于数据量较大，无法一次性将原始点云加载至内存中，须对点云进行降采样处理．本项目DEM格网为1 m，根据Nyquist采样定律，将原始点云划分至0.5 m大小的格网中，通过在每个格网内选取质心作为保留点，实现点云数据的降采样．为验证降采样地面滤波对点云精度的影响，对降采样前和降采样后的点云数据分别进行同样的CSF地面滤波处理，滤波之后再进行高程精度检核，得出结果两者之间的高程中误差在毫米级，因此降采样并不会对后续地面滤波的操作带来负面影响．CSF地面滤波效果如图7所示．10.13245/j.hust.220912.F007图7地面滤波效果图CSF地面滤波方法对于密集树林可以有效滤除树木，并显露出地面点；对于火山口与熔岩台地等崎岖地表，此方法不会影响地面点的基本形态，并滤除部分噪声．图7右列中CSF滤波在滤除树木的同时，保留了崎岖地表的特征．2.3　地形测量精度分析从重投影误差以及检查点中误差分析两个测区摄影测量的绝对精度．两个测区的像控点平均重投影误差分别为1.2和0.5像素，检查点分别为1.7和0.3像素．影像控制点与检查点的均方根误差(root-mean-square error，RMSE)如表 3所示．10.13245/j.hust.220912.T003表3影像控制点与检查点均方根误差测区点类型数量平面误差/cm高程误差/cm点位误差/cm中心区域控制点810.749.8314.56检查点516.7313.7721.67近区加密控制点69.028.8512.64检查点64.606.247.76从表3可以看出：两区域的检查点RMSE均符合1:2 000地形图平面位置中误差2.5 m，高程中误差1.2 m的精度要求．中心区域的控制点、检查点RMSE均比近区加密测区大，这可能是受到中心区域黑土地面积大，影像匹配困难的影响．图8为两个测区航空摄影测量像控点的误差图，其中误差椭圆表示误差的水平分布，而颜色代表高程误差的大小．图8表明：两个测区控制点与检查点的误差椭圆方向呈随机分布，且测区边缘点的误差相对较大，可以认为观测结果未受到系统误差的影响．10.13245/j.hust.220912.F008图8像控点误差图激光雷达检查点的误差如图9所示，检查点主要分布于道路、停车场等人工平坦地物上．将激光雷达检查点的高程与所在DEM格网的高程进行对比，可得每个点与DEM格网的差值，以不同颜色区分．两个测区的检查点误差分布并无规律，具有随机性．中心区域测区26个检查点的高程中误差为4.6 cm，最大误差为8.9 cm；近区加密测区16个检查点中误差为4.2 cm，最大误差为-7 cm．符合航空重力梯度标定场对DEM精度的需求．10.13245/j.hust.220912.F009图9激光雷达检查点误差(色标单位：m)3 结论介绍了我国首个民用航空重力梯度标定场的精细地形测量的流程与结果，包括传感器选用、任务与控制点规划、数据处理与DEM生成，以及结果精度分析．运用了无人机摄影测量、激光雷达手段获取了面积共约33 km2测区的高分辨率正射影像与高密度激光点云，并使用分块平差、CSF地面滤波的方式生成高精度DEM．结果表明：影像检查点平面、高程RMSE分别为16.73和13.77 cm；DEM高程中误差分别为4.2和4.6 cm，成果满足航空重力梯度标定场建设的分辨率与精度要求．生成的数字地形模型可为后续高精度地形改正的计算与分析提供重要的基础产品．在自然地物为主的场景下进行了大范围无人机精细地形测量，在任务规划、参数设置等方面总结了一些经验．首先，精细DEM的生成必须依靠大量地面点的直接观测，因此观测时点应选在庄稼未种植或已收割、植被较为稀疏，且未降雪的春冬季节．相较于影像，机载激光扫描具有可穿透植被的特性，以较高的点密度进行观测可以获得大量地面点．激光扫描点密度可以通过提升扫描频率与降低飞行速度实现，而这须要综合考虑激光雷达MTA混淆与飞行效率．其次，自然地物往往具有重复纹理、运动的特性，不能提供稳定的影像匹配点，因此须要提升航空摄影测量的重叠度．另外，相较于多旋翼，固定翼无人机平台稳定性较差，且相机多为硬连接，更容易受到风的影响从而降低影像质量，因此相机参数设定与其他平台与场景不同．