最小二乘配置在物理大地测量中的应用

武汉大学

硕士学位论文

最小二乘配置在物理大地测量中的应用

姓名:邹贤才

申请学位级别:硕士

专业:固体地球物理学

指导教师:李建成

20030504

摘要

本文描述了最小二乘配置法的基本原理及其基础理论,详细探讨了经验协方差函数的确定方法。在此基础上,以我国西部高山地区作为实验分析范例,利用实测重力数据、数值地面模型和EGM96地球重力场模型,确定了该地区的大地水准面。在研究如何应用最小二乘配置时,本文在中国和美国两个区域用最小二乘配置方法联合重力数据,GPS水准数据解算大地水准面模型,并与多项式拟合法作了比较,初步探讨了最小二乘配置法在这一领域中应用的实用性。针对最小二乘配置耗费大量计算机时的问题,本文同时讨论了快速最小二乘配置,由于信号与噪声无关的格网数据,其协方差阵是块Toeplitz矩阵,这样,可用快速算法提高计算速度且不损失计算精度,并大大减少了计算机内存资源的需求。本文以此研究了测高数据的处理方法与一般理论,并由Geosat/GM测高数据计算了宽阔海域海平面高,用快速最小二乘配置方法获得了该地区的格网重力异常,讨论了阶方差模型对配置结果的影响。

【关键词】.最小二乘配置;大地水准面:协方差函数

Abs仃act

TheelementarytheoryofLeast。SquaresCollocation(LSC)inPhysicalGeodesyisdescribedindetailinthisthesis,aswellasthemethoddefiningthecovariancefunction,whichisthekeyproblemwhenutilizingLSCinphysicalgeodesy.Basedonthesepremises,anexperimentisimplementedbycombingtheobservedgravitydata.D下MandEGM96GlobalGeopotentialModelt0constructthelocalgriddedgeoidmodel。ThentheapproachofcombinggravityanomaliesandGPS/LeVe¨ingdataisdiscussed.ThecomputationsarecarriedoutbothinChinaandinUnitedStatesforascertainingtheconclusions.Whencomparedwiththemethodmorecommonlyusednow—polynomialfitting,thefeasibilityandvalidityofLSCmethodarepresented.Inordertoovercomethemostpnmarynumer]calproblemthatLSCisaCPUtime?consumingmethod,fastcollocationisdiscussedsubsequentlyduetothespecialstructureofcovariancematrixthatiftheinputobservationsaregridded、Sigal-Noiseirrelevant,thenthecovariancematrixoftheobservationsisablock/7blockToeplitzmatrix,thespenttimeandmemorycanbereducedsharplywithouttheexplicitlossofaccuracy.Thegeneraltheo吖ofprocessingaltimetrydataisthendiscussed.AfterprocessingtheGEOSAT/GMaltimetrydata,theSeaSurfaceHeightsintheopenseaarederivedandthenusedtopredictgfiddedgravityanomalieswiththefastcollocation;theeffectsofdegreevariancemodeIonLSCresultsarediscussedparticularly.

【KeyWords]:Least。Squarescollocation;geoid;covariancefunction

?U?

最小二乘配置在物理太地铡置中的应甩

第一章绪论

§1.1引言

大地测量学是描述地球形状的-I'q基础科学,其描述方法可以从几何和物理两个方面获取点位的空间信息、物理信息和时间信息。这两种方法菩有侧重,但有着紧密的联系。物理大地测量学研究的重点是地球形状的物理属性一地球形状及(外部)重力场,它在地学研究,特别是在大地测量学的研究中.起着十分基础的作用,是构成现代人地测量学的重要支柱。地球作为人类居住和生活的场所,它的存在与运动形式时刻制约、影响着人类活动。地球弓l力场,特别是外部引力场.作为对地球物质分布与运动的一种直接响应.更是成为了地球的物理属性中最重耍的一个方面,受到了广泛的关注:物理大地测量学的研究,成为地学领域中的~项基础性研究.并日蓝加强其在学术研究中的地位。

牛顿万有引力公式的发现,使得人们第一次能够从严格的物理意义上,用数学方法对地球重力场进行定性,定量的细致研究。但在很长的一段时间里.这一研究仅仅局限于理论层面,其贡献、成就主要体现在理论上。上世纪60年代.第颗人造地球卫星升空.人类跨入了空闻时代,重力场研究的紧迫性日益凸显,其实用价值,越来越受到人们的关注。高空飞行器飞行弹道设计、空间飞行器的精密定轨,空间受力分析.都急切需要姗道地球;I力场的知识和地球外部重力场的空间分布细节。这些使得物理大地测量成为丁一门应用性很强的学科。同时,空间科学以及以计算机为代表的信息技术的发展,也极大地促进了物理大地测量学的研究:物理大地测量的进步,不仅体现在实际的成某(如:高阶敬地球位模型)上:方法上,也在这些外因的影响下,发生了深刻的变革,产生了新的理论问题:重力一测高边值问题,混合边值问题,超定边值问题,离散型边值问题等等。这门刨立于19世纪的学科,也得到了全面的发展。具体来讲.随着多颗带有地学任务的卫星实验.人们获得了大量的对在轨人造天体的观测资料,由于这些观测信息包含着丰富的地球弓l力场信息.对它们的处理使人们在地球重力场研究中获得了具有跨时代意义韵进步,Kaula于1966年综台利用卫星轨道摄动资料和地面重力测量资料,导出了一个8阶地球位模型。从今天的观点来看.这一成果所采用的方法实际上更具有意义。Kaula(19661将其成果写入《卫星大地测量学》一书,成为卫星重力学和卫星大地测量学的代表性著作。晟新一代的全球重力场模型EGM-96。己达到.360阶汝,分辨率50kra,而甚高阶,超高阶局部

第一章鳍论

将为研究活动构造带地壳运动提供重要的地球内部构造信息。在经典大地测量中,大地水准面只作为大地测量观测值的归算和研究地球的外部形状的依据,徐菊生、赖锡安等1997年研究了华北地区大地水准面起伏及其在地质构造上的意义后,指出华北地区大地水准面起伏的总体形态和该地区的地壳厚度呈镜像关系。这是现代大地测量对地学研究的一项新贡献。物理大地测量学的结果与地震学结果相互验证.必将摧动我们对地球魂部构造理解的深化,随着重力测量技术以及多种观测手段的使用,我们将获得更为丰富的重力资料,重力场模型的精度与分辨率将会得到不断的提高,特别是在局部重力场方面,这一进度将会更快。李建成等处理了大量的卫星测高资料、水准资料和GPS观测数据后,得到了多个高精度、高分辨率的局部似大地水准面。在我国许多地区,在较大区域精度可以达10aIl以内,在小区域(城市)可达2~3aIl。因此,大地水准面将不再仅仅局限于作为归算面,而会在地学研究中起到越来越大的作用,也可以预见,会有新的应用方向产生。

地球重力场模型的重要性,将使物理大地测量学的发展得到更大的推动。传统的获取重力资料的方法是进行地面重力测量和天文观测。这一方法仍然是现在获取重力资料的重要手段,但其精度已经得到了巨大的提高。目前.绝对重力测量的精度已经达到了1HGal,比如新一代的FG5绝对重力仪。在电子技术的支持下,uGaI级的便携式相对重力仪已经铍,。泛地应用于大范围重力测量。重力探测技术的另一重要进展是卫星重力技术的应用,新的卫星重力探测计划.如德国的CHAMP计划,NASA的GRACE计划和欧空局(ESA)的GOCE计划.目的都在于以前所未育的精度和分辨率确定重力场,这些计划如果能够实现,那么在重力场模型研究中将获得重大突破,预计全球大地水准面精度可以达到1~2c【【l,分辨率达到80h。

当然,地球重力场模型的研制不能仅仅依靠数据,它还必须有坚实的理论支持。因此,大地测量边值问题和重力场逼近技术是研究的基础。

§1.2物理大地测量边值问题及重力场逼近方法

物理大地测量边值问题(GBVP)是地球重力场模型研究的基础,也是我们进行实际计算的理论根据a这一边值问题可以表述为:在大地水准面或自然表面给定边值条件及其边值,确定该边界面及其外部引力位,使其满足边值条件。物理大地测量边值问题的发展过程可以大致分为Stokes边值问题。M010densky边值问题和超定边值问题三个阶段。

Stokes定理理论上的意义非常重大,同时也能以很高的精度满足研究地球形状的要求。但造.它本身有缺陷一大地水准面外部无质量这一条件无法实现。虽然入们讨论了各种重力归算方法,提出了许多模型.但是,这些方法都要改变地球质量的真实分布。这样一来就改变了地球外部重力场。Molodensky于1945年以地球的自然表面为边界面.提出了新的边值问题--Molodcnsky边值问题。这一理论克服了Stokes问题要求大地水准面外无质量的困难,并不需要重力归算,无须对地球岩石圈的密度作假设,(注:实际计算中,仍然有重力赔算这一步骤,但是,在Molodensky问题中的重力归算与Stokes问题中的蓬力归掉目的、意义已经不同了),解算出来的是地球的真正形状。这样,大地测量边值问题的研究重点很快就转到了研究Mdodensky问题这一方向上来。但Molodensky问题却是一个

-2?

最小二乘配置在物理大地测盘中的应用

异常复杂的数学问题一它的边界面是地球自然表面,Molodensky问题是一个高度非线性化自由边值问题,,解算非常困难。Bjerhammer于1964年提出了一种新方法,他以埋藏在地球某一深度的虚拟球面为边界面,将地面重力值延拓到这一边界面上,使由新边值导出的外部位不变。这一方法,解是以Stokes公式的形式给出,而虚拟球面上的重力异常由地面重力异常向下延拓而得.数学上是一个积分公式,我国学者许厚泽、朱灼文将这一方法发展为虚拟单层密度法.并证明了两种方法的等价性。但是.向下延拓在数学上是非适定的;因此,在应用上.这种理论也受到了限制。

随着观测手段的不断更新,人们获得了各种不同的边值:在广阔的海洋上,通过卫星测高技术,可以直接得到大地水准面的观测值,通过梯度测量,可以得到梯度边值;这样,在边值问题中.就获得了比数学上,严格求解所必需的边值更多的边值条件,而这些信息都是非常有用的。如何综合运用这些观测信息?这些信息都是有观测误差的,尽管它{f3的场源相同,但实际上这些观测值却不相容。为解决这些问题,人们提出了超定边值问题。

对大地测量边值问题的研究,有助于我们对重力场的了解.为地球外部重力场的逼近提供理论基础,并保证解的稳定性,控制解的精度。在这些理论的支持下,重力场模型研制进展迅速。由于Stokes问题本身的缺陷以及精度上的原因,当今的重力场模型主要建立在Molodensky理论的基础上。早期由于地面重力赍料有限,主要以对卫星轨道的观测,通过轨道摄动分析进行计算:由于卫星轨道只对重力场的长波信息敏感。因而。这一方法只能获得重力场的低阶位系数,研制的单位有美国的史密松天文台,哥达德宇航中心,直到上世纪60年代,其最高阶也不过50阶。高阶地球重力场模型的研制开始于上世纪70年代,美国的俄亥俄州立大学(OSU)率先发表了第一个180阶次的重力场模型(Rapp78),并形成了自己的系列,其最新的360阶次重力场模型为OSU91A1F。德国汉诺威大学发表了GPM系列模型,美国德克萨斯大学的重力场模型系列名为TEG,我国也于1977年首次提出了DQM77A/B(20阶和22阶),而后又有武汉铡绘科技大学推出的WDM系列,最新的模型是由李建成等推出的wDM2000重力场模型。

这些地球重力场模型的研制,大体上采用的思路是:联合由卫星轨道观测值分析得出的低阶位系数与地面重力测量值,获得改进的低阶位系数及地面重力值,然后,根据重力异常与地球外部扰动位位系数的积分关系.计算中高阶位系数;与低阶位系数叠加后形成了重力场模型的位系数集。位系数所能展开的最高项与数据的分辨率有关,可由Nquist频率公式确定。因此,计算上,高精度、高分辨率的模型要有高分辨宰的实测数据,但这一条件不是充分的;理论上,把地球作为一个球体,对于Stokes公式和Molodensky公式都会带来极大的简化,当需要以分米甚至厘米级的水平确定大地水准面时,则需要考虑地球的扁率。OSU86重力场模型以后的研制,已经开始采用Molodensky理论,运用K.rarup给出的严密线性化边值条件,进行解算。Pavils(Rapp,1988)推导了考虑扁率级的边值条件,

a7’月71a"o。.

给出了三项椭球改正,即矗(导方向化为=e二-L方向的改正),£,(詈方向化为.方向.v_L

的改正)和£。(将重力异常投影到椭球法线方向的改正)的计算公式,讨论了将正常重

?3?

第一章绪论

力值展开到关于高程的二次项,将地面重力值延拓到参考椭球面,顾及大气影响并考虑了不同区域高程基准的差别等这些在全球重力场模型解算中遇到的主要问题。

§1.3本文的研究方法和主要内容

1.3.1最小二乘配置方法

最小二乘配置方法是上世纪60年代以来大地测量学界在局部重力场逼近理论方面所取得的最重要的一项进展。在此以前的理论。共同点之一是它们所建立的已知量到未知量的泛函关系都是从一类重力场量到另一类重力场量的映射.它们的关系是严格确定的。而我们早就认识到,确定性和随机性,必然性和偶然性并没有严格的区分。晟小=乘配理论的产生,其重要意义就在于第一次在重力场理论中引入了随机性概念。虽然它有很多不完善的地方.但其优点也是非常明显的。最小二乘配置方法可以联合使用多种观测值,统一解算,这也是Krarup1973年提出整体大地测量概念的基础。随着观测技术的不断进步和观测手段的多样化,我们可以获得多种高精度的观测值,这些观测值所包含的丰富信息则是进一步精化重力场模型所必须的.这样就产生了如何综合运用多种观测数据的问题。从静止的观点上看(对于重力场的变化,就现在的目的而言已足够精确),地球重力场是唯一的:因而,以其为场源的重力边值应该是相容的;但由于观测误差,边值理论指出,超定边值问题是无解的。为了运用多种数据,现行的方法中有采取转换的方法。如:运用逆Stokes公式将由卫星测高推求的大地水准面高转换成重力异常,加入到重力异常数据中参与重力场模型计算,那么这一转换过程是否会造成某些重力场信息的遗失呢?从这个角度上讲。需要有一种能够综合运用多种观测数据的方法。最d'---乘配置方法采用了统计的观点来分析数据,从而具备这种特性,在重力场模型计算中显示了它独具的优点。运用这种方法处理数据,其关键是需要拟台出一个能充分表达研究区域范围内重力场特性的协方差函数。其优良特性是若某一类数据可以由对扰动位施以算予而得,则其协方差函数可以通过对扰动位协方差函数施以相应的算子而得,对互协方差也是如此。这样,如果获得了扰动位的协方差函数模型.那么对于重力场观测量.由于它们均与扰动位存在泛函关系,因此就可以获得任意观测量之间的协方差,从而对多种观测数据一并加以解算。同时,我们还可以对所获得的结果进行精度评定。从测量的角度来看,这也是一个非常好的特性。Tscheraiag/Rapp于1974年提出了一种获取协方差函数的方法,印认为可以用函数模型来拟合地球位模型的阶方差,然后,拟合出这个函数模型的参数:具体来讲,就是用函数模型来表达重力场的整体变化特性、一般变化规律,而用参数来表达实际观测数据的特性。他们在合作中是以全球模型为研究对象的。对于局部地区,则可以把在一个局部范围内为常量的长波部分的影响剔除,而用中高阶部分模型来逼近实际数据计算出的协方差,即将低阶部分截断。Knudsen于1987年进一步发展了这种方法.他引入了先验位模型的误差阶方差来表达位模型在某一地区的拟合程度,高阶部分则仍采用函数拟合的方法。这些方法,其主要工作就是要得到一个能表达(局部)重力场特性的协方差函数,其计算则以由观测数据计算出的协方差为基础。现今.国际上采用最小二乘配置方法大体上都遵循这种思路。

‘4’

虽d'--乘配置在物理大地测量中的应用

本文的工作也以此为基础展开。

在已开展的最小二乘配置研究中,主要以Tscheming和Rapp教授的工作具有代表性。由Tscheming教授提出的阶方差函数模型可以导出闭台的协方差函数表达式;这样,在大规模计算中运用最,J,--乘配置就变得可行。这一内容在文献[45】中有详细的阐述;这之后的很多工作,都是以此为基础或稍作改进。Rapp教授则是多次在测高问题中运用最小二乘配置方法:主要以1983年,他利用全球Seasat卫星测高数据,用最小二乘配置在10个地区计算了重力异常:获得的1。×r重力异常标准差为±5.1magl,5。×5。重力异常标准差为±2.7magl。在欧洲许多国家,也采用过最小二乘配置方法计算局部大地水准面。但这一方法没有大规模应用开来。主要原因是协方差函数的确定是最小二乘配置中的一个难点:另外,做最小二乘配置需要解算一个维数与观测数据个数相等的方阵的逆。当观测数据很多时,计算工作量巨大,因此,对计算机的软硬件环境提出了很高的要求。这使得最d,22乘配置主要应用在局部重力场逼近中,即便这样,最小二乘配置计算也是非常耗费计算机时的。

1.3.2本文的研究内容

本文的主要研究目的是通过大规模的数值分析,掌握晟d'-乘配置在物理大地测量中应用的几个方面。主要结构如下:

1.系统阐述晟小二乘配置理论及其在局部重力场逼近中的应用:

2.探讨快速最小二乘配置方法并实现算法:

3运用最小二乘配置方法计算局部大地水准面,通过与FFT方法的比较来讨论这一方法的效果:

4.运用最d'--乘配置方法做GPS水准和重力异常的联合配置,并与多项式拟合法进行比较:

5.推估海洋重力异常并讨论阶方差模型的选择对海洋重力异常估算的影响。

?5。

第二章最小二乘配置理论

第二章最4"---乘配置理论

§2.1最小二乘预估及协方差函数在重力场中的初步应用

用解析逼近值去拟合一定数量的线性泛函以确定一个函数的方法称为配置法:当选取最4"--乘准则时Ⅲq称这种方法为最小二乘配置方1法。在重力场逼近中,最小二乘配置方法一个很大的优点就是可以同时利用多种重力观测值来确定地球外部重力场,运用这一方法的关键是要获得~个符合实际情况的协方差函数。下面,以最小二乘预估为例.指出协方差函数在重力场中的应用。

设有待估信号s”,观测信号为11,s和I都是列向量,且期望值为O:即:

E{l}=0,E{s}=0(2.1、这两点假设是符合重力场中残差信号的一般情况的。

重力场中的各类元素,由于场源相同.它们之间存在著相关性,因此s”=[s。.52,,….,s。]1和14=[,I,,2,.…,f4】1之间存在函数关系,可记为:s“=s“U9J,

从简化问题的角度可以假设s“和11之间用线性关系表示已具有相当的精度,则s“=珊-,H是线性算子。这是一种逼近,给予不同的评价标准,H会有不同的表达形式,或者说,

最优逼近是在某~前提下的最优:设观测值11的自协方差C“以及s”和11的互协方差C。都己知,且l9的协方差阵C。正定,在待估信号的误差方差为虽小的准则下,则可以求出H

的表达式。因此可设:

§=Hl(2.2)且CⅣ=cov(1,Z),C,,=COV(S,,)已知;CⅣ正定。令s=;一s,这里:§是s”关于19的线性最小方差估值,占为信号估计误差。它的协方差阵为:

C。=E{船7)=E{G—s)(i—s)7)(2-3)这是一个方阵,其对角线上的值为信号估值5,(f=1,2,.…,m)的误差方差。按照方差最小原则,经过不是很复杂的推导,可知H=C“CB。【奠里兹,10841,带入(2-2)式.可得矿的线性最小方差为:

;=C“CⅣ“l(2.4)如果观测值14含有带偶然误差特性的观测噪声,则问题的解变为:

‘6‘

最小二乘配置在物理大地测量中的应用

§=C,,(CⅡ+D)。I(2-5)△g,:【c,。.c,:,c,,f暑誊!j暑;1f攀:1。:.。,

△g,如。.c%c。]Ic::G:…c2w卜l弘。,

lC。,C。:…C.Njl△g。j

上式中,c,,(i=1,_,2…N)是待估信号与观测值之间的互协方差,G(id21…2….N)是观测值之间的协方差,分别与C。,和CⅣ的意义对应。在重力场计算中,会经常碰到内插与外推

等情况.如果在计算区域已知某些观测量之间的协方差函数C(d),且C(d)只与两点间的位置有关,那么在已知了某一区域的协方差函数后就可以很方便地计算不同量间的协方差阵,运用公式(2-4)或(2.5)就可以很方便地进行预估。

§2.2空间扰动位协方差函数

重力异常的协方差函数对重力异常的内插与推估是必不可少的,在最4,-乘配置中,扰动位的协方差函数也是一个最基本的条件。设r(P)和,(Q)是球面上两点的扰动位f,在这里引入球面上的“期望”算子M{?},M算子求平均的意义,是指在整个球面上和所有的方位上取得平均值【海斯卡涅,莫里兹,1967]。所以,我们假设P,Q是在一个充分近似参考椭球的球面上,它的半径为R,扰动位r在这个球r=胄外调和,并称算子M{?}是齐性的(在整个球面上求平均值)和各向同性的(在所有的方位上求平均值)a

定义:

扰动位协方差函数K(P,Q)=M{r(|P),丁(Q))

由于算子M的各向同性,则K(P’Q)必定是只与P,Q两点间的距离矿有关。由上面的定义可以写出:

K(P,Q)=足(P)=嘉££。去£r(只A)r(p‘msin甜&删a(2.7)

式中,臼。A是球面坐标,a表示方位角;内层积分表示了各向同性。而外层积分则表示了齐性。由球面三角学知,v满足下式:

?'?

嚣二章最小二乘配置理论

COS∥=cos口cos口+sint9sincos(2一彳)(2?8)

如果令参考椭球的质量等于地球的总质量的球函数展开式中将不含有零阶和一阶项参考椭球质心与地球质心重合,那么,T(O,五)则地球外部扰动位的球谐展开公式:

m,¨,=孚薹㈣mtg[(a"。costa2+瓦sinmA.)/只胁)】}陆。,可得函数K够)的表达式可以写为【莫里兹,1984]:

足(∥)=∑t只(cos∥)

^。2

其中:只(cosg")是勒让德多项式

t。:窆毓+砝),瓦。与’L是完全规格化的外部扰动位位系数。

m=O

(2.10)

上述推导扰动位协方差函数的方法与文献【海斯卡涅,莫里兹,19671推导重力异常协方差函数的过程是一样的,其结果也很相似。

如果要求协方差函数K(P,Q)在球外调和,由于r在球,=R外是调和的,掇据算子M的定义式,K(P,Q)无论作为P的函数还是作为Q的函数。都应该是调和函数。而n阶球函数在球外调和时与向径r的联系关系为,一‘””,因此,可以写出足(P,Q)的形式

础珈砉焉筹洳帕∞,-,,,r分别表示P,Q点的向径,当r=,=R时,芷(尸,Q)就定义在球面上,可得

COnSI=k。;所以,扰动位的空间协方差函数为

靴)=纠玎咖sy)(2.12)

最小二乘配置在物理大地测量中的应用

§2.3协方差传播定律

我们在上一节中导出了外部空间扰动位的协方差函数,但是,并不是每种观测量的协方差函数都必须从定义出发推导:因此,本节将主要讲述协方差函数传播定律。

事实上,重力场观测量,如重力异常△g,垂线偏差《,刁),大地水准面高等都与扰动位,存在(或近似)着线性泛函关系。

△g:一坚一三丁

or,

,1OT

C=~

rgdD

lar叩;一——i

rycos∞叭(2?13)(2?14)(2-15)

因此,对于这些观测值?可以统一地记为:f.=L,T;这里,L,表示联系观测值f,与扰动位r的算子。如果记

B=(2-16)

}c尸,=p一,c,∥’’,c。,巨。暑|j:妻1院1cz?,,,

?9?

第二章最小二乘配暨理论

上式中G为信号的协方差,但我们只有扰动位的协方差函数,因此,要根据扰动位的协方差函数推导其泛函值的协方差;

记f。=Lfr(Q),LiQ表示泛函L。作为Q的函数作用于丁。则有:

C。=吖{r(P),馨r(Q))=砰M{r(P),r(Q))=馨K(P,Q)(2—18)q=M{Le,T(P),碍r(Q)}=£?碍M{r(P),r(Q))=茸碍K(P,Q)(2一19)

其中,£的下标i,J表明算子的类型,L的上标P,Q表明算子作为独立点P或Q的函数作用于丁的协方差函数。上式即为协方差传播律:即知道了扰动位的协方差函数,则对任何重力场量,由它与扰动位的泛函关系,就可以对扰动位协方差函数施以相应的算子求其协方差。

对于重力异常,由(2-13)式可知

L.=一…'K=…一?

r:一旦一三,L口:一旦~.

c。Vc△gc以△gc②,=(一导一书,(一杀一书。KcP,②cz珈,

,/,、“+l

cov以)2rpro妻.-oⅢt㈥嘶sy)(2-21)

同理,我们可以推出重力异常,垂线偏差,大地水准面高的协方差函数。互协方差函数,这样,就可以用与重力异常的外推、预估类似的方法,联合各种数据。计算扰动位。

如果待估信号不是扰动位.而是其它量,如:估计S=ST=S9T(P)时,对式(2.17)估计出的扰动位r施以算子S,即;=S71与由

j=【e,,C:,…,C,。】G’,(2—22)计算的结果相同【莫里兹.1984,P58】。其中,e。=SPC户』;此即为最小二乘配置相对

于线性变换的不变特性。

这里.为方便应用,本文根据协方差传搔定律给出重力场观测量间的协方差函数计算

?10?

晟小二乘配置在物理大地测量中的应用

方法

cov(T(e),r(Q)=K(P,Q)

K(P,Qc。v(△g(尸),△g(Q))=c(P,Q)=。,。,.足(P,Q)+2rD,K(P,Q)+rD,K(P,Q)+,,4

cov(Ag(P),f(Q))=(一D,K(P,Q)一二置(P,Q))二

rv

cov(f(P),f(Q))=K(P,Q)÷

cov(4(P),f(Q))=一D。K(P,Q)二-

cov(q(P),f(Q))=一DaK(P,Q)—上

rCOS妒∥

cov(4(P),孝(Q))=D。D,郴,Q)南

卜(卯),袍))=。j足(P'Q’万未而

lc。V(刁(P),刁(Q))=。五足(P,Q)i了尹i忑1i面

lcov(△g(P),艘))=吖。甩P'Q)+吾础,Q))专

|cov(△g(P)'机Q))=-OX(D≯(P,Q)+2-,K(P,Q))赢1

l^

(2-23)

这里:Ag表示重力异常,善,r1分别为垂线偏差的南北分量,f为高程异常:,(r‘)

表示P(Q)的地心距。,,(,‘)表示P(Q)的正常重力值,D。=三,Xt可为r.p和五。

当需要与大地水准面高有关的协方差时.可直接采用对应的与f有关的协方差函数公式计

算。

第二章最小二乘配雯理论

§2.4最小二乘配置方法的数学本质

2.4.1泛函基础

最小二乘配置方法的数学解释涉及到Hilbcrt空间理论,需要在Hilbert空间进行讨论。为此,先简要地介绍一些泛函分析中要用到的概念。

度量空间:对偶(Xd),其中X是一个集合.d是X上的一个度量,即d是定义在(X×X)上,对所有x,Y,z∈X,满足:

I.d是非负的实值

II.d(x,y)=0§J=Y

Illd(x,Y)=d(y,z)(对称性)

IV,d(x,Y)sd(x,z)+d(z,),)(三角不等式)

赋范空间:用范数定义了度量的矢量空间。

从赋范空间X到赋范空间Y的映射叫做算子.而从赋范空间X到标量域R或C的映射nq做泛函,这里,R指实数集合,C指复数集合。

所谓矢量空间X上的范数,对x∈X.其范数记为lIxlI,满足:

I.『Ixll≥0

II.恻=0营x=0

III.{{删=p㈣

IVIIx+yll<-llxll+Uyll

可证若定义d(x,y)=肛一川,则d是X上的度量,这种情况下.称d是由范数诱导出的度量。对赋范空间.范数是其代数结构,因此,赋范空间X也记为(X,d),简记为X。内积空间:定义了内积的矢量空间X为内积空间,或准希尔伯特空间(Pro.HilbcrtSpace)。内积指由XxX到标量域K上的映射,即。对X中的每一对矢量x,Y。都有一个标量。记为(x,Y),与之对应.且帆,Y,z∈X和标!ta,有下式成立:

I.(x+y,z)=(x,z)+(_y,z)

II.(麟,力=盘(x,力

?12?

晟小二乘配置在物理大地测量中的应用

111.(x,Y)=(y,x)

IV.(z,对≥0。且(z,z)=0§x=0

完备的内积空间即H曲ert空间,记为H。

函数K(P,Q)为H的还原核函数,它满足:

?足(只Q)∈H.当9固定时:

?f(Q)=(,(P),K(P,Q)),,W∈H

还原核函数K(P,Q)也简称为核函数。可以证明:核函数K(P,Q)是对称的

K(P,Q)=K(Q,P);正定的:∑∑丸^K(只,0)≥0【莫里兹,1984]。

£f

对可分的Hilbert空间,核函数K(P,Q)可以用一个完备的标准正交系妒(尸)表示;

K(P,Q)=∑仍(尸)妒,(Q)(2—24)

i’I

正交分解:设H、Hl均为Hilbert空间,HIcH,Vx∈H,必有唯一的y∈HI和:∈H1I.使得工=Y+z成立。

2.4.2在具有再生核的Hilbert空间(砌(HS)中的最,j、--乘平差

如果取H为在BjerhaIⅡmar球外调和的函数的Hilben空间,假设H中韵核函数具有如同扰动位的空间扔方差函数的形式,即

…,;纠玎‘驰∽弦zs,式中:k。≥0:P和Q的球坐标表示为(r,0,五),(r’,o’,名),y为P.Q问的球面距

离,由(2—8)式确定。

取下列函数为基函数:

其中

妒。(尸)=

?13?

_R

.S

V』下J

R—r

R一,

...

..,.匿压

第二章最小二乘配置理论

雠牡^嘣忡)fcos叫m2l。

由加法公式可以推得:∑蛾(P)识(Q)=K(P,Q):记五n*(护,^)为五。,j。(9,兄)

|’I

为S…

对H中的函数,f=∑∑(瓦瓦。+瓦。瓦。).g=∑∑(己。瓦。+恧。L),

内积如下定义;(,,g):乙—2n;_+1J己-%瓦。L+嚣。砭。)。若令g(P):K(e,Q),则有

n=0“月m=0

e-=熹瓦,E。=熹夏。(2-27)

(,(尹),K(P,Q)),=(,,g),=∑∑(互。疋,+瓦L)=,(Q)(2-28)

^lOm=O

这样,扰动位的空间协方差函数可以看作H中的调和核函数。

一方面,由式(2.18)可以将最小二乘解手看作是基函数

妒,(P)=妒.K(P,Q)线性组合,即

于(尸)=∑既仍(P)(2.29)

恤l

如果设E是由这Ⅳ个函数生成,记M={仍)”Ⅲ,则E=span{M},于∈E。

从另外一个角度讲.已知Ⅳ个观测值

‘=厶丁,(i=1,2,...,Ⅳ)(2—30).定义了一个超平面D。在一个有限维Hilbert空间R4中,它实际上定义了一个补Ⅳ维(/,/一Ⅳ维)的子空间,所有与观测值‘相容的r一定满足关系式(2.30).因此,所有可能的解于就一定在这个超平面上。

这两点论述表明最4'---乘配置解是空间D和E的交点。可以证明,D上E,即我们

?14?

最小二熏配置在物理大地测量中的应用

由公式(2—17)导出的线性最小方差估值于(P)与D正交,这样,在所有可能的解f中必有归J{sIl司I;证明如下:

令r’满足工;T1=0,则于=r‘+于

1Ifo=(争+r’,}+r’)=(争,争)+z(},,1)+(r,r’).而(于,r’)=。

.‘恫|2=|币+州闩哪。

即最小二乘配置解于具有最小范数的性质。

仍然采用上文中的记法,

B=【t】,K=K(P,Q),I=[f,】,b=暇】:

记观铡值的协方差阵为CⅣ.预估信号与观测值之间的协方差阵为C。.则

【B(BK)7b=I

1于:(B彤)rb

.’于=(BK)7(B(BJ!{:)7)~I=(BK)7(B(BK)7)~BT=PT(2—31)其中:P=(BK)’(B(BK)7)~B,是正交投影葬子,P:珏哼E

这样,就找到了晟小二乘配置解在Hilbert空间中的几何解释。

?15?

第三章最'b--乘配置在局部重力场逼近中的应用

第三章最'b--乘配置在局部重力场逼近中的应用

由最小二乘配置的基础理论可知地球重力场外部扰动位的空间协方差函数对晟小二乘配置是不可缺少的,根据定义求定这个函数需要知道空间每一点的扰动位值,但如果是这样,我们就已经知道了外部重力场而无须配置了。因此,我们只能从形式上得到一个与扰动位协方羞函数作用相似的协方差函数.在本文以后的内容中,如不特殊说明协方差函数均指与真协方差函数作用相似的解析表达式。

§3.1协方差函数模型

局部重力场中的一般情况是,在我们所需要考虑的计算区域内,有相对来说精度和分辨率比较高的重力数据,而对外区则知道的比较少。在这种情况下,我们所采取的方法是用一个先验的全球重力场位模型模拟外区的数据,而内区则采用实际观测值,在具体实现上,即所谓的“移去一恢复”(Remove.Restore)方法。在建立局部协方差函数模型的过程中,也要用到这样一种思路,以重力异常为例,先讨论其对应的全球协方差函数模型。

对重力异常的全球协方差函数模型,可由扰动位的全球协方差函数模型推得:

c(P,O)=弘争2聃。sP,(3?1)

这里:c。=垒云箬t,引入Bjerhmm”球半径,井令c。=以‘百RB)2n4,其中,R。为Bjerhammar球球半径。则有

c(P,Q)=静争2“驰。s¨

po

进一步作替换s=二娑.重力异常协方差函数可表达为:

C(P,Q)=∑y。s”2Po(cosv)

(3-2)

^?0

说明:由(3-2)知道,对于扰动位协方差函数,k。是外部扰动位展开系数所有第月阶的

?16?

最小二乘配置在物理大地测量中的应用

系数的平方和,是扰动位的第n阶谱,也称为扰动位第/1阶阶方差,同理,对重力异常而

言,e。为重力异常第H阶阶方差。以后不再赘述。

如上式所示-理论上如果知道了以(或重力异常阶方差c。)的函数表达式。刚等于

获得重力异常的协方差函数,但这是一个无穷阶的求和式,在使用中.必须进行适当的改化。~种可选的方法是找到它的截断阶数,用有限阶求和来代替无限求和,并作为它的足

够精确的近似.但是,在实际计算中,所需要计算的协方差量巨大。比如:由n.个观测值推估s1个信号值,刚至少需要计算(nl十lⅫI,2+JInl个协方差值,而且对于勒让得多项式只(cos∥).达到一定的阶数后,计算它的稳定性很差.因此,这种方案没有实甩价值;

另一种方案则是拟合出阶方差的函数表达式,并且使得协方差函数可以由一个闭合的解析表达式表达。因为我们已经知道,最小二乘配置理论实用上总是用在局部重力场逼近中,由于局部协方差函数昂显著的特点是三个基本参数,而具有这三个特性参数的拂方差函数模型都是“合理”的,因此,实际上我们先根据协方差函数(阶方差)的一般趋势。得到阶方羞的函数形式,保证协方差函数的无穷项求和可班用闭合表达式表达,然后,根据局部数据拟合这个函数中的参数,达到求定协方差,并减少计算量的目的。

§3.2预备公式

若{“<l’则有母函数公式∑n-Or只‘cos∽=了i萧Ztq-成立口‘力武常次t

.¨一CnRI,,f‘

1988】。

简单起见?记£=√l一2fcosp,+r2,则∑f1P。(cosc∥)=÷。对于将要用到的协方

^-D

差模型,会大量碰到形如

釉力2薹南∥1№们

(3—3、的式子,其中,爿是一个整数,且

{A>o,凡=0

p50,no=1一A

?17?

相关主题
相关文档
最新文档