考虑地球曲率的高精度DEM航空重力远区地形改正方法与流程

本发明涉及地球物理勘探技术领域，尤其涉及一种考虑地球曲率的高精度dem航空重力远区地形改正方法。

背景技术：

重力调查中的地形改正，是指对大地水准面至自然地形间的质量进行改正，分近中区和远区分别进行。近中区地形改进是在测点附近(半径2km以内)将地球近似为无限平板，消除高于或低于测点平面的物质盈余或亏损的影响，用平面直立方柱体公式计算；远区地形改进是在距离测点较远(半径2km以外166.7km以内)的区域，由区域重力数据库进行，在重力数据库建成前，各单位可自行改正。

我国高山广布、矿产丰富，需要大量开展快速、高效的航空重力测量，但地形变化对重力值的影响非常大，需要高精度航空重力远区地形改正的方法。并且我国《区域重力调查规范》(dz/t0082-2006)规定地形改正范围统一到166.7km，为了与地面重力对比和联合解释，航空重力的地形改正范围也应统一为166.7km。这个计算范围在不考虑地球曲率影响时计算的重力异常精度无法满足重力解释的要求，和地面观测的重力值相差较大。故需要一种考虑地球曲率的高精度dem航空重力远区地形改正方法。

目前国内重力测量中采用rgis计算远区地形改正值，使用1km×1km的节点高程计算出5′×5′扇形块平均高程，地形改正精度约±0.2×10^-5m/s²。但是rgis软件有两个缺陷：一是不能使用高精度dem数据，难以满足目前高精度重力测量的要求；二是要求测点位于地面，无法满足航空重力测量的要求。对于国内航空重力测量仅处于试验和初步推广阶段。目前国外航空重力测量主要在平原区实施，虽然采用了高精度dem数据计算重力地改值，但地形改正范围不大，计算时也未考虑地球曲率的影响。

又由于航空重力测量获得的是海量数据，目前采用的计算方法计算速度较慢，难以满足生产需要。因此，亟需一种基于高精度dem数据、考虑地球曲率影响、快速高效的计算航空重力远区地形改正值的方法。

技术实现要素：

(一)要解决的技术问题

为了解决现有技术的上述问题，本发明提供了一种考虑地球曲率的高精度dem航空重力远区地形改正方法。大大提高重力远区地形改正精度，同时计算量小，计算速度快；长期工作可拼接不同时期完成的地改高程数据库，以更新和维护该高程数据库。

(二)技术方案

为了达到上述目的，本发明采用的主要技术方案包括：

一种考虑地球曲率的高精度dem航空重力远区地形改正方法，包括以下步骤：

步骤s1、获取测点q的位置信息。

步骤s2、依据测点q的位置信息，从预先构建的地改高程数据库中提取经坐标转换、旋转至测点q位于地球北极时，测点q地改范围内所有dem单元的高程平均值

步骤s3、根据测点q地改范围内所有dem单元的高程平均值结合预先设定的测点q地形改正值计算公式，获取测点q的地形改正值。

步骤s4、根据测点q的地形改正值，对测点q的重力测量进行精确改正。

作为本发明方法的一种改进，预先设定的测点q地形改正值计算公式，包括

其中，σ为地改密度；rq＝r+z，r为地球平均半径，z为测点q的高程；e(θ,r)＝(2-r²-rcosθ-3cos²θ)+3sin²θcosθln(r-cosθ+l)，r为距地心的距离，θ为纬度；θ1、θ2代表测点q地改范围的纬度范围。

作为本发明方法的一种改进，在步骤s1之后，步骤s2之前，还包括：

s11、获取数字高程模型dem。

s12、进行坐标转换和坐标旋转，使测点q位于地球北极，获取dem单元顶面中点p的新空间大地坐标。

s13、根据dem单元顶面中点p的新空间大地坐标，计算测点q地改范围内所有dem单元的高程平均值

s14、建立测点q位置信息与所述高程平均值关联的地改高程数据库。

作为本发明方法的一种改进，步骤s12包括：

s121、将投影坐标系描述的测点q和dem单元顶面中点p转换为空间大地坐标系描述的坐标。

s122、将空间大地坐标系描述的测点q和dem单元顶面中点p转换为空间直角坐标系描述的坐标。

s123、将空间直角坐标系旋转至测点q位于地球北极，获得新空间直角坐标系，并计算空间直角坐标系的旋转角。

s124、依据dem单元顶面中点p的空间直角坐标和空间直角坐标系的旋转角，计算dem单元顶面中点p的新空间直角坐标。

s125、依据dem单元顶面中点p的新空间直角坐标，计算dem单元顶面中点p的新空间大地坐标。

作为本发明方法的一种改进，空间直角坐标系的坐标原点为地心，z轴过地球北极。

作为本发明方法的一种改进，步骤s13中，计算测点q地改范围内所有dem单元的高程平均值包括：计算测点q经度范围为0-2π、纬度范围是距离北极弧长20-166.7km的范围内所有dem单元的高程平均值

(三)有益效果

本发明的有益效果是：

1、基于高精度dem，同时考虑考虑地球曲率对计算结果的影响，使地面重力远区地形改正精度有较大的提高。

2、地改高程数据只对应测点坐标，与测点高程无关，适合测点位于不同高程的重力测量工作。

3、绕开直接计算地改值引起的巨大计算工作量，采用首先形成或使用地改范围的地改高程数据库，通过读取该数据库中的平均高程，使用球壳扇形块计算重力远区地改值，达到快速计算的目的。

4、长期工作可拼接不同时期完成的地改高程数据库，以更新和维护该高程数据库。

附图说明

本发明借助于以下附图进行描述：

图1为本发明具体实施方式中地球近似为球体的示意图；

图2为本发明具体实施方式中球体坐标转换后椭球积分示意图；

图3为本发明具体实施方式中考虑地球曲率的高精度dem航空重力远区地形改正方法的流程图。

具体实施方式

为了更好的解释本发明，以便于理解，下面结合附图，通过具体实施方式，对本发明作详细描述。

如图1所示，将地球近似为球体，地心为球心。如图2所示，为考虑地球曲率和基于高精度dem数据，计算航空重力远区地形改正值，申请人通过坐标转换和椭球积分计算点扇形柱体对航空重力测点的重力影响值。

令球心为坐标原点，z轴过球体顶端(北极)，航空重力测点q在z轴上。把从大地水准面到地球自由表面、航空重力测点周围20～166.7km的区域，沿子午线方向以θ＝θ0,θ1,θ2,…,θn、纬线圈以划分为m×n个地壳扇形块。顶面中点为p的扇形块(i,j)的纬度、经度和高程范围为航空重力测点q的地形影响值为：

其中，r为地球平均半径，取6371.025km；z为航空重力测点q的高程；h1和h2分别为扇形块(i,j)的底面高程和顶面高程；g为引力常数6.67×10^-8cm³/(g·s²)；σij为扇形块(i,j)的平均密度，一般取2.67g/cm³；r为距地心的距离。

式(1)的积分结果为：

其中，e(θ,r)＝(2-r²-rcosθ-3cos²θ)+3sin²θcosθln(r-cosθ+l)，

从式(2)中可以看出，如果航空重力测点位于地球北极，只要求出每个dem在新坐标系中的经纬度和高程，即可计算出其对测点处的重力影响值。

基于上述考虑地球曲率的高精度dem航空重力远区地形改正值的直接计算方法，对实际工作中的航空重力测点q的重力影响值进行计算。具体包括以下步骤：

步骤s1、获取测点q的位置信息和数字高程模型dem，进行坐标转换和坐标旋转，使测点q位于地球北极，获得dem单元顶面中点p的新空间大地坐标。

上述步骤s1包括：

步骤s11、平面投影坐标系转换为空间大地坐标系：将投影坐标系描述的测点q和dem单元顶面中点p转换为空间大地坐标系描述的坐标。使得q和p的坐标都用经度纬度θ和海拔高程h表示。

这是由于在实际工作中，测点q坐标一般使用平面投影坐标，如高斯6°带投影坐标等，高精度dem单元的坐标数据有的使用平面投影坐标系描述，有的使用空间大地坐标系描述。在进行重力远区地形改正值计算时，首先要进行坐标转换，将坐标系统一，然后才能进行重力地形改正计算。本发明将高斯6°带投影坐标统一转换为空间大地坐标系。高斯投影反算公式为：

其中，l0为中央子午线经度；mf＝a(1-e²)(1-e²sin²bf)^-3/2；tf＝tanbf；bf为底点纬度，也就是当x＝x时的子午线弧长所对应的纬度，按子午线弧长公式迭代计算：

步骤s12、空间大地坐标系转换为空间直角坐标系：将以经度纬度θ、海拔高程h表示的测点q和dem顶面中点p转换为空间直角坐标系坐标(x,y,z)。使得q和p的坐标都用x、y和z表示。优选地，所述空间直角坐标系的坐标原点为地心，z轴过地球北极。

由于用大地坐标系进行三维坐标旋转难度较大，将空间大地坐标系转换为空间直角坐标系，为坐标旋转做准备。

步骤s13、将空间直角坐标系旋转至测点q位于地球北极，获得新空间直角坐标系，并计算空间直角坐标系的旋转角。简化了重力计算。包括：

将空间直角坐标系x,y,z轴分别旋转α,β,γ后得到新的空间直角坐标系x′,y′,z′，原测点q(x,y,z)在新坐标系中的坐标变为q(x′,y′,z′)。在进行坐标系旋转时，首先绕z轴旋转γ后，q(x,y,z)的坐标为(x1,y1,z)：

然后坐标系绕x轴旋转α后，q点的坐标为(x1,y′,z2)：

最后坐标系绕y轴旋转β后，q点的坐标为(x′,y′,z′)：

令x′＝y′＝0，解上述方程组(7)、(8)、(9)，获得旋转角α,β,γ分别为：

s14、依据dem单元顶面中点p的空间直角坐标和空间直角坐标系的旋转角，计算dem单元顶面中点p的新空间直角坐标。

dem顶面中点的空间直角坐标p(x,y,z)在新空间直角坐标系的坐标为p(x′,y′,z′)：

s15、依据dem单元顶面中点p的新空间直角坐标，计算dem单元顶面中点p的新空间大地坐标。

其中，为dem单元顶面中点p在新空间大地坐标系中的经度，θ′为纬度，r′为距地心的距离。

步骤s2、根据测点q地改范围内所有dem单元顶面中点p的新空间大地坐标，结合公式(2)，计算测点q的重力影响值。

上述对于测点处重力影响值的计算，虽然是基于高精度dem数据，并考虑了地球曲率的影响，但是直接计算测点的地改值，引发了巨大的计算工作量，计算速度慢，难以满足生产需求。

为此，申请人对上述考虑地球曲率的高精度dem航空重力远区地形改正值的直接计算方法中的步骤s2进行优化，优化后的步骤s2包括：

步骤s21、根据dem单元顶面中点p的新空间大地坐标，计算测点q地改范围内所有dem单元的高程平均值

高程平均值是通过每个扇形块计算的重力值累加，得出测点的远区地形改正值∑δgi,j。然后用公式(1)在地改范围内反算该地改值∑δgi,j对应的高程h2(一般取h1＝0)。这里的h2就是测点地改范围内的平均高程。由于高程平均值不是直接用扇形块的高程计算的，而是通过计算的地改值解非线性方程得到的，误差主要来源于坐标转换、积分计算和解非线性方程，计算速度快，误差可控，精度高。

步骤s22、根据测点q地改范围内所有dem单元顶面中点p的高程平均值结合式(2)的变形式(13)，计算测点q的重力影响值。

其中，σ为地改密度，一般取2.67g/cm³；rq＝r+z；θ1、θ2代表测点q地改范围的纬度范围。

上述根据高程平均值计算测点q的重力影响值，计算工作量相对较小，计算速度较快。

申请人对上述考虑地球曲率的高精度dem航空重力远区地形改正值的方法进行进一步优化，提出的航空重力远区地形改正方法，如图3所示，包括：

步骤s1、构建地改高程数据库。包括以下步骤：

步骤s11、获取测点q的位置信息和数字高程模型dem，进行坐标转换和坐标旋转，使测点q位于地球北极，获得dem单元顶面中点p的新空间大地坐标。

步骤s12、根据dem单元顶面中点p的新空间大地坐标，计算测点q地改范围内所有dem单元的高程平均值

步骤s13、建立测点q位置信息和高程平均值关联的地改高程数据库。

步骤s2、依据测点q的位置信息，从构建的地改高程数据库中提取经坐标转换、旋转至测点q位于地球北极时，测点q地改范围内所有dem单元的高程平均值

步骤s3、根据测点q地改范围内所有dem单元的高程平均值结合式(13)，计算测点q的重力影响值。

步骤s4、根据测点q的重力影响值，对测点q的重力测量进行精确改正。

优选地，步骤s12中计算测点q地改范围内所有dem单元的高程平均值包括计算测点q经度范围为0-2π、纬度范围是距离北极弧长20-166.7km的范围内所有dem单元的高程平均值

上述优化后的方法，基于高精度dem，同时考虑地球曲率对计算结果的影响；地改高程数据只对应测点坐标，与测点高程无关，适合测点位于不同高程的重力测量工作；绕开直接计算地改值引起的巨大计算工作量，采用首先形成或使用地改范围的地改高程数据库，通过读取该数据库中的平均高程，使用球壳扇形块计算重力远区地改值，达到快速计算的目的；长期工作可拼接不同时期完成的地改高程数据库，以更新和维护该高程数据库。

综上，本发明提供的方法填补了考虑地球曲率的航空重力远区地形改正的空白；改善了rgis不能使用高精度dem的弊端，使地面重力远区地形改正精度有较大的提高；通过预先计算地改数据库，在保证计算精度的前提下可大大提高地形改正计算的速度。

本申请上述方法是通过计算机程序实现的，运行主体实质上是一个控制终端。

需要理解的是，以上对本发明的具体实施例进行的描述只是为了说明本发明的技术路线和特点，其目的在于让本领域内的技术人员能够了解本发明的内容并据以实施，但本发明并不限于上述特定实施方式。凡是在本发明权利要求的范围内做出的各种变化或修饰，都应涵盖在本发明的保护范围内。