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

一种考虑地球曲率的高精度DEM航空重力远区地形改正方法，包括：S1、获取测点Q的位置信息。S2、依据测点Q的位置信息，从预先构建的地改高程数据库中提取经坐标转换、旋转至测点Q位于地球北极时，测点Q地改范围内所有DEM单元的高程平均值

Terrain correction method of high accuracy DEM aerogravity far area considering earth curvature

A terrain correction method of high accuracy DEM aerogravity far area considering the curvature of the earth, including S1, obtaining the position information of measuring point Q. S2. According to the location information of survey point Q, extract the average elevation values of all DEM units within the range of survey point Q when survey point q is located in the north pole of the earth through coordinate transformation and rotation from the pre built geoelevation database

【技术实现步骤摘要】
考虑地球曲率的高精度DEM航空重力远区地形改正方法
本专利技术涉及地球物理勘探
，尤其涉及一种考虑地球曲率的高精度DEM航空重力远区地形改正方法。
技术介绍
重力调查中的地形改正，是指对大地水准面至自然地形间的质量进行改正，分近中区和远区分别进行。近中区地形改进是在测点附近(半径2km以内)将地球近似为无限平板，消除高于或低于测点平面的物质盈余或亏损的影响，用平面直立方柱体公式计算；远区地形改进是在距离测点较远(半径2km以外166.7km以内)的区域，由区域重力数据库进行，在重力数据库建成前，各单位可自行改正。我国高山广布、矿产丰富，需要大量开展快速、高效的航空重力测量，但地形变化对重力值的影响非常大，需要高精度航空重力远区地形改正的方法。并且我国《区域重力调查规范》(DZ/T0082-2006)规定地形改正范围统一到166.7km，为了与地面重力对比和联合解释，航空重力的地形改正范围也应统一为166.7km。这个计算范围在不考虑地球曲率影响时计算的重力异常精度无法满足重力解释的要求，和地面观测的重力值相差较大。故需要一种考虑地球曲率的高精度DEM航空重力远区地形改正方法。目前国内重力测量中采用RGIS计算远区地形改正值，使用1km×1km的节点高程计算出5′×5′扇形块平均高程，地形改正精度约±0.2×10-5m/s2。但是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-r2-rcosθ-3cos2θ)+3sin2θ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-8cm3/(g·s2)；σij为扇形块(i,j)的平均密度，一般取2.67g/cm3；r为距地心的距离。式(1)的积分结果为：其中，E(θ,r)＝(2-r2-rcosθ-3cos2θ)+3sin2θcosθln(r-cosθ+l)，从式(2)中可以看出，如果航空重力测点位于地球北极，只要求出每个DEM在新坐标系中的经纬度和高程，即可计算出其对测点处的重力影响值。基于上述考虑地球曲率的高精度DEM航空重力远区地形改正值的直接计算方法，对实际工作中的航空重力测点Q的重力影响值进行计算。具体包括以下步骤：步骤S1、获取测点Q的位置信息和数字高程模型DEM，进行坐标转换和坐标旋转，使测点Q位于地球北极，获得DEM单元顶面中点P的新空间大地坐标。上述步骤S1包括：步骤S11、平面投影坐标系转换为空本文档来自技高网...

【技术保护点】
1.一种考虑地球曲率的高精度DEM航空重力远区地形改正方法，其特征在于，包括以下步骤：步骤S1、获取测点Q的位置信息；步骤S2、依据测点Q的位置信息，从预先构建的地改高程数据库中提取经坐标转换、旋转至测点Q位于地球北极时，测点Q地改范围内所有DEM单元的高程平均值

【技术特征摘要】
1.一种考虑地球曲率的高精度DEM航空重力远区地形改正方法，其特征在于，包括以下步骤：步骤S1、获取测点Q的位置信息；步骤S2、依据测点Q的位置信息，从预先构建的地改高程数据库中提取经坐标转换、旋转至测点Q位于地球北极时，测点Q地改范围内所有DEM单元的高程平均值步骤S3、根据所述测点Q地改范围内所有DEM单元的高程平均值结合预先设定的测点Q地形改正值计算公式，获取测点Q的地形改正值；步骤S4、根据测点Q的地形改正值，对测点Q的重力测量进行精确改正。2.根据权利要求1所述的方法，其特征在于，所述预先设定的测点Q地形改正值计算公式，包括：其中，σ为地改密度；Rq＝R+Z，R为地球平均半径，Z为测点Q的高程；E(θ,r)＝(2-r2-rcosθ-3cos2θ)+3sin2θcosθln(r-cosθ+l)，r为距地心的距离，θ为纬度；θ1、θ2代表测点Q地改范围的纬度范围。3.根据权利要求1所述的方法，其特征在于，在步骤S1之后，步骤S2之前，还包括：S11、获取数字高程模型DEM；S12、进行坐标转换和坐标旋转，使测点Q位于地球北极，获取DEM单元顶面中点P的新空间大地坐标；S13、根据DEM单元顶...

【专利技术属性】
技术研发人员：冯治汉，郭培虹，熊盛青，
申请(专利权)人：中国地质调查局西安地质调查中心，
类型：发明
国别省市：陕西,61