CN108873091B - 卫星重力梯度全张量恢复地球重力场的确定方法及系统 - Google Patents
卫星重力梯度全张量恢复地球重力场的确定方法及系统 Download PDFInfo
- Publication number
- CN108873091B CN108873091B CN201810695075.4A CN201810695075A CN108873091B CN 108873091 B CN108873091 B CN 108873091B CN 201810695075 A CN201810695075 A CN 201810695075A CN 108873091 B CN108873091 B CN 108873091B
- Authority
- CN
- China
- Prior art keywords
- grid
- gravity gradient
- full tensor
- disturbance
- data
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 230000005484 gravity Effects 0.000 title claims abstract description 232
- 238000000034 method Methods 0.000 title claims abstract description 41
- 238000012545 processing Methods 0.000 claims abstract description 37
- 230000014509 gene expression Effects 0.000 claims abstract description 33
- 238000011084 recovery Methods 0.000 claims abstract description 18
- 238000004364 calculation method Methods 0.000 claims description 36
- 238000005259 measurement Methods 0.000 claims description 22
- 239000011159 matrix material Substances 0.000 claims description 18
- 238000007781 pre-processing Methods 0.000 claims description 11
- 238000012937 correction Methods 0.000 claims description 9
- 230000008859 change Effects 0.000 claims description 8
- 230000000295 complement effect Effects 0.000 claims description 7
- 238000009499 grossing Methods 0.000 claims description 5
- MYVIATVLJGTBFV-UHFFFAOYSA-M thiamine(1+) chloride Chemical compound [Cl-].CC1=C(CCO)SC=[N+]1CC1=CN=C(C)N=C1N MYVIATVLJGTBFV-UHFFFAOYSA-M 0.000 claims description 4
- 230000008030 elimination Effects 0.000 claims description 3
- 238000003379 elimination reaction Methods 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 3
- 230000008569 process Effects 0.000 abstract description 6
- 230000006870 function Effects 0.000 description 39
- 230000010354 integration Effects 0.000 description 10
- 238000004422 calculation algorithm Methods 0.000 description 9
- 238000009795 derivation Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V7/00—Measuring gravitational fields or waves; Gravimetric prospecting or detecting
- G01V7/02—Details
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V7/00—Measuring gravitational fields or waves; Gravimetric prospecting or detecting
Landscapes
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Complex Calculations (AREA)
Abstract
本发明公开一种卫星重力梯度全张量恢复地球重力场的确定方法及确定系统。确定方法包括：对卫星重力梯度测量数据进行预处理，获得沿轨扰动引力梯度数据；对沿轨扰动引力梯度数据进行格网化处理，获得格网点值扰动引力梯度数据，建立其全张量表达式；对格网点值扰动引力梯度数据进行积分处理，获得格网均值扰动引力梯度数据，建立其全张量表达式；根据格网均值扰动引力梯度数据及其全张量表达式恢复地球重力场。本发明提供的确定方法及确定系统，能够充分挖掘各类观测数据包含的重力场信息，充分利用各类观测信息最大限度地恢复地球重力场，从而获取更加精确的地球重力场模型，为格网均值扰动引力梯度全张量数据联合反演提供了技术支撑。
Description
技术领域
本发明涉及物理大地测量学领域，特别是涉及一种卫星重力梯度全张量恢复地球重力场的确定方法及系统。
背景技术
当代物理大地测量学领域，构建重力场模型依然是研究热点。球谐函数在重力场模型构建过程中具有举足轻重的作用，地球重力场中各元素可以用球谐函数的级数形式来表达，以此建立起观测量与表征地球重力场信息的球谐系数(或椭球谐系数)的联系，通过重力场反演技术建立法方程，求解得到重力场模型。
卫星重力梯度测量数据可以用于重力场模型的构建，由于格网数据形成的法方程具有块对角结构，能够显著减少内存和时间消耗，特别对于海量的卫星重力梯度测量数据效果更明显。在构建重力场模型之前，首先，需要将沿轨扰动引力梯度数据转化为格网均值扰动引力梯度数据，然后，利用球谐函数积分式建立格网均值扰动引力梯度全张量表达式。目前，现有面球函数积分式只能用于建立部分扰动引力梯度分量的表达式，不能充分利用所有观测信息，无法获得更准确的地球重力场模型。因此，有必要研究重力梯度全张量联合反演地球重力场的方法。
发明内容
本发明的目的是提供一种卫星重力梯度全张量恢复地球重力场的确定方法及确定系统，能够充分利用所有观测信息，获得更准确的地球重力场模型，为格网均值扰动引力梯度全张量联合反演提供技术支撑。
为实现上述目的，本发明提供了如下方案：
一种卫星重力梯度全张量恢复地球重力场的确定方法，所述确定方法包括：
获取卫星重力梯度测量数据；
对所述卫星重力梯度测量数据进行预处理，获得沿轨扰动引力梯度数据；
对所述沿轨扰动引力梯度数据进行格网化处理，获得格网点值扰动引力梯度全张量；
对所述格网点值扰动引力梯度全张量进行积分处理，获得格网均值扰动引力梯度全张量；
根据所述格网均值扰动引力梯度全张量恢复地球重力场。
可选的，所述预处理具体包括：粗差剔除处理、坐标变换处理和正常场改正处理。
可选的，所述对所述沿轨扰动引力梯度数据进行格网化处理具体包括：
以格网计算点为中心，确定搜索半径范围内的沿轨扰动引力梯度数据；
利用扰动引力梯度水平变化率和扰动引力梯度垂直变化率将所述沿轨扰动引力梯度数据归算到格网计算点；
按各个所述计算点的贡献程度确定各个所述计算点的权值；
根据各个所述权值确定格网点值扰动引力梯度数据。
可选的，所述格网点值扰动引力梯度全张量的表达式为：
其中，Tpq(pq＝x′x′,y′y′,z′z′,x′y′,x′z′,y′z′)为格网点值扰动引力梯度数据的分量，x′、y′、z′为局部指北坐标系坐标，GM为地心引力常数，a为正常椭球长半轴，和为完全规格化地球引力位系数，n为阶数，m为次数，nmax为最大阶数，r、θ、λ为计算点的球坐标，r为计算点的地心向径，θ为计算点的余纬，λ为计算点的经度，为完全规格化缔合勒让德组合函数，x＝cosθ。
可选的，所述对所述格网点值扰动引力梯度全张量进行积分处理具体包括：
对所述格网点值重力梯度数据进行辛普森插值处理，获得插值多项式；
沿经度和纬度方向对所述插值多项式进行积分，获得积分多项式；
利用所述积分多项式对所述格网点值重力梯度数据进行积分，获得格网均值重力梯度数据。
可选的，所述格网均值扰动引力梯度全张量的表达式为：
其中，为格网均值扰动引力梯度数据的分量，x′、y′、z′为局部指北坐标系坐标，GM为地心引力常数，a为正常椭球长半轴，和为完全规格化地球引力位系数，n为阶数，m为次数，nmax为最大阶数，r、θ、λ为计算点的球坐标，r为地心向径，θ为余纬，λ为经度，βm为平滑因子，为完全规格化缔合勒让德组合函数的定积分，x1＝cos(θ-Δθ/2)，x2＝cos(θ+Δθ/2)，Δθ为格网数据余纬方向的分辨率。
可选的，所述根据所述格网均值扰动引力梯度全张量恢复地球重力场具体包括：
根据所述格网均值扰动引力梯度全张量构建地球重力全张量观测方程，所述地球重力全张量观测方程为：
AiX＝Li+εi，其中，(i＝x′x′,y′y′,z′z′,x′y′,x′z′,y′z′)，Ai为位系数矩阵，X为位系数参数向量，εi为观测误差向量，Li为格网均值扰动引力梯度数据的分量；
根据所述地球重力全张量观测方程构建全张量误差方程，所述全张量误差方程为：Aix＝li+Vi，其中，i＝x′x′,y′y′,z′z′,x′y′,x′z′,y′z′，x＝X-X0，li＝Li-AiX0，Ai为位系数矩阵，x为位系数的平差改正量，X0为位系数参数向量的估计值，li为格网均值扰动引力梯度数据的分量的残差向量，Vi为Li的平差改正数向量；根据所述全张量误差方程构建全张量法方程，所述全张量法方程：Nix＝Bi，i＝x′x′,y′y′,z′z′,x′y′,x′z′,y′z′，其中，Ni为法方程矩阵，Pi为格网均值扰动引力梯度数据的分量的权阵，Bi为法方程右向量；
根据所述全张量法方程恢复地球重力场，其中，位系数参数向量为：其中，i＝x′x′,y′y′,z′z′,x′y′,x′z′,y′z′，为法方程矩阵Ni的逆矩阵，σi为全张量法方程的权重。
一种卫星重力梯度全张量恢复地球重力场的确定系统，所述确定系统包括：
测量数据获取模块，用于获取卫星重力梯度测量数据；
预处理模块，用于对所述卫星重力梯度测量数据进行预处理，获得沿轨扰动引力梯度数据；
格网化处理模块，用于对所述沿轨扰动引力梯度数据进行格网化处理，获得格网点值扰动引力梯度全张量；
积分处理模块，用于对所述格网点值扰动引力梯度全张量进行积分处理，获得格网均值扰动引力梯度全张量；
地球重力场恢复模块，用于根据所述格网均值扰动引力梯度全张量恢复地球重力场。
根据本发明提供的具体实施例，本发明公开了以下技术效果：
本发明提供的一种卫星重力梯度全张量恢复地球重力场的确定方法及确定系统，根据格网均值扰动引力梯度数据及其全张量表达式进行重力场反演，能够充分挖掘各类观测数据包含的重力场信息，充分利用各类观测数据最大限度地恢复地球重力场，从而获取更精确的地球重力场模型，为格网均值扰动引力梯度全张量联合反演提供了技术支撑，具有重要的理论和实际应用价值。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图。
图1为本发明实施例1提供的卫星重力梯度全张量恢复地球重力场的确定方法的流程图；
图2为本发明实施例2提供的卫星重力梯度全张量恢复地球重力场的确定系统的结构框图；
图3为本发明实施例3提供的利用卫星重力梯度测量数据进行全张量地球重力场反演的流程图；
图4为本发明实施例3提供的建立格网均值扰动引力梯度全张量表达式的流程图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明的目的是提供一种卫星重力梯度全张量恢复地球重力场的确定方法及确定系统，能够充分利用所有观测信息，获得更准确的地球重力场模型，为格网均值扰动引力梯度全张量联合反演提供技术支撑。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
实施例1：
图1为本发明实施例1提供的卫星重力梯度全张量恢复地球重力场的确定方法的流程图。如图1所示，一种卫星重力梯度全张量恢复地球重力场的确定方法，所述确定方法包括：
步骤11：获取卫星重力梯度测量数据。
步骤12：对所述卫星重力梯度测量数据进行预处理，获得沿轨扰动引力梯度数据。所述预处理具体包括：粗差剔除处理、坐标变换处理和正常场改正处理。
步骤13：对所述沿轨扰动引力梯度数据进行格网化处理，获得格网点值扰动引力梯度全张量。
具体地，所述步骤13：对所述沿轨扰动引力梯度数据进行格网化处理具体包括：
以格网计算点为中心，确定搜索半径范围内的沿轨扰动引力梯度数据；
利用扰动引力梯度水平变化率和扰动引力梯度垂直变化率将所述沿轨扰动引力梯度数据归算到格网计算点；
按各个所述计算点的贡献程度确定各个所述计算点的权值；
根据各个所述权值确定格网点值扰动引力梯度数据。
最终获得的格网点值扰动引力梯度全张量的表达式为：
其中，Tpq(pq＝x′x′,y′y′,z′z′,x′y′,x′z′,y′z′)为格网点值扰动引力梯度数据的分量，x′、y′、z′为局部指北坐标系坐标，GM为地心引力常数，a为正常椭球长半轴，和为完全规格化地球引力位系数，n为阶数，m为次数，nmax为最大阶数，r、θ、λ为计算点的球坐标，r为计算点的地心向径，θ为计算点的余纬，λ为计算点的经度，为完全规格化缔合勒让德组合函数，x＝cosθ。
步骤14：对所述格网点值扰动引力梯度全张量进行积分处理，获得格网均值扰动引力梯度全张量。
具体地，所述步骤14：对所述格网点值扰动引力梯度全张量进行积分处理具体包括：
对所述格网点值重力梯度数据进行辛普森插值处理，获得插值多项式；
沿经度和纬度方向对所述插值多项式进行积分，获得积分多项式；
利用所述积分多项式对所述格网点值重力梯度数据进行积分，获得格网均值重力梯度数据。
最终确定的格网均值扰动引力梯度全张量的表达式为：
其中，为格网均值扰动引力梯度数据的分量，x′、y′、z′为局部指北坐标系坐标，GM为地心引力常数，a为正常椭球长半轴，和为完全规格化地球引力位系数，n为阶数，m为次数，nmax为最大阶数，r、θ、λ为计算点的球坐标，r为地心向径，θ为余纬，λ为经度，βm为平滑因子，为完全规格化缔合勒让德组合函数的定积分，x1＝cos(θ-Δθ/2)，x2＝cos(θ+Δθ/2)，Δθ为格网数据余纬方向的分辨率。
步骤15：根据所述格网均值扰动引力梯度全张量恢复地球重力场。
具体地，所述步骤15：根据所述格网均值扰动引力梯度全张量恢复地球重力场具体包括：
根据所述格网均值扰动引力梯度全张量构建地球重力全张量观测方程，所述地球重力全张量观测方程为：
AiX＝Li+εi，其中，(i＝x′x′,y′y′,z′z′,x′y′,x′z′,y′z′)，Ai为位系数矩阵，X为位系数参数向量，εi为观测误差向量，Li为格网均值扰动引力梯度数据的分量，即：
根据所述地球重力全张量观测方程构建全张量误差方程，所述全张量误差方程为：Aix＝li+Vi，其中，i＝x′x′,y′y′,z′z′,x′y′,x′z′,y′z′，x＝X-X0，li＝Li-AiX0，Ai为位系数矩阵，x为位系数的平差改正量，X0为位系数参数向量的估计值，li为格网均值扰动引力梯度数据的分量的残差向量，Vi为Li的平差改正数(残差)向量；
根据所述全张量误差方程构建全张量法方程，所述全张量法方程：Nix＝Bi，i＝x′x′,y′y′,z′z′,x′y′,x′z′,y′z′，其中，Ni为法方程矩阵，Pi为格网均值扰动引力梯度数据的分量的权阵，Bi为法方程右向量；
根据所述全张量法方程恢复地球重力场，其中，位系数参数向量为：其中，i＝x′x′,y′y′,z′z′,x′y′,x′z′,y′z′，为法方程矩阵Ni的逆矩阵，σi为全张量法方程的权重。
其中，Tpq(pq＝x′x′,y′y′,z′z′,x′y′,x′z′,y′z′)为格网点值扰动引力梯度数据的分量，x′、y′、z′为局部指北坐标系坐标，GM为地心引力常数，a为正常椭球长半轴，和为完全规格化地球引力位系数，n为阶数，m为次数，nmax为最大阶数，r、θ、λ为计算点的球坐标，r为地心向径，θ为余纬，λ为经度，为完全规格化缔合勒让德组合函数，x＝cosθ。
本发明提供的方法，首先对卫星重力梯度测量数据进行预处理，获得沿轨扰动引力梯度数据；然后对沿轨扰动引力梯度数据进行格网化处理，获得格网点值扰动引力梯度数据，建立其全张量表达式，并对对格网点值扰动引力梯度数据进行积分处理，获得格网均值扰动引力梯度数据，建立其全张量表达式；最后根据格网均值扰动引力梯度数据及其全张量表达式恢复地球重力场，能够充分挖掘各类观测数据包含的重力场信息，充分利用各类观测信息最大限度地恢复地球重力场，从而获取更加精确的地球重力场模型，为格网均值扰动引力梯度全张量数据联合反演提供了技术支撑，具有重要的理论和实际应用价值。
实施例2：
图2为本发明实施例2提供的卫星重力梯度全张量恢复地球重力场的确定系统的结构框图。如图2所示，一种卫星重力梯度全张量恢复地球重力场的确定系统，所述确定系统包括：
测量数据获取模块21，用于获取卫星重力梯度测量数据；
预处理模块22，用于对所述卫星重力梯度测量数据进行预处理，获得沿轨扰动引力梯度数据；
格网化处理模块23，用于对所述沿轨扰动引力梯度数据进行格网化处理，获得格网点值扰动引力梯度全张量；
积分处理模块24，用于对所述格网点值扰动引力梯度全张量进行积分处理，获得格网均值扰动引力梯度全张量；
地球重力场恢复模块25，用于根据所述格网均值扰动引力梯度全张量恢复地球重力场。
实施例3：
如图3所示，卫星重力梯度测量数据经过预处理和格网化处理后得到格网点值扰动引力梯度数据，对其进行积分处理，获得格网均值扰动引力梯度数据，重力场反演时需要建立格网均值扰动引力梯度全张量表达式，而两类新面球函数积分式的递推算法在这一过程中起到了关键性作用，通过该算法可实现格网均值扰动引力梯度全张量表达式的建立，结合观测数据进行重力场反演，最终得到更准确的地球重力场模型。
全张量地球重力场反演是采用所有扰动引力梯度数据(Tx′x′、Ty′y′、Tz′z′、Tx′y′、Tx′z′和Ty′z′)，以相应权重融合后进行全张量重力场反演，得到的模型包含的重力场信息更全面。采用本发明提供的确定方法及确定系统为格网均值扰动引力梯度全张量联合反演提供了技术支撑。
如图4所示，建立格网均值扰动引力梯度全张量表达式的流程为：首先，建立格网点值扰动引力梯度全张量表达式，其次，沿经、纬度方向积分，然后，引入完全规格化缔合勒让德组合函数Hnm的定积分Onm，最后，利用两类新面球函数积分式的递推算法计算Onm，从而建立格网均值扰动引力梯度全张量表达式，其具体实施步骤如下：
在局部指北坐标系中，将格网点值扰动引力梯度表示为：
其中，Tpq(pq＝x′x′,y′y′,z′z′,x′y′,x′z′,y′z′)为格网点值扰动引力梯度数据的分量，x′、y′、z′为局部指北坐标系坐标，GM为地心引力常数，a为正常椭球长半轴，和为完全规格化地球引力位系数，n为阶数，m为次数，nmax为最大阶数，r、θ、λ为计算点的球坐标，r为地心向径，θ为余纬，λ为经度，为完全规格化缔合勒让德组合函数，其表达式见公式(2)：
其中，为完全规格化缔合勒让德函数，x＝cosθ，y＝sinθ； 是与n、m相关的系数，表达式如下：
其中， 为克罗内克函数。
在重力场反演研究中，经常要用到扰动引力梯度平均值，如果已知格网中每一点的扰动引力梯度，则格网均值扰动引力梯度可以通过公式(3)所示的积分式求得：
式中，Tpq(pq＝x′x′,y′y′,z′z′,x′y′,x′z′,y′z′)为格网点值扰动引力梯度数据的分量，x′、y′、z′为局部指北坐标系坐标，Δσ是格网单元面积，σ为积分变量，Tpq为格网均值扰动引力梯度数据的分量，是格网内所有点值扰动引力梯度的积分中数。
因此，将位系数计算格网点值扰动引力梯度的表达式(1)代入式(3)，则得格网均值扰动引力梯度的表达式：
其中，为格网均值扰动引力梯度数据的分量，x′、y′、z′为局部指北坐标系坐标，GM为地心引力常数，a为正常椭球长半轴，和为完全规格化地球引力位系数，n为阶数，m为次数，nmax为最大阶数，r、θ、λ为计算点的球坐标，r为地心向径，θ为余纬，λ为经度，βm为平滑因子，x1＝cos(θ-Δθ/2)，x2＝cos(θ+Δθ/2)，Δθ为格网数据余纬方向的分辨率，为完全规格化缔合勒让德组合函数的定积分，定义如下：
其中，Δx＝x1-x2，x为积分变量。
因此，对公式(2)两边关于x同时进行积分，可得：
其中，和为第一类新面球函数积分式 为第二类新面球函数积分式 为完全规格化缔合勒让德函数，x＝cosθ，y＝sinθ，是与n和m相关的系数，其定义与表达式中的相同。
利用本发明获得的两类新面球函数积分式和可计算获得带入(4)式可计算得到格网均值扰动引力梯度全张量。
本发明获得的两类新面球函数积分式和 积分式中，当k＝0时，即为常规形式
完全规格化缔合勒让德函数与x、y的组合形式为 其积分形式为两类新面球函数积分式实质上与该式等价。在该积分式中，当k1为偶数时，可转换为当k1为奇数时，可转换为两类新面球函数积分式的递推算法足以解决重力场反演中的现有问题，可见，本发明提出的算法具有普遍性和全面性，该算法针对重力异常、扰动重力、扰动引力梯度等物理量的格网均值，给出无奇异点严密计算公式，在利用卫星重力梯度观测数据恢复地球重力场中，为格网均值重力梯度全张量联合反演提供了技术支撑。
本发明以缔合勒让德函数满足的基本递推公式为出发点，通过严密的理论推导，得到了两类新面球函数积分式和的递推关系式，x＝cosθ，y＝sinθ，x1＝cosθ1，x2＝cosθ2，k∈N，其推导过程如下：
首先，缔合勒让德函数满足下列关系式：
Pn,m+1(x)-Pn-2,m+1(x)＝(2n-1)yPn-1,m(x) (8)
(n-m)Pnm(x)+(n+m-1)Pn-2,m(x)＝(2n-1)xPn-1,m(x) (10)
其中，Pnm(x)为缔合勒让德函数，Pn(x)为勒让德函数，dmPn(x)/dxm为勒让德函数的m阶导数，n为地球引力位系数的阶数，m为地球引力位系数的次数，x＝cosθ，y＝sinθ。
当m＝n时，缔合勒让德函数的递推关系式满足：
Pnn(x)＝(2n-1)yPn-1,n-1(x) (11)
展开为关于y的形式：
Pnn(x)＝(2n-1)(2n-3)…3yn (12)
完全规格化缔合勒让德函数和缔合勒让德函数Pnm(x)的关系为：
当m＝n时，完全规格化缔合勒让德函数的递推关系式满足：
其中，
将(14)式展开为关于y的形式：
(1)第一类新面球函数积分式的递推算法
定义第一类新面球函数积分式为：
式(13)两边同乘以yk，然后带入公式(16)可得：
其中，
下面先推导的递推公式，然后根据(17)式可以得到的递推公式。定义积分式：
由于dy＝-(x/y)dx，故上式转化为：两边同乘以m+k+1，利用分部积分法进行积分，整理得：
式(9)两边同时乘以yk，然后关于x积分，则有：
用(20)式减去(19)式，可得：
将(8)式代入(10)式，整理得：
式(22)两边同时乘以yk，然后关于x积分，有：
将(21)式代入(23)式，可得：
当m＝n时，(24)式会产生奇异性，故对的递推关系式另外推导。式(12)两边同时乘以yk，然后关于x积分，可得：
由(17)、(24)和(25)式得m≠n时：
其中，
式(15)两边同时乘以yk，然后关于x积分，可得：
综上所述，可得：
其中，第一类新面球函数积分式递推公式的初值为：
其中，和的初值为：
(2)第二类新面球函数积分式的递推算法
定义第二类新面球函数积分式为：
式(13)两边同乘以xyk，然后带入公式(31)可得：
其中，
下面先推导的递推公式，然后根据(32)式可以得到的递推公式。定义积分式：
由于dy＝-(x/y)dx，故上式转化为：两边同乘以m+k+1，利用分部积分法进行积分，整理得：
式(9)两边同时乘以xyk，然后关于x积分，则有：
用(35)式减去(34)式，可得：
式(22)两边同时乘以xyk，然后关于x积分，有：
将(36)式代入(37)式，可得：
当m＝n时，(38)式会产生奇异性，故对的递推关系式另外推导。
式(12)两边同时乘以xyk，然后关于x积分，可得：
由(32)、(38)和(39)式得m≠n时：
其中，
式(15)两边同时乘以xyk，然后关于x积分，可得：
综上所述，可得：
其中，
第二类新面球函数积分式递推公式的初值为：
综上所述，本发明以缔合勒让德函数满足的基本递推公式为出发点，通过严密的理论推导，得到了两类新面球函数积分式和的递推算法：
本发明将两类新面球函数积分式的递推算法应用于卫星重力梯度测量中，解决了利用所有格网均值重力梯度全张量观测数据反演地球重力场的问题，充分利用各类观测值，最大限度地恢复地球重力场，实现准确反演地球重力场的效果。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。
Claims (5)
1.一种卫星重力梯度全张量恢复地球重力场的确定方法，其特征在于，所述确定方法包括：
获取卫星重力梯度测量数据；
对所述卫星重力梯度测量数据进行预处理，获得沿轨扰动引力梯度数据；
对所述沿轨扰动引力梯度数据进行格网化处理，获得格网点值扰动引力梯度全张量；具体包括：以格网计算点为中心，确定搜索半径范围内的沿轨扰动引力梯度数据；
利用扰动引力梯度水平变化率和扰动引力梯度垂直变化率将所述沿轨扰动引力梯度数据归算到格网计算点；
按各个所述计算点的贡献程度确定各个所述计算点的权值；
根据各个所述权值确定格网点值扰动引力梯度数据；其中，所述格网点值扰动引力梯度全张量的表达式为：
其中，Tpq为格网点值扰动引力梯度数据的分量，pq＝x′x′,y′y′,z′z′,x′y′,x′z′,y′z′，x′、y′、z′为局部指北坐标系坐标，GM为地心引力常数，a为正常椭球长半轴，和为完全规格化地球引力位系数，n为阶数，m为次数，nmax为最大阶数，r、θ、λ为计算点的球坐标，r为计算点的地心向径，θ为计算点的余纬，λ为计算点的经度，为完全规格化缔合勒让德组合函数，x＝cosθ；
对所述格网点值扰动引力梯度全张量进行积分处理，获得格网均值扰动引力梯度全张量；所述格网均值扰动引力梯度全张量的表达式为：
其中，为格网均值扰动引力梯度数据的分量，pq＝x′x′,y′y′,z′z′,x′y′,x′z′,y′z′，x′、y′、z′为局部指北坐标系坐标，GM为地心引力常数，a为正常椭球长半轴，和为完全规格化地球引力位系数，n为阶数，m为次数，nmax为最大阶数，r、θ、λ为计算点的球坐标，r为地心向径，θ为余纬，λ为经度，βm为平滑因子，为完全规格化缔合勒让德组合函数的定积分，x1＝cos(θ-Δθ/2)，x2＝cos(θ+Δθ/2)，Δθ为格网数据余纬方向的分辨率；
根据所述格网均值扰动引力梯度全张量恢复地球重力场。
2.根据权利要求1所述的确定方法，其特征在于，所述预处理具体包括：粗差剔除处理、坐标变换处理和正常场改正处理。
3.根据权利要求1所述的确定方法，其特征在于，所述对所述格网点值扰动引力梯度全张量进行积分处理具体包括：
对所述格网点值重力梯度数据进行辛普森插值处理，获得插值多项式；
沿经度和纬度方向对所述插值多项式进行积分，获得积分多项式；
利用所述积分多项式对所述格网点值重力梯度数据进行积分，获得格网均值重力梯度数据。
4.根据权利要求1所述的确定方法，其特征在于，所述根据所述格网均值扰动引力梯度全张量恢复地球重力场具体包括：
根据所述格网均值扰动引力梯度全张量构建地球重力全张量观测方程，所述地球重力全张量观测方程为：
AiX＝Li+εi，其中，i＝x′x′,y′y′,z′z′,x′y′,x′z′,y′z′，Ai为位系数矩阵，X为位系数参数向量，εi为观测误差向量，Li为格网均值扰动引力梯度数据的分量；
根据所述地球重力全张量观测方程构建全张量误差方程，所述全张量误差方程为：Aix＝li+Vi，其中，i＝x′x′,y′y′,z′z′,x′y′,x′z′,y′z′，x＝X-X0，li＝Li-AiX0，Ai为位系数矩阵，x为位系数的平差改正量，X0为位系数参数向量的估计值，li为格网均值扰动引力梯度数据的分量的残差向量，Vi为Li的平差改正数向量；
根据所述全张量误差方程构建全张量法方程，所述全张量法方程：Nix＝Bi，i＝x′x′,y′y′,z′z′,x′y′,x′z′,y′z′，其中，Ni为法方程矩阵，Pi为格网均值扰动引力梯度数据的分量的权阵，Bi为法方程右向量；
根据所述全张量法方程恢复地球重力场，其中，位系数参数向量为：其中，i＝x′x′,y′y′,z′z′,x′y′,x′z′,y′z′，为法方程矩阵Ni的逆矩阵，σi为全张量法方程的权重。
5.一种卫星重力梯度全张量恢复地球重力场的确定系统，其特征在于，所述确定系统包括：
测量数据获取模块，用于获取卫星重力梯度测量数据；
预处理模块，用于对所述卫星重力梯度测量数据进行预处理，获得沿轨扰动引力梯度数据；
格网化处理模块，用于对所述沿轨扰动引力梯度数据进行格网化处理，获得格网点值扰动引力梯度全张量；具体包括：以格网计算点为中心，确定搜索半径范围内的沿轨扰动引力梯度数据；
利用扰动引力梯度水平变化率和扰动引力梯度垂直变化率将所述沿轨扰动引力梯度数据归算到格网计算点；
按各个所述计算点的贡献程度确定各个所述计算点的权值；
根据各个所述权值确定格网点值扰动引力梯度数据；其中，所述格网点值扰动引力梯度全张量的表达式为：
其中，Tpq为格网点值扰动引力梯度数据的分量，pq＝x′x′,y′y′,z′z′,x′y′,x′z′,y′z′，x′、y′、z′为局部指北坐标系坐标，GM为地心引力常数，a为正常椭球长半轴，和为完全规格化地球引力位系数，n为阶数，m为次数，nmax为最大阶数，r、θ、λ为计算点的球坐标，r为计算点的地心向径，θ为计算点的余纬，λ为计算点的经度，为完全规格化缔合勒让德组合函数，x＝cosθ；
积分处理模块，用于对所述格网点值扰动引力梯度全张量进行积分处理，获得格网均值扰动引力梯度全张量；所述格网均值扰动引力梯度全张量的表达式为：
其中，为格网均值扰动引力梯度数据的分量，pq＝x′x′,y′y′,z′z′,x′y′,x′z′,y′z′，x′、y′、z′为局部指北坐标系坐标，GM为地心引力常数，a为正常椭球长半轴，和为完全规格化地球引力位系数，n为阶数，m为次数，nmax为最大阶数，r、θ、λ为计算点的球坐标，r为地心向径，θ为余纬，λ为经度，βm为平滑因子，为完全规格化缔合勒让德组合函数的定积分，x1＝cos(θ-Δθ/2)，x2＝cos(θ+Δθ/2)，Δθ为格网数据余纬方向的分辨率；
地球重力场恢复模块，用于根据所述格网均值扰动引力梯度全张量恢复地球重力场。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810695075.4A CN108873091B (zh) | 2018-06-29 | 2018-06-29 | 卫星重力梯度全张量恢复地球重力场的确定方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810695075.4A CN108873091B (zh) | 2018-06-29 | 2018-06-29 | 卫星重力梯度全张量恢复地球重力场的确定方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108873091A CN108873091A (zh) | 2018-11-23 |
| CN108873091B true CN108873091B (zh) | 2019-10-25 |
Family
ID=64296877
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810695075.4A Active CN108873091B (zh) | 2018-06-29 | 2018-06-29 | 卫星重力梯度全张量恢复地球重力场的确定方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108873091B (zh) |
Families Citing this family (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110532585B (zh) * | 2019-06-13 | 2023-06-06 | 中国测绘科学研究院 | 快速解算GOCE卫星重力场模型的Torus方法和系统 |
| CN112965124B (zh) * | 2021-02-08 | 2022-10-11 | 中国人民解放军92859部队 | 一种顾及局域保障条件计算外部重力异常垂直梯度的方法 |
| CN112965127B (zh) * | 2021-02-08 | 2022-05-31 | 中国人民解放军92859部队 | 一种基于重力异常计算外部扰动重力径向分量的方法 |
| CN113885101B (zh) * | 2021-09-28 | 2023-12-12 | 中国船舶重工集团公司第七0七研究所 | 一种基于椭球模型构建重力梯度基准图方法 |
| CN114089432B (zh) * | 2021-11-10 | 2023-04-18 | 中国地质大学（北京） | 一种利用卫星测高数据反演海洋重力梯度的频率域方法 |
| CN115327653B (zh) * | 2022-08-15 | 2023-04-04 | 自然资源部国土卫星遥感应用中心 | 一种基于张量不变理论的卫星引力梯度粗差探测方法 |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CA2488511C (en) * | 2002-06-28 | 2012-07-03 | Gedex Inc. | System and method for surveying underground density distributions |
| CN102567627A (zh) * | 2011-12-12 | 2012-07-11 | 中国人民解放军92859部队 | 基于卫星重力梯度观测数据的圆环面调和分析方法 |
| JP6468606B2 (ja) * | 2016-02-01 | 2019-02-13 | 滋樹 水谷 | 重力偏差データの表層密度値推定方法 |
| CN108020866B (zh) * | 2017-11-20 | 2019-11-12 | 中国空间技术研究院 | 一种星体重力场反演的方法和系统、以及处理器 |
-
2018
- 2018-06-29 CN CN201810695075.4A patent/CN108873091B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN108873091A (zh) | 2018-11-23 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108873091B (zh) | 卫星重力梯度全张量恢复地球重力场的确定方法及系统 | |
| Weiss et al. | High‐resolution surface velocities and strain for Anatolia from Sentinel‐1 InSAR and GNSS data | |
| Wang et al. | Present‐day crustal deformation of continental China derived from GPS and its tectonic implications | |
| CN110058236B (zh) | 一种面向三维地表形变估计的InSAR和GNSS定权方法 | |
| Hu et al. | Three-dimensional surface displacements from InSAR and GPS measurements with variance component estimation | |
| Wang et al. | Slip model of the 2015 Mw 7.8 Gorkha (Nepal) earthquake from inversions of ALOS‐2 and GPS data | |
| Abbondanza et al. | JTRF2014, the JPL Kalman filter and smoother realization of the International Terrestrial Reference System | |
| Teza et al. | Characterization of landslide ground surface kinematics from terrestrial laser scanning and strain field computation | |
| Li et al. | The contribution of the GRAV‐D airborne gravity to geoid determination in the Great Lakes region | |
| Cai et al. | A new algorithm for landslide dynamic monitoring with high temporal resolution by Kalman filter integration of multiplatform time-series InSAR processing | |
| CN108845366A (zh) | 卫星重力梯度张量对角线三分量反演地球重力场的调和分析法模型及其建模方法 | |
| Wang et al. | Crustal density structure, lithosphere flexure mechanism, and isostatic state throughout the Qinling Orogen revealed by in situ dense gravity observations | |
| Luo et al. | A new baseline linear combination algorithm for generating urban digital elevation models with multitemporal InSAR observations | |
| Wang et al. | Long-term continuously updated deformation time series from multisensor InSAR in Xi'an, China from 2007 to 2021 | |
| Li et al. | Inversion of GNSS vertical displacements for terrestrial water storage changes using Slepian basis functions | |
| CN115270608A (zh) | 基于arima与lstm的海岸带地面沉降预测方法 | |
| Feng et al. | A hierarchical network densification approach for reconstruction of historical ice velocity fields in East Antarctica | |
| Liu et al. | Dynamic estimation of multi-dimensional deformation time series from Insar based on Kalman filter and strain model | |
| CN103093101A (zh) | 基于重力梯度误差模型原理的卫星重力反演方法 | |
| Deng et al. | Transfer of height datum across seas using GPS leveling, gravimetric geoid and corrections based on a polynomial surface | |
| Bagherbandi et al. | Deflection of vertical effect on direct georeferencing in aerial mobile mapping systems: A case study in Sweden | |
| Chen et al. | Measuring aquifer specific yields with absolute gravimetry: Result in the Choushui River Alluvial Fan and Mingchu Basin, central Taiwan | |
| Xu et al. | Temporal and spatial movement characteristics of the Altyn Tagh fault inferred from 21 years of InSAR observations | |
| Tao et al. | The performance of LS and SVD methods for SBAS InSAR deformation model solutions | |
| CN114924270A (zh) | 基于GNSS的InSAR形变监测基准建立方法及装置 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |