CN111257956A - 一种基于Matlab的区域似大地水准面精化方法 - Google Patents

Abstract

本申请公开一种基于Matlab的区域似大地水准面精化方法，步骤如下：采集待测区域的重力数据G；对所述G进行粗差剔除，将经过粗差剔除后得到的重力数据进行预处理，得到格网化法耶重力异常G_pre；利用移去恢复法处理所述G_pre,得到待测区域的大地水准面N；利用KTH法处理所述G_pre，得到大地水准面N；根据高程异常和大地水准面N，完成与国家高程基准的融合，得到融合结果R₃；根据水准数据提供的高程异常和大地水准面N，完成与国家高程基准的融合，得到融合结果R₄；利用水准点对R₃和R₄进行精度验证，选择精度高的融合结果作为最终待测区域似大地水准面。本发明基于Matlab实现RCR法和KTH法的似大地水准面精化过程，实现了大地高到正常高的转换。

Description

一种基于Matlab的区域似大地水准面精化方法

技术领域

本发明涉及物理大地测量技术领域，特别是涉及一种基于Matlab的区域似大地水准面精化方法。

背景技术

随着GPS定位技术的飞速发展，可以实现在地表的任意位置测得椭球的经纬度和大地高，能够提供高分辨率的控制网。在我国采用正常高系统，以似大地水准面为基准面，地面点到似大地水准面的垂直距离被称为正常高，大地高与正常高之间的差异称为高程异常。为获得地面点高精度正常高，传统方法通常是通过高精度水准测量，即以位于青岛的国家水准原点为起算点进行每站的高差观测累加求和得到任意点的高程。水准测量工作量大且效率较低，为克服该缺点，如何利用求解大地高和高程异常的差来获取正常高获得更多的关注，GPS具有很高的效率并且不依赖观测条件，通过与IGS跟踪站联测或利用连续运行参考系统(CORS)能够快速有效地获取区域内的高分辨率大地高基准。因此，如何获取高精度高分辨率区域高程异常成为了急需解决的问题。

当下计算高程异常(大地水准面差距)的方法主要包括：地球重力场模型法、Stokes积分法、GPS高程拟合法、最小二乘配置法和卫星无线电测高法等。其中，重力场模型法和Stokes积分法是被广泛使用的方法。重力场模型法提供了全球范围内的(似)大地水准面模型，能够快速得到地球任意一点的大地水准面差距，但由于全球重力数据分布的极不均匀，导致大地水准面模型在不同区域精度差异较大。Stokes积分法理论严密，在实测重力数据的支持下能够达到很高的模型精度，但Stokes理论要求对全球重力异常积分，实际上很难实现。因此，结合二者优势，能够达到较好的效果。当下，RCR法和KTH法是应用较广泛的精化方法，但由于重力数据质量和区域地形特点的不同，两种方法计算的的结果不尽相同，精化效果常出现较大差异。

Matlab是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。该平台具有高效的数值计算及符号计算功能、完备的图形处理能力和接近数学表达式的自然化语言等优点，能够克服RCR法和KTH法计算量大的问题，能够大幅提高计算效率，非常适合应用于似大地水准面精化工作中。

如何充分利用已有计算机软件提供的快速数据处理能力，以选取最合适的精化方法进而达到最优的精化效果，以节省人力、物力及财力，是目前急需解决的技术问题。

发明内容

本发明的目的是提供一种基于Matlab的区域似大地水准面精化方法，以解决上述现有技术存在的问题，充分利用计算机软件提供的快速数据处理能力以选取最合适的精化方法进而达到最优的精化效果。

为实现上述目的，本发明提供了如下方案：本发明提供一种电解铜箔的制作方法，包括如下步骤：

S1.采集待测区域的重力数据G；对所述重力数据G进行粗差剔除，然后将经过粗差剔除后得到的重力数据进行预处理，得到格网化法耶重力异常G_pre；

S2.利用移去恢复(RCR)法对所述格网化法耶重力异常G_pre进行处理,得到待测区域的大地水准面N；利用KTH法对所述格网化法耶重力异常G_pre进行处理，得到待测区域的大地水准面N；

S3.根据水准数据提供的高程异常和所述大地水准面N，完成与国家高程基准的融合，得到融合结果R₃；根据水准数据提供的高程异常和所述大地水准面N，完成与国家高程基准的融合，得到融合结果R₄；

S4.利用水准点对所述融合结果R₃和所述融合结果R₄进行精度验证，选择精度高的融合结果作为最终的待测区域似大地水准面。

优选地，步骤S1的具体内容为：

利用重力仪测量待测区域的重力数据G，所述重力数据G包括待测区域的绝对重力值G_abs和正常重力值G_normal，利用公式ξ^G＝G_abs-G_normal求取待测区域的重力异常ξ^G；根据空间改正方法对所述重力异常ξ^G进行处理得到待测区域的空间重力异常β^G；对所述空间重力异常β^G经过层间改正、局部地形改正和均衡改正后得到待测区域的均衡重力异常ω^G，将所述均衡重力异常ω^G减去层间改正、局部地形改正和均衡改正，得到格网化法耶重力异常G_pre。

优选地，步骤S2中利用RCR法计算待测量区域的重力大地水准面的具体内容为：根据重力场模型计算大地水准面差距N_GM，计算地形间接影响N_T，计算残差大地水准面N_RES，则最终的大地水准面N＝N_GM+N_T+N_RES。

优选地，步骤S2中利用KTH法计算待测量区域的重力大地水准面的具体内容为：对大地水准面进行预估，得到大地水准面预估值

经过垂线延拓校正δN_DWC、综合地形校正

综合大气校正

和大地水准面球近似的椭球校正δN_e后得到如下式所示的重力大地水准面

本发明公开了以下技术效果：本发明获取的结果准确性高，将本发明方法应用于实际工程中，具有较高的经济价值，并且充分利用Matlab平台的优势，通过选取最合适的精化方法以达到最优的精化效果，能够节约大量的计算时间，克服由单个精化方法带来的局限性，提高模型的精度，实现了大地高到正常高的转换，节省了大量财力、物力和人力，并且针对不同区域的特点及重力数据的差异性，灵活应用最佳的精化方法，以达到最优的精化效果。运用FFT技术大幅提升了计算速度，提高了精化效率。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明基于Matlab的区域似大地水准面精化方法的流程图；

图2为本发明所测量区域RCR法的地形间接影响及残差大地水准面；

图3为本发明具体实施方式中某区域KTH法的四项改正；

图4为本发明具体实施方式中某区域RCR和KTH两种方法计算的重力大地水准面对比图；

图5为说明本发明局部地形改正所用示意图。

图3中，a代表垂线延拓校正，b代表综合地形校正,c代表综合大气校正，d代表大地水准面球近似的椭球校正。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

下面首先对本发明基于Matlab的区域似大地水准面精化方法的整体设计思想进行说明，以便于本领域相关技术人员对本发明的理解。本发明首先利用绝对重力仪对所测量区域的绝对重力值进行测量，得到高精度的重力观测数据，然后根据多种粗差探测方法对重力观测数据进行粗差剔除，并通过相关公式完成对原始数据的预处理操作，获取高精度格网法耶重力异常和空间重力异常；其次，基于Matlab平台，利用RCR法和KTH法分别编程实现重力大地水准面的计算，并结合已有的重力场模型完成大地水准面精化；再次，以GPS水准点为已知数据，利用二次曲面拟合方法完成与国家高程系统的融合；最后，通过高精度的GPS水准数据来对比两种方法的精化效果，选取最优方法，完成最终的区域似大地水准面精化。下面结合附图1-4对本发明进行详细说明。

如图1-5所示，本发明提供一种基于Matlab的区域似大地水准面精化方法，具体内容如下：

S1.重力数据采集、预处理与粗差探测

面是由点组成的，如果有足够的高程异常点，则可以组成大地水准面，但是测量的重力点总是有限的，为了用网格状的数值来表示似大地水准面上的变化状况，则需对地面上离散点的重力值进行归算。基于此思想，首先将待测区域离散化，分成一系列网格。利用重力仪获取每个网格内的绝对重力值与正常重力值，然后作差求得每个网格中计算点的重力异常，并根据空间改正技术得到其空间重力异常。其中，空间改正技术具体内容如下：

空间改正技术是将地面重力异常按照空间垂直向下延拓到大地水准面上，基于正常重力值随着距离地球椭球表面高度的增加而按照一定比例减少的理论，可以得到空间改正公式，经过空间改正后的重力异常称为空间重力异常：

g₁＝g₀+0.3086*h

其中，h是计算点正高，g₁是空间重力异常，g₀是地面重力异常。

由于重力数据在获取时总会伴随有粗差存在，计算前剔除粗差会大幅度提高模型的精度。粗差数据提出方法有很多，在数据总体精度较高且粗差较少时，可以根据粗差数据在栅格图像中显示的明显跳跃手动剔除，在粗差较多时，采用均值漂移模型粗差探测，即根据重力点周围分布点均值变化探测和剔除粗差。本实施例中利用均值漂移模型对采集的的重力数据进行粗差探测，保留精度高的数据。

然后利用层间改正、局部地形改正和均衡改正，获得高平滑度的地形均衡重力异常，以便于进行内插处理，形成平均地形均衡重力异常的基础格网数据，其中层间改正、局部地形改正和均衡改正的具体内容如下：

层间改正：空间改正在归算时是假设地面和大地水准面之间没有质量，实际情况并非如此，因此需要将此部分质量产生的重力影响移去。在计算时，将地表与大地水准面的地形看做是一个光滑的平面，且令层面的半径为无限大，称作布格片。布格片具有一定密度，一般认为该密度值为2670kg/m³。层间改正计算公式为：

g₂＝g₁-0.1119*h

式中，g₂是不完全布格异常，h是计算点正高。

局部地形改正：层间改正和局部空间改正仅假设地面相当于一个光滑的平面，大地水准面与地面点之间的质量可以用一个简单的模型代替。实际情况下，地表并不是一个光滑的平面，因而需要将平面与真实地形之间的差异考虑进去，使得结果更精确，这种方法称为地形改正。地形改正的计算方法一般有两种，一是严格积分法，二是谱方法。两种方法各有利弊，严格积分法精度较高，但是积分速度很慢；谱方法积分速度快，但在中央区的计算精度有限。考虑到二者的优缺点，采用组合法计算，即在中央区使用严格积分法，在不同地区近区的定义不同，积分半径一般为4′×6′(平坦区)、6′×8′(山区)、10′×10′(高山区)，在中央区以外区域采用谱方法。两种方法综合使用既能达到很好的精度又能提高运算效率。地形改正的计算公式通过牛顿引力公式获得，其中，严格积分法按照如下公式计算：

δg＝-Gρ(x₂lnA₁+x₁lnA₂+y₂lnB₁+y₁lnB₂+z₂arctanC₁+z₁arctanC₂)

式中，G代表万有引力常数，取6.67259×10N·m/kg，ρ是地形密度，一般取常数2670kg/m³，i，j，k取1或2，(x₁,y₁,z₁)和(x₂,y₂,z₂)代表棱柱距离中心点最邻近点和最远点的空间直角坐标。

如图5所示，中心点就是待求点，即图5里的原点o，(x₁,y₁,z₁)和(x₂,y₂,z₂)分别是离原点最近和最远的两个点的坐标。

令z₁＝0，z₂＝H-H₀＝Δh，(Δh是计算点与积分点的高差)，就是局部地形改正，公式为：

δg_TC＝-Gρ(x₂lnA₁+x₁lnA₂+y₂lnB₁+y₁lnB₂+Δh*arctanC₁)

谱方法地形改正的计算公式如下所示：

其中，G为万有引力常数，ρ为地形密度一般取2670kg/m³，Δh和l分别代表计算点和流动点之间的高差和平面距离。

均衡改正：经过空间改正、层间改正以及局部地形改正后，理论上会将不均匀的重力场去掉，使得布格异常围绕零值波动。但经过检验，布格异常在山区的计算值为负值，并且数值较大，平均每增加1000米，布格异常值就增加100毫伽，这说明大地水准面下的地壳质量有所亏损，地形的质量以某种方式被补偿，这种现象可以通过地壳均衡学说来解释。本次吉林省大地水准面采用的是应用较为广泛的Airy-Heiskanen地形均衡模型。

Airy-Heiskanen地形均衡模型认为山体“漂浮”在密度较大的底层上，山体越高沉下去的部分越深，该部分的地形密度为补偿Δρ₁＝0.6g/cm³，该项计算时根据沉下去的部分对应的亏损地形物质引力来求得计算点影响，改正值为正。沉下去的深度d可以根据漂浮条件求计算：

式中，h_p为地形高。令抵偿密度为Δρ，地形均衡改正的积分为：

δg_TC＝-GΔρ(x₂lnA₁+x₁lnA₂+y₂lnB₁+y₁lnB₂+z₂arctanC₁+z₁arctanC₂)

式中，

Δρ＝-0.6g/cm³

z₁＝-T，z₂＝-(T+4.45H)

对于海洋上的均衡改正，原理与陆地部分相同，只是在海洋区域地下质量有所溢出需要移去。其中，

Δρ＝1.64g/cm³，z₁＝-T，z₂＝-(T-2.73H)

式中，Δρ为补偿的密度差，T表示抵偿深度，一般取30公里，H为计算点正高。

通过空间改正、层间改正、局部地形改正和均衡改正后，获得了高平滑度的地形均衡重力异常，再对每个格网的地形均衡异常按地面重力归算的逆过程，分别减去层间改正、局部地形改正和均衡改正，恢复基础格网地面平均空间异常。

重力数据的格网化经过对地面重力数据的一系列改正后，实现了地面到大地水准面的归算，并且最终计算得到平滑的均衡异常，适合外推和内插。在计算大地水准面时，需要用到格网形式的数据，采取的格网化方法会对最终计算结果精度产生直接影响。数据的格网化方法分为移动拟合法、多面函数法和Shepard内插法等几种，下面分别对其进行具体介绍。

移动拟合法：移动拟合法是一种局部函数内插法，永远以待定点为中心，用它周围的已知数据定义一个函数，应用时首先将坐标原点移到待定点P，平移后数据点i的坐标为：

X_i＝x_i-x_p，Y_i＝y_i-y_p

线性移动拟合法的内插模型为：

Δg＝a+bX_i+cY_i

根据已知点的坐标,组成观测值的误差方程式，按最小二乘法求解待定系数，对待求点P，X_i＝0，Y_i＝0，故待定点的拟合值为：

Δg_p＝a

为提高精度，可根据平面距离来定权。

多面函数法：多面函数法的基本思想是在每一个数据点上，过其余所有数据点建立一个函数，通过将这些多面函数的值迭加起来，以达到最佳的表面拟合值。多面函数拟合法的表达式为：

Δg＝K₁Q(xyx_iy₁)+K₂Q(xyx₂y₂)+…+K_nQ(xyx_ny_n)

其中，K为迭加系数，Q为核函数，n为结点数。多面函数的不同核函数形式，对拟合精度影响显著，多面函数的表达式及核函数形式为：

Δg_P＝(Q_P1,Q_P2,Q_P3,…,Q_Pn)

其中，光滑因子δ＝S_ijmin为j点与距其最近的K点之间平面距离，Δx_ij、Δy_ij、Δh_ij分别为j点距i点X、Y坐标差及高差，Δg_i是第i点的均衡重力异常，i＝1,2,…,n。

Shepard内插法：Shepard内插法是以计算点为中心，取拟合半径R以内已知函数值的权中数，数据点上的权按距中心点的不同范围采用不同的权函数确定，使靠近中心点的权增大，远离中心点的权迅速减少。在Shepard局部内插模型中，选R＝0.5°，并规定：

内插函数的模型为：

式中，r_i＝((x-x_i)²+(y-y_i ²)²)^1/2，Δg是待求点均衡异常，(x_i,y_i)和(x,y)分别代表已知点和待求点的平面坐标。

S2.基于Matlab语言，运用FFT技术，利用RCR法和KTH法分别计算大地水准面。

FFT技术：二维快速傅里叶变换(2D-FFT)能够大幅提升积分速度，是离散傅里叶变换(DFT)的快速算法，基本思想是将时域转换到频域中进行计算和处理，二维快速傅里叶变换将计算的复杂度由幂级数转换为对数形式，使计算过程中的乘法次数大幅减少，尤其在需要变换的抽样点数越多时计算量节省越显著。二维函数f(x,y)的二维傅里叶变换F(u,v)的通用公式为：

相应的逆变换为：

二维函数f(x,y)和g(x,y)的卷积表达式定义为：

如此，通过对(x',y')反号使得其中一个函数关于原点旋转180°，经过位移再乘以另一个函数，最后对乘积进行积分来获得卷积积分在特定位移处的值。

利用RCR法计算重力大地水准面：地面某一点的大地水准面差距可分为三部分：

N＝N_GM+N_T+N_RES

其中,N_GM代表重力场模型分项，N_T代表地形间接影响，N_RES代表残差大地水准面。

N_GM是EGM2008重力场模型计算得到的模型大地水准面差距，

式中，GM为地心引力常数；γ为计算点的正常重力；a是参考椭球的长半径；

λ和r分别是计算点的地心纬度、经度和向径；

和

为完全规格化系数；

是完全规格化缔合Legendre函数；max＝2190。N_T为地形间接影响，按下式计算：

式中，h和h_p分别为流动点和计算点的高程，高程数据可利用SRTM数据代替。

N_RES为残差大地水准面，由剩余重力异常按stokes积分公式计算：

式中，S为stokes核函数，

是流动点坐标，

是计算点坐标，纬度值均为地心纬度。δg为剩余重力异常：

δg＝Δg_B+A_B+c-Δg_GM

式中，

式中，Δg_B是布格重力异常，A_B是布格片改正，c为地形改正，Δg_GM是模型重力异常，表达式为

斯托克斯积分时，采用格网化的平均重力异常。由以上公式可计算出大地水准面差距的三个分项，相加后求得一个大地水准面差距值，即为重力(似)大地水准面。

然后利用KTH法计算重力大地水准面：KTH法计算大地水准面，是直接利用由地面重力数据预处理后的格网空间重力异常和地球重力场模型(Global gravitationalmodel，GGM)确定大地水准面高度的预估值

再加上所有的改正项来获取最终的大地水准面差距。

其中，

是综合地形校正，包括地形对大地水准面直接和间接影响，δN_DWC是垂线延拓校正，

是综合大气校正，包括直接和间接大气影响，δN_e是大地水准面球近似的椭球校正。与RCR法不同的是，KTH法分析并利用了模型重力异常和实测重力异常的误差，计算改正参数并对Stokes核函数进行修正，可得到大地水准面预估值：

KTH法附加改正中，地形影响指地形对大地水准面的直接影响和间接影响，表达式为：

其中，H为计算点高程。

垂线延拓改正能够避免重力异常向下归算的过程，其计算公式为：

式中，

和

分别为布格片和地形的影响，表达式为：

其中，P、Q点分别为计算点和流动点，r_P＝R+H_P为计算点球半径，

为计算点大地水准面预估值。

大气改正是大气对大地水准面直接和间接影响的总和：

其中，ρ₀为海平面大气密度。

椭球改正计算了Stokes公式中大地水准面球近似带来的影响。

表达式为：

式中，

为莫洛金斯基截断系数。

第三步，根据GPS水准数据提供的高程异常和两种方法分别计算的大地水准面差距，利用二次曲面拟合实现与国家高程基准的融合。

二次曲面拟合是多项式拟合的一种形式，该方法计算简单，便于编程实现，其拟合模型为：

ΔN＝a₀+a₁x+a₂y+a₃xy+a₄x²+a₅y²

式中，ΔN代表真实值与估计值的残差；a₀～a₅为模型的待定系数；(x,y)是待定点的平面坐标。根据若干已知点列出误差方程式，利用最小二乘法求解六个待定系数，并代入上式能求得任意坐标的拟合值。二次曲面拟合法拟合的模型较平滑，依据少量的已知点，通过二次曲面拟合进行拟合就可以得到较好的效果。利用二次曲面拟合，能够实现大地水准面到似大地水准面的转换，即完成与国家高程基准的融合。

第四步，利用GPS水准点进行RCR和KTH两种精化方法的精度验证，选择最优的结果作为最终的区域似大地水准面。

在本发明的描述中，需要理解的是，术语“纵向”、“横向”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。

以上所述的实施例仅是对本发明的优选方式进行描述，并非对本发明的范围进行限定，在不脱离本发明设计精神的前提下，本领域普通技术人员对本发明的技术方案做出的各种变形和改进，均应落入本发明权利要求书确定的保护范围内。

Claims (4)

1.一种基于Matlab的区域似大地水准面精化方法，其特征在于，包括如下步骤：

S1.采集待测区域的重力数据G；对所述重力数据G进行粗差剔除，然后将经过粗差剔除后得到的重力数据进行预处理，得到格网化法耶重力异常G_pre；

S2.利用移去恢复法对所述格网化法耶重力异常G_pre进行处理,得到待测区域的大地水准面N；利用KTH法对所述格网化法耶重力异常G_pre进行处理，得到待测区域的大地水准面N；

S3.根据水准数据提供的高程异常和所述大地水准面N，完成与国家高程基准的融合，得到融合结果R₃；根据水准数据提供的高程异常和所述大地水准面N，完成与国家高程基准的融合，得到融合结果R₄；

S4.利用水准点对所述融合结果R₃和所述融合结果R₄进行精度验证，选择精度高的融合结果作为最终的待测区域似大地水准面。

2.根据权利要求1所述的基于Matlab的区域似大地水准面精化方法，其特征在于，步骤S1的具体内容为：

利用重力仪测量待测区域的重力数据G，所述重力数据G包括待测区域的绝对重力值G_abs和正常重力值G_normal，利用公式ξ^G＝G_abs-G_normal求取待测区域的重力异常ξ^G；根据空间改正方法对所述重力异常ξ^G进行处理得到待测区域的空间重力异常β^G；对所述空间重力异常β^G经过层间改正、局部地形改正和均衡改正后得到待测区域的均衡重力异常ω^G，将所述均衡重力异常ω^G减去层间改正、局部地形改正和均衡改正，得到格网化法耶重力异常G_pre。

3.根据权利要求1所述的基于Matlab的区域似大地水准面精化方法，其特征在于，步骤S2中利用RCR法计算待测量区域的重力大地水准面的具体内容为：根据重力场模型计算大地水准面差距N_GM，计算地形间接影响N_T，计算残差大地水准面N_RES，则最终的大地水准面N＝N_GM+N_T+N_RES。

4.根据权利要求1所述的基于Matlab的区域似大地水准面精化方法，其特征在于，步骤S2中利用KTH法计算待测量区域的重力大地水准面的具体内容为：对大地水准面进行预估，得到大地水准面预估值

经过垂线延拓校正δN_DWC、综合地形校正

综合大气校正

和大地水准面球近似的椭球校正δN_e后得到如下式所示的重力大地水准面

Images