基于激光雷达数字高程模型的重力勘探地形改正方法.pdf

2023 年 10 月第 58 卷 第 5 期基于激光雷达数字高程模型的重力勘探地形改正方法邱隆君1，2，孙诚业*1，2，杨亚斌1，2，吴新刚1，2（1.中国地质科学院地球物理地球化学勘查研究所，河北廊坊 065000；2.国家现代地质勘查工程技术研究中心，河北廊坊 065000）摘要：目前基于激光雷达（LiDAR）数据分析陆地重力近区、中区地形改正方法及其影响因素的文献较少。以往有关地形改正方法的精度评价都是基于某一个数学模型的假设条件，只能讨论高程测量误差对地形改正的影响，不能分析由扇形锥、扇形柱等数学模型产生的误差。为此，采用美国地调局三维高程计划（3DEP）的 1 m1 m间距 LiDAR 高程数据，选择 12个数据区作为研究对象，以直棱柱模型重力解析式计算结果为地形改正参考值，从统计角度对比以往近区和中区地形改正数学模型与计算方法的精度主要创新点；在此基础上，分析不同计算方法、不同地形网格间距的差异，并且对比不同方法的计算效率。结果表明：近区地形改正范围与计算模型是影响近区地形改正计算的主要因素。混合模型方案可有效降低计算误差，当近区地形改正半径为 20 m 时，采用方案可以保证 K 区以外的 11个区的平均相对误差小于 20%；当近区地形改正半径为 50 m 时，采用方案可以将全区的平均相对误差控制在 13%以内，是计算精度最高的方案。地形起伏较平缓的丘陵（K区）的差异最大值与中山（J区）相当；在复杂地表的露天铜矿（L 区），7种计算方案的差异最大值处于相同量级。补角地形改正的修正公式不适用于丘陵，区域重力调查规范中的内接口补角地形改正的简化公式的适用性更强。影响中区地形改正计算结果的因素包括高程数据网格间距、中区地形改正范围和计算方法，其中网格间距是主要因素。综合计算精度和计算时间两个因素，质量线模型法为中区地形改正的首选计算方法。关键词：LiDAR数字高程模型，地形改正，网格间距，计算方法，计算精度中图分类号：P631 文献标识码：A doi：10.13810/J.cnki.issn.1000-7210.2023.05.022Gravity exploration terrain correction methodsbased on LiDARderived digital elevation modelQIU Longjun1，2，SUN Chengye1，2，YANG Yabin1，2，WU Xingang1，2（1.Institute of Geophysical and Geochemical Exploration，Chinese Academy of Geological Sciences，Langfang，Hebei 065000，China；2.National Center for Geological Exploration Technology，Langfang，Hebei 065000，China）Abstract：There are few references on terrestrial gravity nearand mediumzone terrains correction and their influencing factors based on LiDAR data systems.Previous accuracy evaluation related to terrain correction me thods is based on the assumption of a mathematical model，and it can only discuss the influence of the elevation measurement errors on terrain correction and fails to analyze the errors generated by mathematical models such as a sectorial cone or sectorial cylinder.To this end，this paper employs the 1 m1 m LiDAR elevation data from the U.S.Geological Survey threedimensional elevation program（3DEP）and selects 12 data areas as research objects.Meanwhile，it adopts the calculation results of gravity analytical formulas of the rectangular prism model as the reference value of terrain correction and statistically compares the accuracy of previous mathematical models and calculation methods for nearand mediumzone terrain correction，which is the main innovation of this paper.On this basis，this paper analyzes the differences in the calculation methods and grid spacing in various terrains and compares the calculation efficiency of different methods.The results are as follows.The correction range of the nearzone 非地震 文章编号：1000-7210（2023）05-1255-14*河北省廊坊市广阳区金光道 84号中国地质科学院地球物理地球化学勘查研究所重磁探测研究室，065000。Email：本文于 2022年 10月 17日收到，最终修改稿于 2023年 7月 18日收到。本项研究受中国地质科学院物化探所基本科研项目（现代地质勘查工程技术集成与创新）与中国地质调查局项目“柴达木盆地盐湖区物探综合调查”（DD20230298）联合资助。石 油 地 球 物 理 勘 探2023 年terrain and the computation model are the main factors affecting the calculation results，and the hybrid model can reduce the calculation error.When the correction radius of the nearzone terrain is 20 m，Scheme V can ensure that the average relative error of the 11 zones except the K zone is less than 20%.When the correction radius is 50 m，scheme VII can control the average relative error of all the zones below 13%，which indicates the highest calculation accuracy.The maximum difference in the hill area with gently undulating terrain（Zone K）is comparable to that in the middle mountain（Zone J），and in the open pit copper mine zone with complex terrain（Zone L），the order of magnitude of the maximum difference for the seven schemes is nearly the same.The correction formula for compensating angle terrain correction is not suitable for hilly areas，while the simplified formula for compensating angle terrain correction at the internal interface in the regional gravity survey specification has stronger applicability.The mediumzone terrain correction is mainly affected by three factors，including the grid spacing of elevation data，the terrain correction range for the medium zone，and calculation methods，with the most important factor being the grid spacing.By considering the calculation accuracy and calculation time，the mass line model method is the preferred choice for mediumzone terrain correction.Keywords：LiDARderived digital elevation model，terrain correction，grid spacing，calculation method，calculation accuracy邱隆君，孙诚业，杨亚斌，等.基于激光雷达数字高程模型的重力勘探地形改正方法 J.石油地球物理勘探，2023，58（5）：12551268.QIU Longjun，SUN Chengye，YANG Yabin，et al.Gravity exploration terrain correction methods based on LiDARderived digital elevation model J.Oil Geophysical Prospecting，2023，58（5）：12551268.0引言重力勘探是一种传统的地球物理勘探方法，是对实测重力数据进行自由空气校正、中间层校正、地形校正、正常场校正等预处理获得布格重力异常12，再用场分离技术和重力反演方法获得地下密度异常体的位置信息或者密度信息34 的一种间接勘探方法。在数据预处理阶段，为了消除起伏地形密度体的干扰，需要进行地形改正。它的计算过程主要考虑两个方面：第一，对起伏地形选择恰当的密度数值。既可选择能代表研究区地层特征的密度值，也可以按照有关规范设定岩石密度、海水密度、淡水密度56。第二，结合实际的地形改正范围选择合适的地形数据库。一般近区和中区地形改正范围为 02000 m，采用空间分辨率和精度较高的地形数据库，如应用高精度差分卫星定位测地技术（GPSRTK）或航摄获得的高精度数字高程模型（DEM）数据7；远区地形改正范围为 2000166700 m，需要覆盖率较高的地形数据库，如航天飞机雷达地形测绘使命（SRTM）或高级星载热发射和反射辐射仪（ASTER）得到的DEM数据。地形数据库的空间分辨率、精度、覆盖率都会影响计算结果，分析某单一变量难以清晰地论述这些关系。因此，地形数据库是全面、深入地研究地形改正的各种影响因素及其影响程度的关键要素。全站仪或GPSRTK 方式获取的高程数据精度最高，但成本高，工作效率低；航天遥感技术与空间测量分析技术的快速发展将地面测量获得DEM转变为由主动遥感方法自动生成 DEM8，其中激光雷达（LiDAR）在获取大范围、高精度、多种 DEM 方面的效率较高，明显提升了大比例尺地形测绘和勘察效率9。一般假设地形体的密度为常数，基于 DEM 地形改正的主要研究内容包括地形数据水平分辨率与计算精度之间的关系1012、不同数据源的地形改正精度要求7。由于单一源的高分辨率数据的覆盖范围有限，地形改正方法的精度评价都是基于某一个数学模型的假设条件，只能讨论高程测量误差对地形改正的影响，不能分析由扇形锥、扇形柱等数学模型近似实际地形产生的误差。目前，利用 LiDAR数据分析陆地重力近区地形改正数学模型产生的近似误差以及中区地形改正的影响因素的文献较少。本文采用美国地调局三维高程计划（3DEP）的LiDAR高程数据，选择12个数据区作为研究对象，利用 LiDAR高程数据为参考数据，以直棱柱模型重力解析式计算结果为地形改正参考值，从统计角度对比近区和中区地形改正的数学模型和计算方法的精度，这是本文的主要创新点；在此基础上，分析不同计算方法、不同地形网格间距的差异，并且对比不同方法的计算效率。1数据说明美国地质调查局从2016年至2023年实施3DEP，目前可公开访问的高程数据的最高精度为1 m1 m，1256第 58 卷 第 5 期邱隆君，等：基于激光雷达数字高程模型的重力勘探地形改正方法即栅格尺寸为1 m1 m，数据维度为1001210012，高程值为由机载LiDAR扫描获得的初始数据经滤波得到的裸露地表高程13。本文选用 3DEP 高程数据集中地形起伏特征明显的12个DEM为试验数据区，其地理位置分布见图1。在生产过程中会融合其他数据，因此最终数据产品的精度并不统一。基于地貌划分原则，试验区的地貌包括丘陵（海拔500 m，相对高度500 m）、高山（海拔3500 m，相对高度500 m），详细的统计数据见表 1，数据下载网址为：https：www.usgs.gov/3delevationprogram。3DEP提供的原始数据的水平分辨率为1 m1 m，表 1机载激光雷达 DEM 统计表试验区ABCDEFGHIJKL试验区地名Lyman Hill（喀斯喀特山脉）Mount Shasta（喀斯喀特山脉）Gannett Peak（落基山脉）Blue Mountain（落基山脉）Blanca Peak（落基山脉）Elbert Mountain（落基山脉）Longs Peak（落基山脉）Mount Catskill（阿巴拉契亚山脉）Mount Great Smoky（阿巴拉契亚山脉）Blue Ridge（阿巴拉契亚山脉）Texas Hill（德克萨斯丘陵区）Bingham Canyon Mine（宾汉峡谷铜矿区）地貌中山中山高山中山高山高山高山中山中山中山丘陵中山高程最小值m121.1714.02447.881.12522.22957.62702.8197.7605.4277.8262.31472.3高程最大值m1312.81684.74208.11305.54275.64268.84341.51257.91577.9816.9386.63081.2高程平均值m636.71177.73381.5475.72936.33603.53496.7570.31056.7444.5306.32059.4高程标准差m315.2208.6365.9292.4273.6265.5267.9213.4197.6124.736.2358.410002000 kmDEM!#NDAHLJIKGEFCB图 1DEM 试验区的地理位置分布图（背景地形图为 Stamen Terrain地形图）每个试验区对应的栅格数据的范围均为 10012 m10012 m，栅格像元尺寸均为 1 m1 m。后续数据处理中，将全部试验区数据的原始维度从 1001210012截断到 1000010000，经过粗网格处理获得较低分辨率的地形数据。1257石 油 地 球 物 理 勘 探2023 年为平面投影后数据，与 ASTER GDEM（全球 DEM）等以经纬度表示的数据不同，更适合作为重力近区和中区地形改正的试验参考数据。为分析不同网格间距对地形改正精度的影响，将 1 m1 m 栅格数据进行“粗网格”处理，处理过程的本质是“平滑”原始数据，使地形更“平缓”。以栅格数据存储 DEM 的数据空间排列都是有规律的，有些情况需要平面坐标变换，如 ASTER GDEM 高程数据，具体方法见下文。每个栅格像元的数值表示该像元覆盖面积内的平均高程，这与节点高程数据只表示单一位置高程信息的本质不同。2近区与中区地形改正方法2.1近区地形改正方法根据区域重力调查规范，1 5万比例尺近区地形改正半径为20 m，其他比例尺的近区地形改正半径为50 m1415。采用八方位圆域法计算时，内环一般采用扇形锥模型（图2a），外环采用扇形柱模型（图2b）。扇形锥模型的地形改正为 16 gc=2GRDnc(1-cos c)(1)式中：G、RD、nc、c分别为万有引力常数、近区地形改正半径、方位数、扇形锥顶面中心线与其平面投影线间的夹角；为模型密度。扇形柱模型的地形改正为gp=2Gnp(R2i+h2+Ri+1-R2i+1+h2-Ri)(2)式中：Ri 为第i环的半径；h为扇形柱相对于测点的平均高差；np为方位数。对于方域近区地形改正，重力规范中一般采用由测点 P、角点 B以及边长中点 C构成的斜顶面三角棱柱模型（图2e和图2f）计算地形改正值，模型整体与扇形锥模型类似，不同的是该模型高度主要受角点和边长中点的控制。基于该模型的方域近区地形改正为gs=j=18GRs ln(1+2)-1b2j+1ln ajbj+b2j+1+()a2j+1()b2j+1+2ajbj()b2j+1+()b2j+12ajbj+(aj 2+1)(b2j+1)(3)其中 aj=h1，jRsbj=h2，j-h1，jRs(4)式中：h1，j、h2，j（j=1，2，8）分别为第j个斜顶面三角棱柱单元的中点和角点相对于测点的高差；Rs为方域近区地形改正半径。2.2圆域数据与方域数据接口处理分环、分方位进行重力近区地形改正要将实际地RiRi+1hADCBPRm0?25.0?5.0?5.0?5.0?5.0?5.0?25.0?5.0?5.0?5.0 5.0?5.0?25.0?25.0?5.0?5.0?5.0?5.0?5.0?5.0?5.0?5.0?5.0?5.111111 0?75.0?5.0?75.11 0?5.0?5.11 0?75.0?5.0?75.111111xzyBCPCBPBCADEFGHRDRsRsh2h1(e)(f)(d)(c)(a)(b)C图 2近区和中区地形改正采用的数学模型（a）扇形锥模型；（b）扇形柱模型；（c）内接口（灰色充填部分）；（d）梯形积分系数分布；（e）斜顶面三角棱柱模型平面投影；（f）斜顶面三角棱柱模型，h1，h2分别表示方域四边中点和斜顶面三角棱柱角点相对于测点的高差。1258第 58 卷 第 5 期邱隆君，等：基于激光雷达数字高程模型的重力勘探地形改正方法形数据按照圆域数据处理，而高精度重力中区地形改正通常采用规则网格的高程数据。这些数据通常以附带平面投影坐标或者附带大地坐标信息的栅格数据存储，在平面上坐标数据通常以方域显示。在调用这些方域数据对某个测点进行中区地形改正时，还需要计算近区圆域数据与中区方域数据的间隙产生的地形改正值（图2c的灰色充填部分），这一部分质量产生的地形改正值要累加到近区圆域地形改正值中数据接口处理，在重力调查规范中也称为内接口补角计算。区域重力调查规范给出的内接口补角地形改正的简化公式为gm=G()1-0.25 Rm1.105 1-1.105Rm()hc2+()1.105Rm2 (5)式中Rm、hc分别为近区地形改正半径、补角重心高程与测点高程差，而hc=m=14wmHmm=14wm-HP(6)其中wm=1()xm-xP2+()ym-yP2(7)为权重系数。式中：xm、ym和Hm分别为补角范围内代表高程特征的点的平面坐标和高程；xP、yP和HP分别为测点的平面坐标和高程值，通常计算hc需要四个高程数值和对应的权重系数，如图 2c 的 A、B、C及 D四点的高程数据及其权重系数。由式（7）可知，wA=wC=wD=1/R2m，wB=0.5/R2m。夏江海17 给出补角地形改正的修正公式gm=2G 4 R2m+()hc2-Rm+Rmln()2+1-Rm2ln 2R2m+()hc2+Rm2R2m+()hc2-Rm+hc arcsin 2R2m+()hc22R2m+()hc2-2hc(8)2.3中区地形改正方法中区地形改正方法的实现取决于采用的数学模型。通常采用直棱柱模型解析式求解，或采用二重梯形积分快速计算，还可以采用倾斜顶面棱柱模型，利用表面积分法近似计算。理想模型数值模拟表明，高斯表面积分法的计算精度高于二重梯形积分公式。可将棱柱体质量压缩至中心垂线上质量线模型，这种简化模型的计算表达式简洁、计算速度较快。以下讨论直棱柱模型解析式的离散形式、质量线模型方法、直棱柱模型的二重梯形积分方法及基于三角面元的表面积分方法。2.3.1直棱柱模型解析式根据区域重力调查规范1415，中区地形改正采用平顶面棱柱模型计算公式，某点的水平坐标和高程为(xk，yq，hkq)，即表示x方向第k个、y方向第q个直棱柱顶面中心坐标及其高程；测点的水平坐标和高程值为(xP，yP，hP)。假设该测点进行中区地形改正的范围共包含J I个矩形棱柱，那么全部直棱柱产生的总地形改正为|gtp=Gk=0I-1q=0J-1 xlny+r()x，y，z+ylnx+r()x，y，z-z arc tan xyzr()x，y，zxP-()xk+x2xP-()xk-x2yP-()yq+y2yP-()yq-y20h(9)1259石 油 地 球 物 理 勘 探2023 年式中：(x，y，z)为相对坐标，坐标原点位于直棱柱顶面中心；r为直棱柱顶面中心到测点的相对距离；h=hkq-hP，表示直棱柱顶面中心与测点的高程差；x、y为直棱柱的水平尺寸。应用上述公式计算测点的中区地形改正时，需要排除近区范围的直棱柱模型的影响，避免重复计算。2.3.2质量线模型方法基于质量线模型的计算是将三重积分进一步简化为离散的二重积分形式，即gtm=Gk=0I-1q=0J-1()xkq()ykq()1rkq-1r2kq+h2(10)式 中：rkq为质 量 线 到 测 点 的 水 平 距 离；(x)kq、(y)kq为某一个直棱柱的水平尺寸，由于(x)kq、(y)kq在累加求和符号内部，因此能够对不同尺寸模型单元进行计算。2.3.3直棱柱模型的二重梯形积分方法利用高斯公式将地形质量引力场的三重积分转化为二重表面积分，在测点数千米范围内进行地形改正可忽略地球表面弯曲影响。由于地形体的外边界的侧面与高程基准面互相垂直，其法线与z轴垂直，因此在侧表面的积分结果为零，其余为地形高程基准面（大地水准面）以及起伏地形表面的积分结果，即gtp=Gzdxdydz()x2+y2+z232=G(Sddxdyrkq-Stdxdyr2kq+h2)(11)式中：rkq为直棱柱顶面中心到测点的水平距离；Sd表示高程基准面；St表示地形起伏表面。当用直棱柱单元划分地形时，式（11）近似为梯形积分表达式gT Gxyk=0I-1q=0J-1 Ckqrkq1-()1+()hrkq2-1(12)式中Ckq为选用的梯形积分系数，其中内节点系数为1，边缘点系数为0.5，外角点系数为0.25，内角点系数为0.75，系数具体分布见图2d。当近区和中区间的边界通过原始栅格时，需要测点附近的四个高程节点（高程数值仍为测点的高程值）的地形改正值，再进行双线性插值，从而获得测点的地形改正近似计算结果，此方法称为共用点法。2.3.4基于三角面元的表面积分方法基于散度定理，将式（11）的体积分转换为起伏地形和地形高程基准面包围区域的表面积分，这两部分通过三角面元离散，再利用高斯积分等数值积分方法估算每个三角面元的积分近似结果，最后累加求和。地形改正近似表示为18 gtp GL22n=1NEu=1NQ()WuRnu-WuRnu(13)其中Rnu=()xun2+()yun2(141)Rnu=()Rnu2+()zun2xun=x(u)1+x(u)2+x(u)3yun=y(u)1+y(u)2+y(u)3zun=z(u)1+z(u)2+z(u)3(142)式中：(x，y，z)、(x，y，z)以及(x，y，z)分别为 第n个 三 角 面 元 的 三 个 顶 点 的 空 间 坐 标；(xun，yun，zun)为 第u个 节 点 的 空 间 坐 标；(u)l(l=1，2，3；u=1，2，NQ)为节点坐标，NQ 为使用的高斯节点数目；Wu为求积权重系数；L为规则数字化网格的网格间距，L2/2为三角面元在水平面上的投影面积；NE为离散的三角面元总数。本文采用三点高斯积分公式19，节点的面积坐标分别为(0，0.5，0.5)、(0.5，0，0.5)、(0.5，0.5，0)，求 积 权 重 系 数 分 别 为1/3、1/3、1/3。2.4大地坐标高程数据的处理目前，公开访问的 ASTER GDEM、SRTM 等高程数据都采用WGS84大地基准，由理论地球位场模1260第 58 卷 第 5 期邱隆君，等：基于激光雷达数字高程模型的重力勘探地形改正方法型 EGM96推导的大地水准面为高程基准，计算栅格像元的沿经线方向尺寸为Mx，与沿纬线方向尺寸为My，即Mx=recos eLMy=ree (15)式中：e、e以及L分别为栅格像元的地心纬度、地心纬度跨度以及经度跨度；re为栅格像元到地心的距离。实际计算中，全部栅格像元水平尺寸的转换过程是一项繁琐的计算过程，由于中区范围一般在测点周围2 km以内，因此可以利用测点位置地形模型单元的水平尺寸代替测点周围全部栅格像元的水平尺寸，可极大地减小计算量。某一测点位置附近的栅格像元尺寸为 Mx()rE+NEGM96+HPcos eLMy()rE+NEGM96+HPB (16)式中：B为栅格像元纬度跨度；NEGM96为利用WGS84椭球参数和 EGM96地球重力位模型计算的大地水准面高度；e可由测点的大地纬度B换算；rE为WGS84椭球面的半径，它是e与第二偏心率e的函数，即rE=a1+e2sin2 e(17)式中 e2=(a2-b2)/b2，其中参数 a、b 分别为参考椭球的长半轴和短半轴。e与 B满足e=arctan()b2a2tan B(18)综上所述，近区地形改正方法既可采用扇形锥与扇形柱的混合模型，也可仅采用扇形锥模型，而斜顶面三角棱柱模型可避免圆域与方域数据接口问题。中区地形改正方法中直棱柱模型解析式的精度最高；梯形积分与三角面元积分都是三维体积分的二维化简形式，其中三角面元积分通过调整积分阶数提高计算精度，是较灵活的计算方法；质量线模型方法是将模型化简的方法。3数值试验3.1误差分析与数值试验环境重力勘探近区和中区地形改正的计算精度及其改善方法一直以来都是关注重点。图 3为重力中区地形改正的工作流程。从实际地形通过某种测量手段（水准测量、机载激光雷达等）获得数字高程数据，或者直接从测绘局申请数字高程数据，该过程存在测量误差（包括点位误差和高程误差）和 DEM 处理误差。在获得研究区完整 DEM 后，利用某种数学模型进行重力近区和中区地形改正，该过程会引入模型近似误差及其计算误差，这两种误差通常难以分离。本文默认模型近似误差包括计算误差，作为主要研究对象。区域重力规范1415 一般都假定密度为常数，基于某种数学模型及其公式讨论高程精度以及点位误差对重力地形改正精度的影响，如杨亚斌等20 利用扇形、锥形和扇形柱公式计算理论误差。本文以机载 LiDAR高程数据为参考，衡量不同数学模型产生的模型近似误差以及网格间距对地形改正精度的影响。假定地形密度为常数，忽略实际地层密度不均匀产生的影响，全部数值试验均在 Intel Core i711800H CPU（8核、基准频率为 2.30 GHz）、64 GB RAM的移动工作站、Windows 10 操作系统的计 算 环 境 下 完 成，采 用 基 于 MEX 文 件 的 Matlab 2020b 并行计算方法。!#$%&()%&*+,%&-./0123456DEM/0789:3;%&图 3重力中区地形改正的工作流程1261石 油 地 球 物 理 勘 探2023 年3.2近区地形改正的模型近似误差以往评价近区地形改正精度时以扇形锥、扇形柱以及三角棱柱等理论模型获得的地形改正值为参考值20，忽略数学模型近似产生的差异。本文以 1 m1 m间距的高程数据为参考数据、直棱柱模型重力解析式计算结果为参考值，评价常用数学模型的近区地形改正精度。采用八方位圆域或者方域三角棱柱计算方法，近区地形改正的范围分为020 m与050 m两种，形成7种计算方案（表2）。可见，既可以单独使用扇形锥或者三角棱柱单一模型进行近区地形改正（、与），也可以内环采用扇形锥、中环和外环采用扇形柱的混合模型方案（、与）。为模拟野外情况做如下约定：基于 1 m1 m网格数据，以间隔 45读取 1 个高程数据。如方案中采用 3环，环的外边缘到测点距离分别为 10、25和 50 m，利用计算机程序读取原始地形每个方位的半径为 10、17.5、37.5 m 圆上的高程节点值21，内环利用扇形锥、中环和外环利用扇形柱模型计算后累加求和。近区地形改正值试验设定如下：每个试验区中央 6000 m6000 m 范围内，平均分布间距为 100 m的3721(61 61)个测点，测点高程为1 m 1 m网格的高程值，12个试验区一共有44652个测点。由于近区范围为圆域，在近区分界线的直棱柱模型只有一部分处于近区范围，为保证计算严谨性，将这些模型单元细分为更小的棱柱，利用质量线模型方法计算小棱柱体在测点的地形改正值。计算每种计算方案的计算值与参考值的差异，利用该差异的最大值（10-8 m/s2）和平均相对误差进行评价。差异最大值反映了每一个试验区3721个测点在近区地形改正时由模型近似产生的最大误差；平均相对误差定义为试验区内各测点差异值的绝对值的总和与各测点参考值的总和的比值（用百分比表示），是整体评价不同方案计算精度的指标。图 4为近区地形改正不同计算方案的差异最大值对数曲线。由图可见：不同计算方案得到的差异最大值与近区地形改正半径有直接关系，差异最大值通常集中在和，即使八方位的最大差异值在某些试验区仍高于、和。单一模型方案（、和）的最大差异值较大，混合模型方案有效降低了最大差异值。差异最大值与地貌密切相关，如、和中，大于550 10-8 m/s2的差异最大值出现在高山（C区和G区）。值得注意的是，虽然K区（丘陵）的高程数据标准差只有36.2 m，但是其差异最大值与 J区（中山）相当。由于人工开采，L区形成环状阶梯式台阶，高程范围为 1472 m，3081 m，标准差为358.4 m，是相对复杂的地貌区，然而不同计算方案的差异最大值却处于一个量级（最小为13510-8 m/s2，最大为187 10-8 m/s2)，即不同方案的差异最大值变化不大。以上分析说明，不能单独以高程数据的范围和标准差判定地貌对近区地形改正的影响，如地形起伏较平缓的丘陵的差异最大值与中山相当，而地形复杂的露天铜矿区的差异最大值大致处于同一量级。图 5为近区地形改正不同计算方案的平均相对表 2近区地形改正的七种计算方案方案计算方法方域三角棱柱方域三角棱柱圆域一环圆域一环圆域二环圆域二环圆域三环地形改正范围m020050020050020050050内环范围数学模型（020 m）三角棱柱（050 m）三角棱柱（020 m）扇形锥（050 m）扇形锥(010 m)扇形锥(020 m)扇形锥(010 m)扇形锥中环范围数学模型(1025 m)扇形柱外环范围数学模型(1020 m)扇形柱(2050 m)扇形柱(2550 m)扇形柱ABCDEFGHIJKL!#$%&101001000()*+,/(10 m)-8s-2图 4近区地形改正不同计算方案的差异最大值对数曲线1262第 58 卷 第 5 期邱隆君，等：基于激光雷达数字高程模型的重力勘探地形改正方法误差曲线。由图可见：、的平均相对误差曲线的变化趋势一致，计算精度从高到低的排序为、，特别是可以保证 K区以外的 11个试验区的平均相对误差小于20%。、的平均相对误差曲线的变化趋势相近，计算精度从高到低的排序为、，其中八方位使全区的平均相对误差小于13%，计算精度最高。在相同计算方案情况下，K区的平均相对误差高于其他区；L区的平均相对误差总体上小于15%。除了差异最大值和平均相对误差外，也可以使用绝对误差值小于50 10-8 m/s2的测点数与总测点数的比值（频率）评价近区地形改正的精度。图6为A区、和的绝对误差直方图及绝对误差地形标准差散点图。由图可见：相对于，和的绝对误ABCDEFGH!#$%&010203040()*+,-%/IJKL图 5近区地形改正不同计算方案的平均相对误差曲线0 5.0 4.0 3.0 2.0 1.0!102030400#$%&/m()*/s-2(10?m)-8-200-1000100200()*/(10?ms?)-82-150-100-500150100500 5.0 4.0 3.0 2.0 1.0!102030400#$%&/m()*/s-2(10?m)-8-200-1000100200()*/(10?ms?)-82-150-100-500150100500 5.0 4.0 3.0 2.0 1.0!102030400#$%&/m()*/s-2(10?m)-8-200-1000100200()*/(10?ms?)-82-150-100-50015010050图 6A区 （上）、（中）和（下）的绝对误差直方图（左）及绝对误差地形标准差散点图（右）模型绝对误差分布直方图总测点数目为 37211263石 油 地 球 物 理 勘 探2023 年差更趋向于零误差区间分布；、和的频率分别为 0.554、0.772和 0.926，因此八方位是非常有效的计算方案，将使约 90%测点的模型绝对误差小于50 10-8 m/s2。地形标准差小于20 m的测点数占总测点数的 98.1%，地形标准差小于 10 m 的测点数占总测点数的 45.2%；由于和引入扇形柱体模型，因此误差散点逐渐向零误差线收敛、聚集。文中将、和用于 12个试验区并统计了频率（图 7）。可见，除了 C区与 G区，在其他试验区的频率都大于 0.97；与相比，明显提升了 12个试验区的频率，即频率提升最小值为 0.059、最大值为 0.539、平均值为0.274。3.3内接口补角计算方法对比根据式（5）与式（8），对 12个试验区的全部测点的东北补角进行计算，采用质量线公式计算补角地形改正的参考值。对比每个试验区3721个测点的计算结果与参考值，并统计每个试验区全部测点的差异最大值（10-8 m/s2）以及平均相对误差（图 8）。可见：式（8）降低了补角地形改正的差异最大值，幅值降低了约(2 8)10-8 m/s2(图8a)；若不计入 J 区与 K区的计算结果，由式（8）得到的平均相对误差都小于 20%（图 8b）。表 1表明，J区与 K区为地形最平缓的地区，因此补角地形改正值很小，一般都在10-11 m/s2量级。值得注意的是，K 区作为典型的丘陵区，地形相对起伏只有 36.2 m，但由式（8）得到的平均相对误差大于 90%，因此推断式（8）不适用于丘陵地区。相比而言，式（5）简洁、计算稳定性更好，适用性更强。3.4高程数据网格间距和计算方法对中区地形改正的影响重力勘探中区地形改正主要分为圆域计算或方域计算，前者利用扇形柱体模型，多用于早期手工数据成图阶段或圆域