【技术实现步骤摘要】
本专利技术涉及航空重力测量，具体是一种航空重力数据向下延拓的优化计算方法。
技术介绍
1、在大地测量方面，地球重力场用于确定地球形状和高程基准。重力测量是对连续重力场的离散采样，其中航空重力测量在地面重力测量难以达到的困难地区(如高山地区、极地)，以及在陆海交接重力空白区具有明显优势，可以快速获取中短波重力信息。航空重力数据是在飞行测线上的重力值，对于大多数大地测量应用来说，通常需要将空中的重力数据向下延拓到地面或其它参考面，用于多源重力数据融合或直接计算大地水准面。
2、航空重力向下延拓属于不适定反问题，延拓过程放大观测误差，将难以获得稳定可靠的延拓解。目前，解决这一问题的主要方法是求解逆poisson(泊松)积分方程。在poisson积分方程的基础上，进行变换后可以得到迭代法延拓模型，但迭代使观测噪声影响不断累积，被直接求逆的最小二乘延拓模型所替代。
3、为了改善延拓的稳定性，利用位系数模型的“移去-恢复”来改变重力数据的频谱特征。通过对poisson核函数进行修正，可以避免对长波的扭曲，减少远区效应的影响，改善延拓计算精度。常用的修正poisson核函数方法包含基于w-g思想的球面poisson核、基于molodensky(莫洛金斯基)方法的修正球面poisson核，其中球面poisson核的截断系数大于修正球面poisson核，从而修正球面poisson核更具优势。
4、航空重力测量数据受到各种误差的影响，当数据间隔较小或延拓高度较大时，法方程病态，最小二乘解算时需要进行正则化处理。选
5、综上，由于航空重力向下延拓观测方程病态，当数据间隔较小或延拓高度较大时，病态程度加剧，难以获得稳定可靠的延拓解。此外，poisson核截断系数计算通常采用数值积分法，当积分步长减少时，数值计算精度提高，但计算速度变慢。
技术实现思路
1、本专利技术的目的在于提供一种航空重力数据向下延拓的优化计算方法，建立poisson核截断系数严密计算方法，以及基于核修正、正则化、奇异值分解的航空重力向下延拓优化方法，以提高航空重力向下延拓数据的精度，并兼顾计算效率。
2、为实现上述目的，本专利技术提供如下技术方案：
3、航空重力数据向下延拓的优化计算方法包括以下步骤：
4、步骤一：根据poisson核截断系数严密公式计算poisson核截断系数；
5、步骤二：利用poisson核截断系数求解修正未知参数，对poisson核函数进行修正，计算远区影响，离散化近区后组成观测方程；
6、步骤三：对观测方程进行正则化处理，并基于奇异值分解计算延拓结果。
7、上述方法中，步骤一的实现方法为：
8、先使用公式(1)计算中间量rn,k(t)和rn,n(t)：
9、
10、式中，ψ0表示球冠积分半径，令t＝cosψ0表示球冠积分半径余弦，pn(t)表示n阶勒让德多项式，pk(t)表示k阶勒让德多项式，n表示待计算的poisson核截断系数阶数，k表示循环索引，取值范围从0到n；
11、再利用递推公式(4)计算中间量un(t,h)和vn(t,h)：
12、初值为公式(2)和公式(3)：
13、
14、递推公式(4)为：
15、
16、式中，h＝r/r，r为地球半径，r为计算点地心向径；
17、最后，根据公式(5)直接计算poisson核截断系数qn(r,ψ0,r)：
18、
19、所述步骤二的实现方法为：
20、首先计算修正poisson核函数：通过公式(6)可以形成l+1个线性方程组，求解方程组计算修正未知参数tm(r,ψ0,r)：
21、
22、式中，l表示修正阶数，n表示待计算的poisson核截断系数阶数，取值范围从0到l，m表示修正未知参数索引，取值范围从0到l；
23、修正poisson核函数km(r,ψ,r)表示为公式(7)
24、
25、式中，ψ表示计算点和流动点的球面角距；
26、然后，计算远区影响：
27、修正poisson核截断系数可表示为公式(8)：
28、
29、重力异常n阶拉普拉斯谐函数δgn(θ,λ)由位系数模型公式(9)计算：
30、
31、式中，gm为地球引力常数，pnm(cosθ)为完全规格化连带勒让德函数，和为完全规格化扰动位系数；
32、远区影响计算的级数形式为公式(10)：
33、
34、式中，nmax代表模型最大阶，θ表示计算点极距，λ表示计算点经度；
35、最后，组成观测方程：
36、系数矩阵的非对角线第i行第j列元素aij和对角线第i行第i列元素aii分别为公式(11)和公式(12)：
37、
38、其中，
39、
40、式中，dm(ri,ψ0,r)表示中间临时变量，nc表示积分半径内待求点个数；
41、观测方程为公式(14)：
42、
43、式中，δgi(r,θ,λ)表示航空重力观测值，δgj(r,θ′,λ′)表示待计算的延拓值，
44、表示远区影响，n表示待计算延拓值个数；
45、设有m个航空重力观测值，那么根据公式(14)，形成矩阵形式即公式(15)：
46、b＝ax (15)
47、式中，b表示公式(14)左端m个航空重力观测值减远区影响形成的列向量，a表示根据公式(11)、(12)、(13)计算形成的m行n列系数矩阵，x表示n个待计算延拓值组成的列向量；
48、所述步骤三的实现方法为：
49、首先，利用公式(16)对系数矩阵a进行奇异值分解svd：
50、a＝uλvt (16)
51、式中，u、v分别由左、右特征向量组成，u＝[u1 u2…um]，v＝[v1 v2…vn]，λ＝diag(λ1,λ2,…,λn)为对角矩阵，λi为奇异值；
52、然后，利用广义交叉验证gcv法确定正则化参数。设置正则化参数α选取范围和步长，进行遍历计算。航空重力数据向下延拓的正则化解为公式(17)：
53、
54、式中，表示正则化参数取α时解算的待计算延拓值，ui表示第i个左特征向量；
55、使用公式(18)计算中间量fα
56、
57、判定准则gα为公式(19)：
58、
<本文档来自技高网...【技术保护点】
1.一种航空重力数据向下延拓的优化计算方法，其特征在于，包括以下步骤：
2.根据权利要求1所述的航空重力数据向下延拓的优化计算方法，其特征在于，步骤一的实现方法为：
【技术特征摘要】
1.一种航空重力数据向下延拓的优化计算方法，其特征在于，包括以下步骤：
2....
【专利技术属性】
技术研发人员：赵辉，耿晓燕，王斌，王文超，孟静娟，陈俊英，
申请(专利权)人：自然资源部大地测量数据处理中心，
类型：发明
国别省市：
还没有人留言评论。发表了对其他浏览者有用的留言会获得科技券。