CN108415879B - 基于向上延拓的航空重力最小二乘向下延拓解析方法 - Google Patents

Abstract

本发明涉及一种基于向上延拓的航空重力最小二乘向下延拓解析方法，其主要技术特点包括：根据泰勒级数展开模型，建立航空重力向上延拓和向下延拓与测线高度上重力异常的各阶垂向偏导数的解析关系；通过Poisson积分式向上延拓航空重力数据，得到各高度上的重力异常值，基于最小二乘理论计算测线高度上重力异常的各阶垂向偏导数，将各阶垂向偏导数代入泰勒级数展开式实现向下延拓稳定解算。本发明将航空重力向下延拓解算过程转换为向上延拓计算和垂向偏导数解算两个步骤，通过第一步的处理有效抑制数据观测噪声对解算结果的干扰，通过第二步的处理成功实现向下延拓反问题的稳定解算，较好地解决了向下延拓解算固有的不适定性问题。

Description

基于向上延拓的航空重力最小二乘向下延拓解析方法

技术领域

本发明属于海洋重力技术领域，尤其是一种基于向上延拓的航空重力最小二乘向下延拓解析方法。

背景技术

航空重力位场向上和向下延拓是地球物理和大地测量数据分析应用中最重要的技术环节之一。航空重力位场向下延拓的目的是将数据观测面外推到更接近于位场源体，通过改善观测数据的信噪比，增强地球内部浅层质量分布异常的映射强度，来凸显地球重力场的局部变化特征，以提高位场数据解释推断的可靠性。但航空重力向下延拓在数学上属于典型的不适定反问题，其延拓算子对高频噪声具有明显的放大作用，很小的观测噪声也会引起延拓解算结果严重偏离真实的问题。

发明内容

本发明的目的在于克服现有技术的不足，提出一种基于向上延拓的航空重力最小二乘向下延拓解析方法，解决航空重力向下延拓的不适定问题。

本发明解决其技术问题是采取以下技术方案实现的：

一种基于向上延拓的航空重力最小二乘向下延拓解析方法，包括以下步骤：

步骤1、根据泰勒级数展开模型，建立航空重力向上延拓和向下延拓与测线高度上重力异常的各阶垂向偏导数的解析关系；

步骤2、通过Poisson积分式向上延拓航空重力数据，得到各高度上的重力异常值，基于最小二乘理论计算测线高度上重力异常的各阶垂向偏导数，将各阶垂向偏导数代入泰勒级数展开式实现向下延拓稳定解算；

所述步骤1的实现方法为：

设海拔高度为h_o的待求地球表面O点的重力异常为Δg_o，已知对应飞行高度面h_p上的空中重力异常为Δg_p，则根据泰勒级数展开模型得到Δg_o和Δg_p关系作为第一泰勒级数展开式，表示如下：

式中，Δh_po＝h_p-h_o代表空间P点相对于地面O点的高度差；

代表空中重力异常Δg_p在P点的n阶垂向偏导数；δΔg_po代表重力异常Δg_p到Δg_o的向下延拓改正数；

设位于飞行高度面上方、海拔高度为h_q处的空中重力异常为Δg_q，令Δh_qp＝h_q-h_p，则根据重力场解析延拓理论，将测线高度重力异常沿垂直方向的各阶偏导数与飞行高度面上方重力异常的解析关系作为第二泰勒级数展开式并表示为:

式中，δΔg_pq代表重力异常Δg_p到Δg_q的向上延拓改正数；

所述步骤2的实现方法为：

首先，在航空重力测量飞行高度面上方，以一定的间隔选取M个高度面Q₁，Q₂，…，Q_M，其相对于飞行高度面的高度差为：Δh₁，Δh₂，…，Δh_M；利用飞行高度面上的重力异常观测量Δg_p，依据如下向上延拓Poisson积分公式分别计算上述M个高度面上的重力异常

式中，r_p＝R+h_p，r_q＝R+h_q，R为地球椭球平均半径；Δg_pq代表与计算点Q相对应的飞行高度面上的重力异常；Δg_p为飞行高度面流动点的重力异常；

为计算点与流动点之间的空间距离，ψ为计算点与流动点之间的球面角距；

然后，将计算得到的M个高度面上的重力异常作为过渡观测量代入第二泰勒级数展开式，得到一系列以重力异常垂向偏导数作为未知数的观测方程；对于飞行高度面上的某个P点，使用处于不同高度面但在同一垂线方向上的M个重力异常

建立M个相对应的观测方程；设第二泰勒级数展开式的最高阶数N＝4，则观测误差方程为：

式中，v_i代表重力异常观测误差和向上延拓计算误差的综合影响；

进行如下设置：

其中i＝1,…,M；j＝1,…,N

A＝[a_ij]_M×N

V＝(v₁,…,v_M)^T

将所述观测误差方程表示为如下矩阵形式：

L+V＝AX

取M＞N，求得上式的最小二乘解为：

X＝(A^TA)^-1A^TL

将上式计算得到的各阶垂向偏导数代入第二泰勒级数展开式中实现向下延拓稳定解算。

本发明的优点和积极效果是：

1、本发明基于航空重力向下与向上延拓之间存在固有的内在联系，依据泰勒级数展开模型，将航空重力向下延拓解算过程转换为向上延拓计算和垂向偏导数解算两个步骤，通过第一步的处理有效抑制数据观测噪声对解算结果的干扰，通过第二步的处理成功实现向下延拓反问题的稳定解算，较好地解决了向下延拓解算固有的不适定性问题。

2、本发明利用超高阶位模型EGM2008建立的模拟标准场数据对最小二乘向下延拓解析模型解算结果的合理性和有效性进行了数值验证，试验证明本发明实用易行，具有较高的应用价值。

具体实施方式

以下通过实施例对本发明做进一步详述。

一种基于向上延拓的航空重力最小二乘向下延拓解析方法，依据泰勒级数展开模型，将航空重力向下延拓解算过程转换为向上延拓计算和垂向偏导数解算，包括以下步骤：

步骤1、依据泰勒级数展开模型，建立航空重力向上延拓和向下延拓与测线高度上重力异常的各阶垂向偏导数的解析关系。

在本步骤中，假设海拔高度为h_o的待求地球表面O点的重力异常为Δg_o，已知对应飞行高度面h_p(P点)上的空中重力异常为Δg_p，则由泰勒级数展开模型可知，Δg_o和Δg_p之间的关系记为第一泰勒级数展开式，表示为:

式中，Δh_po＝h_p-h_o代表空间P点相对于地面O点的高度差；

代表空中重力异常Δg_p在P点的n阶垂向偏导数；δΔg_po代表重力异常Δg_p到Δg_o的向下延拓改正数。由公式(1)可知，实现航空重力向下延拓计算的关键是，准确获取测线高度重力异常沿垂直方向的各阶偏导数。

假设位于飞行高度面上方、海拔高度为h_q(Q点)处的空中重力异常为Δg_q，令Δh_qp＝h_q-h_p，那么，根据重力场解析延拓理论，将Δg_q和Δg_p之间的关系记为记为第二泰勒级数展开式，表示为：

式中，δΔg_pq代表重力异常Δg_p到Δg_q的向上延拓改正数，其它符号意义同前。通过公式(2)可以得到测线高度重力异常沿垂直方向的各阶偏导数与飞行高度面上方重力异常的解析关系，即航空重力向上延拓和向下延拓与测线高度上重力异常的各阶垂向偏导数的解析关系。

步骤2、通过Poisson积分式向上延拓航空重力数据，得到各高度上的重力异常值，基于最小二乘理论计算测线高度上重力异常的各阶垂向偏导数，将各阶垂向偏导数代入泰勒级数展开式实现向下延拓稳定解算。

在本步骤中，在航空重力测量飞行高度面上方，以一定的间隔选取M个高度面(Q₁，Q₂，…，Q_M)，它们相对于飞行高度面的高度差为：Δh₁，Δh₂，…，Δh_M。然后利用飞行高度面上的重力异常观测量(Δg_p)，依据向上延拓Poisson积分公式分别计算上述M个高度面上的重力异常

Poisson积分公式为：

式中，r_p＝R+h_p，r_q＝R+h_q，R为地球椭球平均半径；Δg_pq代表与计算点Q相对应的飞行高度面上的重力异常；Δg_p为飞行高度面流动点的重力异常；

为计算点与流动点之间的空间距离，ψ为计算点与流动点之间的球面角距。

将计算得到的M个高度面上的重力异常作为过渡观测量代入泰勒级数展开式(2)，可得到一系列以重力异常垂向偏导数作为未知数的观测方程。对于飞行高度面上的某个P点，使用处于不同高度面但在同一垂线方向上的M个重力异常

可建立M个相对应的观测方程。考虑到航空重力测量分辨率的有限性和观测噪声干扰的现实性，这里取泰勒级数展开式(2)的最高阶数N＝4，由此可建立如下形式的观测误差方程：

式中，v_i代表重力异常观测误差和向上延拓计算误差的综合影响。令：

A＝[a_ij]_M×N

V＝(v₁,…,v_M)^T

可将式(4)表示为如下矩阵形式：

L+V＝AX (5)

取M＞N，可求得式(5)的最小二乘解为：

X＝(A^TA)^-1A^TL (6)

将式(6)计算得到各阶垂向偏导数代入泰勒级数展开式(1)实现向下延拓稳定解算。

下面采用超高阶位模型作为标准场开展数值计算检验及分析比较研究。以美国本土西部一个3°×3°区块(

λ:248°E～251°E)作为试验区，选用EGM2008位模型模拟产生不同高度面2′×2′网格重力异常及其垂向偏导数的“真值”。选择美国本土作为试验区，是考虑到EGM2008位模型在美国地区具有更高的逼近度。该区块属于地形变化比较剧烈的大山区，试验效果具有一定的代表性。由EGM2008位模型计算空中网格重力异常的公式为：

式中，GM为地球引力常数，R为地球椭球平均半径，r＝R+h，h为计算点高度。

为完全规格化缔合勒让德函数，

和

为完全规格化位系数。连续对式(7)求r的偏导数，可得到不同阶次的垂向偏导数计算模型。

这里取R＝6371km,h_i＝ikm(i＝0,1，…，10)，使用EGM2008模型(以360阶次作为参考场)分别计算11个对应于r_i＝R+h_i球面上的2′×2′网格重力异常“真值”Δg_i(tru)。表1列出了4个不同高度面上的EGM2008位模型(361～2160阶次)残差重力异常的统计结果。本发明采用的试验方案具体设计为：以5km高度面作为航空重力测量的观测面，即假设Δg₅(tru)是已知的观测量，依次采用不同的向下延拓模型，由Δg₅(tru)分别计算4km、3km、2km、1km和0km高度面上的重力异常

将计算值与相对应的“真值”Δg_i(tru)(i＝0,…，4)作比较，从而获得对应延拓模型解算精度的评估参数。

表1不同高度面位模型残差重力异常统计结果(单位：10^-5m﹒s^-2)

高度(km)	最小值	最大值	平均值	均方根
					0	-107.09	189.19	-0.44	34.28
1	-88.22	162.05	-0.39	30.26
					3	-68.34	120.33	-0.33	24.07
5	-55.22	90.66	-0.28	19.56

为了检验向上延拓Poisson积分方程(3)的计算效果，由Δg₅(tru)依据公式(8)分别计算6km、7km、8km、9km和10km高度面上的重力异常

求计算值与相对应“真值”Δg_i(tru)(i＝6，…，10)的互差，可得到向上延拓计算模型的精度评价，具体结果见表2。其中，积分半径统一取为ψ₀＝30′。为了减小积分边缘效应对评估结果的影响，计算区域边缘30′范围内的数据不参加对比分析，这个区域的数据也不再参加下一步的向下延拓计算检验。

表2向上延拓模型计算精度检核(单位：10^-5m˙s^-2)

高度(km)	最小值	最大值	平均值	均方根
					6	-0.27	0.28	-0.01	0.10
7	-0.21	0.23	-0.01	0.08
					8	-0.17	0.19	-0.01	0.07
9	-0.14	0.16	-0.01	0.05
					10	-0.12	0.13	-0.00	0.04

由表2结果看出，当观测数据不存在噪声干扰时，向上延拓模型计算结果可以达到很高的精度水平，即使是1km的延拓高度差，计算误差也不超过0.3mGal。

为了检验最小二乘向下延拓解析延拓模型的计算效果，首先在5km飞行高度面的上方，以0.5km的间隔选取20个计算高度面(即计算高度延伸到15km)，使用向上延拓积分式(3)，由Δg₅(tru)分别计算对应于上述20个高度面上的重力异常

以每一个2′×2′网格点作为计算单元，将前面计算得到的同一点但在不同面上的20个重力异常值代入式(6)，解算得到相对应的重力异常垂向偏导数，并将其代入式(1)，进一步依据设定的延拓高度差即可求得想要的向下延拓值。求计算值与真值的胡茬，可得到向下延拓计算式(6)的精度评价，具体结果见表3。

表3最小二乘解析模型误差均方根(单位：10^-5m˙s^-2)

高度差(km)	1	2	3	4	5
						N＝1	0.21	0.54	1.00	1.62	2.41
N＝2	0.07	0.16	0.30	0.50	0.80
						N＝3	0.06	0.14	0.26	0.42	0.65
N＝4	0.07	0.18	0.33	0.56	0.87

由表3结果看出，最小二乘解析延拓模型的计算精度随垂向偏导数最高阶次取值的增大而逐步提高，但当最高阶次增大到N＝4时，计算精度反而有所下降，这个结果一方面说明增加高阶项对提高解析延拓模型的计算精度具有重要作用，另一方面也说明增加高阶项越多，对数据观测质量的要求也越高。当延拓高度差小于3km时，该模型的计算精度优于1mGal。

为了进一步考察数据观测噪声对向下延拓计算结果的影响，这里在模拟观测量Δg₅(tru)中分别加入±1mGal、±3mGal和±5mGal的白噪声，对应生成三组带噪声的观测量，并重复前面的计算和比对检核过程。加入观测噪声后，与表2相对应的向上延拓模型计算精度检核结果(这里只列出互差均方根值(rms)，下同)如表4所示，与表3相对应的向下延拓模型计算精度检核结果如表5所示。

表4噪声影响下的向上延拓误差均方根值(单位：10^-5m˙s^-2)

高度(km)	6	7	8	9	10
						噪声＝1mGal	0.26	0.21	0.17	0.14	0.12
噪声＝3mGal	0.70	0.56	0.45	0.37	0.31
						噪声＝5mGal	1.18	0.95	0.78	0.65	0.54

表5噪声影响下的最小二乘向下延拓误差均方根值(单位：10^-5m˙s^-2)

对比表2和表4数值结果可以看出，数据观测噪声对向上延拓计算结果的影响很小，几乎不改变计算参数的变化形态，说明向上延拓模型作为低通滤波器的作用相当明显。进一步对比表3和表5结果可以看出，得益于向上延拓模型抑制高频噪声的优良特性，最小二乘向下延拓模型计算结果也几乎不受数据观测噪声的影响，在延拓高度差5km的高度内，仍可获得2.1mGal左右的计算精度，证明本发明提供的技术方案实用易行，具有较高的应用价值。

需要强调的是，本发明所述的实施例是说明性的，而不是限定性的，因此本发明包括并不限于具体实施方式中所述的实施例，凡是由本领域技术人员根据本发明的技术方案得出的其他实施方式，同样属于本发明保护的范围。

Claims (1)

1.一种基于向上延拓的航空重力最小二乘向下延拓解析方法，其特征在于：包括以下步骤：

步骤1、根据泰勒级数展开模型，建立航空重力向上延拓和向下延拓与测线高度上重力异常的各阶垂向偏导数的解析关系；

步骤2、通过Poisson积分式向上延拓航空重力数据，得到各高度上的重力异常值，基于最小二乘理论计算测线高度上重力异常的各阶垂向偏导数，将各阶垂向偏导数代入泰勒级数展开式实现向下延拓稳定解算；

所述步骤1的实现方法为：

设海拔高度为h_o的待求地球表面O点的重力异常为Δg_o，已知对应飞行高度面h_p上的空中重力异常为Δg_p，则根据泰勒级数展开模型得到Δg_o和Δg_p关系作为第一泰勒级数展开式，表示如下：

式中，Δh_po＝h_p-h_o代表空间P点相对于地面O点的高度差；

代表空中重力异常Δg_p在P点的n阶垂向偏导数；δΔg_po代表重力异常Δg_p到Δg_o的向下延拓改正数；

设位于飞行高度面上方、海拔高度为h_q处的空中重力异常为Δg_q，令Δh_qp＝h_q-h_p，则根据重力场解析延拓理论，将测线高度重力异常沿垂直方向的各阶偏导数与飞行高度面上方重力异常的解析关系作为第二泰勒级数展开式并表示为:

式中，δΔg_pq代表重力异常Δg_p到Δg_q的向上延拓改正数；

所述步骤2的实现方法为：

首先，在航空重力测量飞行高度面上方，以一定的间隔选取M个高度面Q₁，Q₂，…，Q_M，其相对于飞行高度面的高度差为：Δh₁，Δh₂，…，Δh_M；利用飞行高度面上的重力异常观测量Δg_p，依据如下向上延拓Poisson积分公式分别计算上述M个高度面上的重力异常

式中，r_p＝R+h_p，r_q＝R+h_q，R为地球椭球平均半径；Δg_pq代表与计算点Q相对应的飞行高度面上的重力异常；Δg_p为飞行高度面流动点的重力异常；

为计算点与流动点之间的空间距离，ψ为计算点与流动点之间的球面角距；

然后，将计算得到的M个高度面上的重力异常作为过渡观测量代入第二泰勒级数展开式，得到一系列以重力异常垂向偏导数作为未知数的观测方程；对于飞行高度面上的某个P点，使用处于不同高度面但在同一垂线方向上的M个重力异常

建立M个相对应的观测方程；设第二泰勒级数展开式的最高阶数N＝4，则观测误差方程为：

式中，v_i代表重力异常观测误差和向上延拓计算误差的综合影响；

进行如下设置：

其中i＝1,…,M；j＝1,…,N

A＝[a_ij]_M×N

V＝(v₁,…,v_M)^T

将所述观测误差方程表示为如下矩阵形式：

L+V＝AX

取M＞N，求得上式的最小二乘解为：

X＝(A^TA)^-1A^TL

将上式计算得到的各阶垂向偏导数代入第二泰勒级数展开式中实现向下延拓稳定解算。