（大地测量学与测量工程专业论文）区域gps高程转换方法研究.pdf

摘要 尽管g p s f l 皂提供高精度的三维坐标，然而由于g p s 测定的高程是相对于参考椭球面 的大地高，而不是测量中通常所需要的相对于似大地水准面的正常高。因此，如何将大 地高转换为正常高成为人们关注的一个焦点问题。 在采用几何法进行区域高程异常拟合的研究中，利用g p s 水准资料构建空间高程 异常变化模型，最为常用的是函数模型法和基于函数模型与随机模型的综合拟合法，二 者均有较好的拟合效果，但建模精度与函数模型及随机模型的选择有关。支持向量机 ( s v m ) 与最d x - 乘支持向量机( l s s v m ) 理论的提出，为数据拟合研究提供了更为 新颖的方法。基于上述研究背景，本文主要做了以下工作： 1 对函数模型法中的多面函数中心节点的选择做了深入讨论，提出了具有“位置 自适应”匹配特点的k 均值聚类法。实验计算的结果表明，无论高程异常变化平缓区 域还是起伏明显区域，均有较好的拟合效果。 2 基于s v m 理论，在高程异常变化复杂的条件下，对改善s v m 模型拟合精度的 回归算法进行了深入研究，提出了残差修正支持向量机回归算法，通过山区高程异常拟 合实验，证明了该算法的有效性。 3 针对l s s v m 模型参数选择难问题，提出将遗传算法( g a ) 引入到l s s v m 模 型参数的优化选择中，避免了试算过程中的繁琐工作，提高了计算效率。 4 结合l s s v m 对小样本，非线性数据处理的优势，提出了“在线最小二乘支持 向量机回归 算法，以提高模型的解算效率。 关键词：多面函数：k 均值聚类；支持向量机；残差修正；最小二乘支持向量机；遗传 算法； a b s t r a c t a l t h o u g hg p sc a r l p r o v i d eh i g h - p r e c i s i o nt h r e e d i m e n s i o n a lc o o r d i n a t e s ，b u ta sar e s u l t o fg p se l e v a t i o n ，i t sr e l a t i v et ot h ee a r t hr e f e r e n c ee l l i p s o i dh i g h ，r a t h e rt h a nh i g h - n o r m a li s u s u a l l yr e q u i r e di n r e l a t i o nt ot h eq u a s i - g e o i d t h e r e f o r e ，h o ww i l lt h eg p se l e v a t i o n c o n v e r t e dt oh i g h n o r m a lt ob e c o m eaf o c u so fa t t e n t i o n u s i n gt h eg e o m e t r i cm e t h o di nt h eh e i g h ta n o m a l yf i t t i n gr e g i o n a ls t u d y ，u s e i n gg p s a n dt h eh e i g h to fl e v e li n f o r m a t i o n ，t h em o s tc o m m o n l yu s e di saf u n c t i o nm o d e la n d f u n c t i o nm o d e lb a s e do ns t o c h a s t i cm o d e l sa n di n t e g r a t e dt ob el e g a l b o t ha r eg o o df i t t i n g e f f e c t ，b u tt h ea c c u r a c yo fm o d e l i n ga n df u n c t i o nm o d e la n dt h es e l e c t i o no fs t o c h a s t i c m o d e l s s u p p o r tv e c t o rm a c h i n e ( s v m ) a n dl e a s ts q u a r e ss u p p o r tv e c t o rm a c h i n e ( l s s v m ) t h e o r yf o rt h ed a t af i t t i n gp r o v i d eam o r en o v e la p p r o a c h b a s e do nt h ea b o v e r e s e a r c hb a c k g r o u n d ，t h i sp a p e rh a sd o n et h ef o l l o w i n g ： 1 t h ef u n c t i o nm o d e lo ft h em u l t i q u a d r i c sc e n t r a ln o d et oc h o o s et od oa ni n - d e p t h d i s c u s s i o n ，p u tf o r w a r dt h e ”a d a p t i v el o c a t i o n ”t om a t c ht h ec h a r a c t e r i s t i c so ft h ek - m e a n s c l u s t e r i n gm e t h o d e x p e r i m e n t a lc a l c u l a t e dr e s u l t ss h o wt h a tt h ea b n o r m a le l e v a t i o nc h a n g e s ， w h e t h e rf l a to ru n d u l a t i n gr e g i o nm a r k e dr e g i o n ，t h ef i t t i n gr e s u l t sa r eg o o d 2 b a s e do nt h es v m t h e o r y ，a b n o r m a lc h a n g e s i nt h ee l e v a t i o no fc o m p l e xc o n d i t i o n s ， t oi m p r o v et h ea c c u r a c yo fs v mr e g r e s s i o nm o d e lc o n d u c t e di n - d e p t hs t u d yo fa l g o r i t h m s p r o p o s e da m e n d m e n tr e s i d u a la l g o r i t h mf o rs u p p o r tv e c t o rm a c h i n er e g r e s s i o nb yf i t t i n gt h e e x p e r i m e n t a lo fm o u n t a i nh e i g h ta n o m a l y , d e m o n s t r a t et h ee f f e c t i v e n e s so ft h ea l g o r i t h m 3 c h o i c eo fm o d e lp a r a m e t e r sf o rl s s v md i f f i c u l tp r o b l e m ，g e n e t i ca l g o r i t h m ( g a ) i n t ol s s v mo p t i m a lc h o i c eo fm o d e lp a r a m e t e r s ，t h es p r e a d s h e e tt oa v o i dt h ec u m b e r s o m e p r o c e s so fw o r k i n g t oi m p r o v et h ec o m p u t a t i o n a le f f i c i e n c y 4 l s s v mc o m b i n a t i o no fs m a l ls a m p l e ，n o n l i n e a rd a t ap r o c e s s i n ga d v a n t a g e s ，t h e o n l i n el e a s ts q u a r e ss u p p o r tv e c t o rm a c h i n er e g r e s s i o n ”a l g o r i t h m ，i no r d e rt oe n h a n c et h e e f f i c i e n c yo ft h em o d e ls o l v e r k e y w o r d s ：m u l t i q u a d r i c ；k - m e a n sc l u s t e r i n g ；s u p p o r tv e c t o rm a c h i n e ；r e s i d u a le r r o r ；l e a s t s q u a r e ss u p p o r t v e c t o rm a c h i n e s ；g e n e t i ca l g o r i t h m ； i i 论文独创性声明 本人声明：本人所呈交的学位论文是在导师的指导下，独立进行研究工 作所取得的成果。除论文中已经注明引用的内容外，对论文的研究做出重 要贡献的个人和集体，均已在文中以明确方式标明。本论文中不包含任何 未加明确注明的其他个人或集体已经公开发表的成果。 本声明的法律责任由本人承担。 论文作者签名： 汤纸 1 聪卧飞_ 论文知识产权权属声明 本人在导师指导下所完成的论文及相关的职务作品，知识产权归属学 校。学校享有以任何方式发表、复制、公开阅览、借阅以及申请专利等权 利。本人离校后发表或使用学位论文或与该论文直接相关的学术论文或成 果时，署名单位仍然为长安大学。 ( 保密的论文在解密后应遵守此规定) 论文作者签名：稚艮 导师签名：嵇氦奄 长安大学硕士学位论文 第一章绪论 1 1g p s 高程转换的研究意义 全球定位系统( g p s ) 是美国国防部建立的第二代卫星导航定位系统，具有全天候 作业、可提供三维坐标、定位精度高、观测时间短和自动化程度高等诸多优点，已在大 地测量、工程测量、变形监测等众多领域得到了广泛的应用【1 1 。随着数据处理技术的不 断改善，g p s 已具有相当高的定位精度，其短基线平面控制网相对精度可达1 0 一s - l o - 9 ， 垂直分量的精度也已达毫米级1 1 。然而由于g p s 技术测定的是大地高日大地高，而不是工 程建设中常用的正常高常高，尽管大地高与正常高之间存在如下的函数关系 】： h m 常高= k 地高+ 善或善= 厶k 常高一厶k 地高 ( 1 i i ) 其中孝为似大地水准面高度( 也称高程异常) ，即可以通过高程异常进行相互转换。然 而，由于高精度的高程异常值获取难度较高，至使高精度的大地高转换成正常高( 水准 高) 时的精度受损，这在一定程度上抑制了g p s 的应用。于是获取高精度的高程异常， 成为从根本上解决g p s 定位技术无法直接提供高精度正常高的问题，对取代传统的水准 测量方法，提高工作效率，真正实现g p s - - - 维定位功能具有十分重要的实际应用意义。 1 2 高程系统中起算面及相互关系 为了能够更清楚地理解各高程系统之间的关系，有必要首先介绍一下与高程定义相 关的各参考面及高程系统。 1 大地水准面、正高与正高系统 大地水准面是静止海水面向大陆延伸所形成的不规则的封闭曲面，是地球重力场中 的一个等位面，只要给定一点的重力位值，就可以唯一确定过该点的等位面，它是一个 物理曲面，也是与地球最为密合的特殊等位面2 训。 正高是某点沿铅垂线方向到大地水准面的距离，通常用h m 高表示，以大地水准面为 基准面的高程系统称为正高系统。 2 似大地准面、正常高与正常高系统 任一点位的实际重力值与地球内部质量有关，且随着深入地下深度的不同而不同， 而地球内部质量的分布及密度又难以知道，因此很难获取实际重力位值，通常用正常重 1 第一章绪论 力来代替实际重力，由此推求得到的高程即为正常高常高。 似大地水准面则是地面点沿垂线向下量取正常高常高所形成的连续曲面，是一个 与大地水准面极为接近的基准面【2 训，它不是水准面，而是用于描述地球形状的一个几 何曲面，不具有实际物理意义，是用以计算的辅助面。一般来讲，它与大地水准面不完 全吻合，差值为正常高与正高之差。在海洋面上，似大地水准面与大地水准面重合，在 平原和山区，两者的差距与点位的高程有关【4 1 。 我国采用正常高高程系统作为国家高程系统，由传统水准测量测取的地面点高程属 于正常高系统。 3 参考椭球面、大地高与大地高系统 由于地球内部质量分布不均匀，各点位重力( 垂线) 方向不规则变化，使得大地水 准面实际上是一个形状复杂的不规则曲面，难以将地面上的测量结果归算到大地水准面 上进行相关计算。 为了达到有效、方便的目的，需要寻找一个和大地水准面非常相近且可用简单数学 公式表达的几何形体来代替大地水准面。在测量上选用了椭球绕其短轴旋转而成的参考 旋转椭球来近似大地水准面【1 4 】。 大地高坎地高即定义为空间点沿法线方向到参考椭球面的距离。由参考椭球面为起 算面的高程系统为大地高系统。空间坐标系统通常以观测点纬度、经度及观测点到参考 椭球面的距离大地高来表示，即( b ，l ，呶地高) 。 4 正高、正常高与大地高之间的关系 从上面高程系统的定义可知，不论正常高系统还是正高系统，与大地高系统差别主 要在于起算面的不统一，即( 似) 大地水准面和参考椭球面的不一致。图1 1 为各起算面 关系示意图。 似大地水准面 大地水准面 参考椭球面 图1 1 高程系统中各起算面关系示意图 2 长安大学硕士学位论文 大地高、正常高、正高存在转换关系如下： 日大地高= h i e 常高一善h 大地高= h i e 高一n ( 1 2 1 ) 式中善和分别表示似大地水准面和大地水准面到参考椭球面的距离，又称为似大地水 准面和大地水准面高度，如上图所示。 g p s 观测得到的大地高是以参考椭球面为起算基准的大地高系统f 1 】，而日常生产、 生活、服务中，采用的高程基准是以似大地水准面为起算面的正常高系统或者以大地水 准面为起算面的正高系统2 4 1 。 1 3 国内外研究综述 g p s 观测值经数据处理后能够提供高精度的大地高，欲将高精度的大地高应用于测 量生产及人们的日常生活，国内外学者都对g p s 高程转换技术进行了深入的研究。美国、 加拿大采用正高系统，对北美地区的正高与大地高的转换进行深入理论分析，并通过生 产实践建立了北美区域g 9 5 0 1 ( 美国) 、g s d 9 5 ( 加拿大) 大地水准面模型。欧洲国家 采用正常高系统，法国、德国等欧洲国家协同合作，利用其基础观测成果建立了e g g 9 7 欧洲区域似大地水准面模型。在本世纪之初，我国也推出了新一代的c q g 2 0 0 0 似大地水 准面模型【4 1 。 不论大地高转换到哪种高程系统，目前转换方法归结为两类：重力法和几何法【2 4 1 。 重力法中，主要以测点点位的重力异常值作为边值数据，通过斯托克斯积分或莫洛金斯 基积分来求解大地水准面高或似大地水准面高【3 , 4 1 。积分计算过程中考虑不同波长的贡献 大小不同，为了提高拟合精度通常采用移去恢复技术，分别对中长波及短波分量独立求 解，最后叠加得出结果。同时在计算过程中，采用快速傅里叶变换技术，来降低问题求 解复杂度【3 捌。一般来讲，通过重力法确定的似大地水准面，在重力资料分布密集及密 度均匀的条件下，解算精度很高，但求解过程相对复杂。但当重力资料分布不均且密度 不足时，建立的似大地水准面模型精度较低。我国新一代c q g 2 0 0 0 就存在模型东西部精 度不均的缺点。由于重力资料不仅获取成本高，且作为国家保密测量成果，其使用受到 严格限制。因此，尽管重力法有相对较高的转换精度，但难以推广。 几何法即通过若干个采样点，分别测定其大地高与正常高( 水准高) ，建立高程异 常与坐标位置之间的函数关系，以此来确定大地水准面，由此即可对区域范围内的点进 行g p s 高程转换。由于受采样点密度，地形复杂度及函数模型选择的多因素影响，相比 3 第一章绪论 重力法，应用几何法确定全球，大范围的似大地水准面模型精度相对较低，但对中小区 域范围，其拟合精度可以达到g p s 高程转换的精度要求，加之应用几何法拟合原理简单、 建模直观、资料获取成本较低，所受限制条件相对宽松等特点，现已广泛应用到各项工 程领域。相关应有研究主要有：小区域范围的高程基准控制网【5 1 、地表垂直沉降监测6 1 、 代替山区区域四等水准测型7 1 、跨海高程传谢8 1 、水下航道水准测量【9 1 。 根据数据处理手段的不同，几何法拟合空间高程异常变化又可分为数学模拟法和机 器学习法。 1 数学模拟法 根据采用的模型不同，数学模拟法可分为函数模型拟合法和随机模型拟合法【l o 】。其 中函数模型拟合法一般求定高程异常的系统性或规律性变化的趋势，常用的有平面函 数、二次曲面函数、样条函数、多面函数。一般来讲平坦区域建模时，最为简单方便的 就是函数模型法中的曲面函数，具体应用中函数形式通常取平面函数、相关平面函数、 二次曲面函数。在空间高程异常变化复杂的区域，通常可以考虑采用建模精度相对较高 的多面函数模型法、样条函数拟合法等【n 】。区域范围较大时，可以考虑分区拟合方案， 如移动曲面拟合法、有限元方法等区域拟合方案，对整个数据拟合区域拟合结果通过各 种分区方案拼接条件完成【1 2 15 1 。 由于地球质量密度内部分布不均匀，高程异常的空间变化不可能完全通过函数模型 来模拟，可以将统计学中的随机模型法引入到空间高程异常拟合，即统计模型拟合法。 随机模型法是根据各离散数据之间的统计特性对有限观测值进行统计，以确定各数据点 之间的相关性模型求定其随机性变化。随机模型法的特点是，对随机变化的信息拟合较 好，但不能很好的反映系统特性的信息。 通常，空间高程异常变化不仅表现出随机特性，同时有一定的系统性。大地测量与 地球物理界研究者提出了结合两者特点的综合模型法、克里金估值法、最d , - 乘配置法 等【16 1 。最d x - 乘配置法中，不仅确立了系统部分，同时调整信号部分与之相匹配，在数 据拟合推估中取得了较好的拟合效果【1 眈们。 2 机器学习算法 上世纪4 0 年代以来，随着计算机技术发展，模拟人类思维方式的机器学习算法发 展迅速，最具代表的就是神经网络( n e t u r a ln e t w o r k ，n n ) 2 1 】。由于n n 对非线性变 化数据具有较好的拟合效果，因此有学者提出用神经网络进行高程异常拟合【2 2 2 3 1 。考虑 到n n 存在模型的相关参数设置不易，泛化能力弱等问题，相关学者又结合智能优化领 4 长安大学硕士学位论文 域另有算法遗传算法( g a ) ，提出通过g a 选择n n 的参数建立遗传算法神经网络 高程转换模型【2 4 1 ，结果表明采用遗传算法( g a ) 优化选择n n 模型参数是非常有效的。 随着最优化理论、人工智能、计算机技术的进一步发展，v a n p i k 经过多年的研究， 于上世纪九十年代中期完善了他提出的统计学习理论( s t a t i s t i c a ll e a r n i n gt h e o r y , s l t ) ， 并给出了支持向量机( s u p p o r tv e c t o rm a c h i n e ，s v m ) 的概念。由于统计学习理论( s l t ) 给出了控制模型精度的另一标准结构风险最小( s t r u c t u r a lm s km i n i m i z a t i o n , s 刚) ，提高了模型的泛化能力。s v m 最早应用于模式识别领域，后来通过引入占不敏 感损失函数概念，使之应用于函数回归问题。解决非线性问题时，通过统计学习理论 ( s l t ) 中的核函数( k e r n e lf u n c t i o n ) 引入降低了模型求解的难度。s v m 模型不仅对 非线性数据具有很强的回归特性，同时对样本之间相关信息几乎没有要求，具有建模方 便的特点【2 5 。2 8 1 。 构建高程异常空间变化支持向量机( s v m ) 模型时，以若干个联测了水准高程的 g p s 点的大地坐标佃，三) 或平面坐标( x ，y ) 为模型的输入量，对应的高程异常值 孝= 畋地高一鲰准高为模型输出量，建立高程异常空间变化s v m 回归模型，进而通过该 模型内插计算出区域内任意点的高程异制2 9 1 。 在“建模数据 数目增大的条件下，由于s v m 模型的求解是一个凸二次规划问题 ( q p ) ，存在求解效率较低问题。针对s v m 这一缺点，s u y k e n s 提出最小二乘支持向量 机( l e a s ts q u a r e ss u p p o r tv e c t o rm a c h i n e ，l s s v m ) ，他将占不敏感损失函数改化为占2 的 形式，使得问题求解转化为求解一个线性方程纠3 0 3 1 1 。对比s v m 求解q p 问题，求解 线性方程组的l s s v m 计算效率明显提高。 不论s v m 、l s s v m 用于模式识别还是回归问题，合适模型参数的选择都是一个重 点、难点问题。模型参数选择，通常可以引入遗传算法( g a ) 、蚁群算法( a c o ) 、人 工免疫算法( i c s a ) 等智能计算完成参数的优化选择【3 2 】。 1 4 本文研究的内容 显然，函数模型选择不当或对需要转换的数据相关性( 实际不确定先验信息) 描述 不够，都会影响g p s 高程转换精度。本文在对比分析常用的区域高程异常拟合法基础 上，重点研究了多面函数中心节点的选取方法。考虑到高程异常转换本身的复杂性，提 出应用基于小样本学习，具有较强非线性映射能力的支持向量机进行高程异常转换，并 气 第一章绪论 对相关问题进行了研究探讨，具体研究内容概括如下： 1 总结了常用g p s 高程异常拟合方法中的曲面拟合法、多面函数法、最小二乘配 置法。着重分析了多面函数拟合法中的核函数中心节点选择问题，提出了具有“位置自 适应”匹配特点k - m e a n s 聚类法选择多面函数中心节点，并通过实例验证它的有效性。 2 阐述了分类s v m 、回归s v m 的基本原理，及几类具体类型的回归s v m 算法 和相互关系，提出了残差修正回归s v m 算法( r e s i d u a le r r o rs u p p o r tv e c t o r sr e g r e s s i o n ， r es v r ) 。 3 分析最小二乘支持向量机( l s s v m ) 算法和遗传算法( g a ) 基本理论，并将 g a 引入到l s s v m 参数优化选择中。最后，分析了作者提出的“在线l s s v m 算法。 4 讨论回归支持向量机( s v r ) 、残差修正回归支持向量机( r es 己) 、基于遗传 算法( g a ) 的最小二乘回归支持向量机( l s s v r ) 在高程异常拟合中的应用，并通过 实例分析，得出相应的结论。 5 总结本文的研究工作和取得的成果以及存在的不足之处，并对后续工作进行展 望。 6 长安大学硕士学位论文 第二章高程异常拟合常用方法 几何转换法是区域g p s 高程转换常用方法，它是利用一定数目联测了水准高的g p s 点，通过数据拟合方法建立高程异常与平面坐标或大地坐标的关系，从而插求出其它未 知点的高程异常值，进而达到g p s 高程转换的目的。根据拟合时采用的模型不同，可 分为函数模型拟合法和随机模型拟合法。其中函数模型拟合法一般求解某种规律性、趋 势性的成分，如曲面函数拟合法、多面函数拟合法、样条函数拟合法等【l 。随机模型法 则求定随机性变化的信息，典型的有k r i n g i n g 拟合法、最d , 、- - 乘配置法等【l6 1 。本章主 要介绍几种常用的高程异常拟合法。 2 1 曲面拟合法 建立空间高程异常变化模型，最为简便且常用的方法是函数模型法。在高程异常变 化趋缓的小区域，多采用平面模型或二次曲面模型【1 1 3 3 3 4 1 。地形起伏明显，空间数据变 化复杂的情况下，理论上可采用高次多项式拟合，达到数值逼近效果。然而，由于其稳 健性能较差，使得拟合结果可能产生“空间振荡 现象，数值分析中称之为“龙格现象 【3 5 1 ，因此不宜采用高次曲面函数。 2 1 1 曲面拟合基本原理 高程异常拟合中，若任一点高程异常值善，可用曲面函数善= f ( x i ，乃) 表示，则任 一点高程异常值轰，可以用该曲面函数求解。由于受观测误差及函数形式的影响，很难 找到能够正确反映各点高程异常的曲面函数，因此在确定曲面函数时，要求善= f ( x i ，m ) 能够与高程异常变化实际曲面最佳逼近【3 3 3 4 1 。则曲面函数模型可表示为 手= 孝+ ( 2 1 1 ) 确立高程异常善函数形式的关键是求解模型的系数项参数，则上式可表示为 孝= 邛+ a ( 2 1 2 ) 其中孝= 缶，彘己r 为观测值，f 为设计系数矩阵，通常由曲面函数的形式确定，为 待求的曲面函数系数项参数，为高程异常的真误差。误差方程可表示为 v = f l l - 孝 ( 2 1 3 ) 7 第二章高程异常拟合常用方法 根据最小二乘原理，解得 = ( f r ，) 1 ( f r f ) ( 2 1 4 ) 将上式解得的代入( 2 1 2 ) 式，便可获得确定的曲面函数模型，从而插求出其它 未知点的高程异常值。 2 1 2 曲面函数形式的选择 曲面拟合中，函数形式的选择很关键，常用的函数形式主要有平面函数、相关平面 函数、在x 方向或者在y 方向上变化明显的函数形式、二次曲面函数形式及高次曲面函 数形式【3 6 1 。 1 平面函数( 线性内插) f ( x ，y ) = p o + 一+ 及y ( 2 1 5 ) 其中f l i ( i = o ，1 ，2 ，) 为未知参数，此时要求公共点至少为3 个。 2 相关平面函数( 双线性内插) f ( x ，y ) = 风+ 一+ 厦y t 屈砂 ( 2 1 6 ) 其中f l ，( i = 0 , 1 ，2 ，3 ) 为未知参数，此时要求公共点至少为4 个。 3 沿x 或y 方向变化明显的曲面函数 f ( x ，y ) = p o + i l l x + f 1 2 y + f 1 3 x y + f 1 4 x 2 ( 2 1 7 ) f ( x ，y ) = p o + f l , x + f 1 2 y + f 1 3 x y + f 1 4 y 2 ( 2 1 8 ) 其中f l ，( i = 0 , 1 ，2 ，3 ，4 ) 为未知参数，此时要求公共点至少为5 个。 4 二次曲面函数 f ( x ，y ) = 风+ i l l x + f 1 2 y + f 1 3 x y + f 1 4 x 2 + f l s y 2 ( 2 1 9 ) 其中f l ，( i = 0 ，1 ，2 ，3 ，4 ，5 ) 为未知参数，此时要求公共点至少为6 个。 以上几类曲面函数中，在数据变化趋势缓和、空间数据分布规则时，平面函数拟合 效果较好：相关平面函数在数据随x ，y 方向具有波动但变化趋势单一的情形下拟合效果 较好；二次曲面函数则能够拟合变化趋势稍微复杂，但规律性明显的数据；在空间数据 变化复杂、分布不规则时，上述几类形式的曲面函数拟合效果都不是非常理想。 r 长安人学硕士学位论文 2 2 多面函数法 上世纪7 0 年代，美国的h a r d y 教授在解决航天器外形设计时，提出了由多个简单 曲面叠加，来完成对复杂的外形模拟的多面函数法【3 7 l 。他认为任何一个圆滑的数学表面， 总可用一系列有规则的数学表面的总和以任意精度逼近。 2 2 1 多面函数的基本原理 根据数据逼近原理，已知曲面函数f ( x ，y ) ，称满足方程 ( 厂( 薯， ) 一缈( 玉，咒) ) 2 = m i n ( 2 2 1 ) 的函数缈( x ，y ) 为曲面函数f ( x ，y ) 的逼近函数。根据h a r d y 多面函数的原理，逼近函数 伊( x ，y ) ，可表示为1 3 3 8 4 1 1 缈( x ，j ，) = j 兰= l 岛q 凡工，y ) ，( x ，y j ) ) ( 2 2 2 ) 式中，q 为核函数，为待定系数，上限“为核函数个数( 中心节点个数) ，o ，y ，) 为 选取的中心节点。理论上讲，核函数可以取任一简单函数，为了计算方便，一般取对称 型、距离型的核函数，如图2 1 所示。 |、 万2 = 1 0 么夕朱 心 铲= 一夕 否2 = o 钟形核函数 钵形核函数 图2 1 两种类型核函数【3 1 】 根据核函数形状，常用的核函数分为两类：钟形函数和钵形函数。前者随 着h ，m 的增大，o - - + 0 ，后者随着h ，l y i 的增大，q o o 。( 2 2 3 ) 式( 2 2 5 ) 式列 出了空间数据拟合中常用三种核函数具体表达形式【3 8 舶1 钟形核函数： 倒双曲函数 q = ( ( 工一x ，) 2 + ( y y ，) 2 + 万) 一1 7 2 ( 2 2 3 ) 9 第一二章高程异常拟合常用方法 钵形核函数： 正双曲函数 三次曲面型 q = ( ( x - - x + ( y y ) 2 + 万) 1 7 2 ( 2 2 4 ) q = ( ( x - x j ) 2 + ( y - y j ) 2 ) 3 7 2 + 万2 ( 2 2 5 ) 以上三式中，万表示核函数平滑因子，起到调节核函数形状的作用。 设有力个空间观测数据( 置， ，) ，i = 1 , 2 ，刀，则可建立误差方程组 v l2 二层g ( ( 五，乃) ，( 而，_ ，) ) 一 表示成矩阵形式，则有 其中 2 盖乃g ( ( 而，y 2 ) ，( ，y o j ) ) 一乞( 2 2 6 ) 屹= j ，q ，( ( 毛，儿) ，( 而，) ) 一厶 ，穹1 。 v = a p l v = ( k ，圪圪) r = ( 。，：夕。) r 三= ( ，乞，厶) r ( 2 2 7 ) ( 2 2 8 ) ( 2 2 9 ) ( 2 2 1 0 ) fq 。( ( 工l ，y 1 ) ，( x y o ) ) ，q ( ( x i ，y i ) ，( x o py o 。) ) 1 彳：l ；f ( 2 2 1 1 ) xhl i 幺( ( ，y 。) ，( x o ，y 。) ) ，q 。( ( _ ，y 。) ，( p ，y o 。) ) j 式中，u 多面函数核函数中心节点的个数，应用最小二乘原理，可解得 夕= ( 彳r 彳) 一1a r l ( 2 2 1 2 ) 将解得的系数夕代入( 2 2 2 ) 式，则可获得多面函数模型，进而求得插值点对应的 函数值。 2 2 2 多面函数相关参数设置 从多面函数的基本原理可以看出，多面函数模型的建模精度，由核函数形式、平滑 1 0 长安人学硕士学位论文 因子设置，以及及中心节点的选取决定。 1 核函数的选择 多面函数模型中核函数的选择是关键，理论上讲任一简单的函数都可，但核函数选 择通常与数据自身有很强相关性，相关研究也总结出了一些一般性的结论【3 8 枷】。黄立人 等给出的判断核函数优劣原则主要有：首先保证问题有解，其次要对测点具有良好的拟 合效果，并且拟合值与实测值之间不存在系统偏差，再者要求未测点的推估值较为平稳， 不至于在预测点偏离数据点时，内插值有急剧的变化，直观的理解便是在插值时避免数 据振荡现象出现。目前认为具有较好拟合效果的核函数有倒双曲函数( 2 2 3 ) 式、正双 曲核函数( 2 2 4 ) 式、三次曲面核函数( 2 2 5 ) 式【3 8 4 0 】。 2 平滑因子的确定 平滑因子万的作用是改变核函数的形状，核函数确立后，需要考虑平滑因子万及核 函数中心节点的选择。从图2 1 核函数的具体形状可以看出，万设置不同时，虽然函数 大体走向、趋势不会改变，但会对核函数的细节产生明显的影响。经验认为万应取一个 比较小的正数或者零【3 8 - 4 0 ，空间d e m 数据内插实例发现，万在一定程度上可以放大。 3 中心节点的获取 核函数与平滑因子选定后，中心节点的数目及位置选择是一个相对重要的问题。中 心节点位置，一般取分布在研究区域内均匀分布的特征点【3 8 4 0 1 ，数据分布均匀时也可以 选择非严格的“数据格网化 后的网格交叉点。选择均匀分布的特征点时，多面函数的 拟合效果非常理想，但建模过程中，需要事先对数据进行分析，找出这些特征点。一般 情况下，寻找特征点的过程比较繁琐。借鉴模式识别领域k 均值聚类【4 2 - 4 4 1 ( k - m e n a s ) 的算法特性，本文提出基于k - m e n a s 算法选择多面函数中心节点。 k - m e n a s 聚类的主要思想是：将个对象数据集构造k 个分组，给出一种初始分 组方法，通过反复迭代来改变分组，使得每次改进后的分组方案都较前一次的好。 k - m e n a s 聚类属于聚类方法的一种基本划分，聚类过程中，通常采用“误差平方和函数 作为准则函数，其函数形式为 j 。= j 主- 。左。 一朋川2 ( 2 2 1 3 ) 其中，m ，代表第_ ，个聚类中心，代表第类各数据点。 k - m e n a s 聚类过程中，考虑初始点的选择对聚类的速率会产生较大影响【4 3 朋】。为了 提高k 均值聚类的搜索效率，取“格网划分后均匀分布的“格网交叉点为聚类算 第二章高程异常拟合常用方法 法起始搜索点。具体计算步骤如下： 1 根据采样数据最小、最大值确定数据分布范围，按照一定的间隔划分获取格网 交叉点。依照欧氏距离，即二范数条件确定能够“响应”到“建模数据 的交叉点，这 样可以排除孤立点对聚类的影响； 2 以响应到的交叉点作为对应数目个类的中心值，将数据赋予距离其最近的类； 3 计算每个类中数据平均值，以平均值为新的聚类中心； 4 重复执行步骤2 、3 直到各类中心不再发生变化。 2 3 最小二乘配置法 最小二乘配置简称为配置法，由t k r a r u p 和h m o r i t z 先后进行系统的研究，在地 球重力场的研究中，配置法模型由两部分组成：随机部分和非随机的系统部分。高程异 常由于受到空间重力梯度变化、地球物质密度分布影响，在体现系统特性的同时也表现 出一定的随机性 1 6 , 1 9 , 4 5 1 。 2 3 1 最小二乘配置法基本原理 最d - - 乘配置法的一般线性函数模型描述为t 4 5 】 三：彳y + 否y + ( 2 3 1 ) 式中三为n 维列向量，彳为甩“维设计矩阵，x 是u 1 维非随机向量，y 是信号向量， 是咒维观测噪声列向6 o 其中】，= 霎， ，否= b ，。 ，b 是已测点信号的设计矩阵，o 是相应的零矩阵，s 是力个已测点信号向量，s 是任意多个未测点信号向量，可假设为 q x l 维。 很明显，观测量可分为三部分：a x 由一种规律性变化而产生的，称为倾向部分； b y 的变化因随机信号s 、s 的随机扰动引起，这种扰动受周围环境扰动的影响【4 5 1 ；观 测噪声，可以是随机独立也可以非随机独立。与( 2 3 1 ) 式对应的误差方程为 y = 彳它+ 百矿一l ( 2 3 2 ) 式中，y 是观测值残差向量， 夕吲， 雪、已测点和未测点的信号估值。假设观测 噪声和信号都己中心化，则随机先验信息相关关系如下 1 2 e ( ) = 0 = e ( a a r ) e c y ，= 。 y = e c y 】，r ，= 主羔乏善 设信号与观测噪声不相关，则y = 0 ，其中，a 为观测噪声阵，s 、s 、好分别为 已测点、未测点、以及已测点和未测点方差协方差阵，其值由给定的信号协方差函数求 得。 假设单位权方差2 = 1 ，只是观测噪声权矩阵，b 是信号权矩阵，则有， p = 列= 阿兰1 旺3 依照拉格朗日乘数法，构造目标函数如下 q = y r 只y + 矿r b 矿一2 k 7 ( 么岩+ 百p y 一三) = m i n ( 2 3 4 ) 上式分别对矿、j 、矿求一阶偏导数并令等于零，则可解得 y = 一巧1 k = 一a k ( 2 3 5 ) 矿= 巧1 百r k = y 西7 k ， ( 2 3 6 ) a r k ：0 ( 2 3 7 ) 顾及a y = 。，】，； ； ，及否= c 召，叨，将以上三式及c 2 3 2 ，式联立求解得到 j = ( z ：1 彳) 一1 a r ：1 l ( 2 3 8 ) = s b r ：1 ( 一彳克) ( 2 3 9 ) 雪= s s b r ：1 ( 三一詹) ( 2 3 1 0 ) 其中：a + b z s b r 。在( 2 3 1 ) 式中，若彳= o 或者x = o 则简化为不带非随机参数 的推估模型，若百：0 ，或者】，= 0 ，则转化为参数平差模型1 9 4 5 1 。 2 3 2 最小二乘配置法中协方差函数的拟合 配置法中，信号协方差阵s 、的确定是拟合推估的关键，通常它们又是根据协 方差函数计算得到的。因此数据拟合中，协方差函数确定的准确与否关系到算法的稳定 第二章高程异常拟合常用方法 性，进而影响空间数据的拟合精度【4 5 ，4 6 1 。 在空间数据拟合中，随机信号的先验统计信息的相关性表现为欧氏距离特性，即距 离越近，相关性越强；距离越远，相关性越弱；当距离达到一定大小时相关性消失。因 此，可以通过“随距离增加相关程度降低的函数来描述随机信号的这一空间数据特性。 实践中可依据一定的数理实测点来拟合协方差函数，进而获取信号协方差。 对于空间数据拟合，尤其是g p s 高程转换，对高程异常的倾向变化部分选择合适 的倾向表达式，根据最小二乘原理计算倾向参数估值 膏聒= ( 彳r 只么) 一1a r 只三 ( 2 3 1 1 ) 从原始观测数据中减去倾向面，得到残差s 为 a s = 三一衍岱 ( 2 3 1 2 ) 其中a s = s + a ，即残差中包含了信号和观测噪声，通常认为信号与噪声不相关，并把 “倾向面”看作一个常数，因此有工= 砖= + s 。 假设s = 置，是，a 岛a 昌，a s , 】r ，参照统计学，测量平差中相关期望，方差的 定义，可以计算出残差的数学期望估计值趣，方差估计值仃a ，( 0 ) 、协方差估计值 吼p 0 ) 如下 去孥西( 2 3 1 3 ， 1 一 ( o ) 2 刍吾( 岛一仫) 2 ( 2 3 1 4 ) 计算两点之间间隔为f 0 的协方差估值 1广1 ，( ”o ) 2 者善【- ( 昂叫如) ( & i - 。) - - u a s ) j q 3 j 5 其中n ，为两点间间隔为f 的点对个数。 假定是等精度的独立观测，则是一对角阵，其主对角线元素为吼2 ，即 ( o ) = 2 ，( f = 0 ) ( 2 3 1 6 ) 信号协方差阵主对角线匕的元素是信号方差，即 1 4 长安大学硕上学位论文 s ( 0 ) = 2 ，( f = o ) ( 2 3 1 7 ) 顾及工= k = + s ，则有 。( 0 ) = ( o ) + s ( 0 ) ，0 = 0 ) ( 2 3 1 8 ) 。( f ) = 5 ( f ) ，p 0 ) ( 2 3 1 9 ) 根据统计得到的方差仃a j ( o ) 和不同间隔的协方差盯a j ( r o ) ，选择合适协方差函数 模型，如高斯函数模型s ( f ) = z s ( 0 ) e x p ( k 2 f 2 ) ，得到一个方程组，根据最小二乘估计拟 合出协方差函数模型中的估计参数【1 7 】。 2 。4 实验分析 前面对常用高程异常拟合方法中的曲面函数法、多面函数法、最d x - - 乘配置法的理 论体系做了介绍，下面从实验的角度来分析这几类数据拟合方法在区域范围g p s 高程 转换中的应用。介绍实验之前，首先对本文实验结果统计表中公共使用的符号所代表的 含义，及其计算方式作出如下说明： r m s r 、r m s e 分别代表拟合点和检核点的均方根误差； r m s r = k 宰k 】，其中k 表示拟合误差向量，啊为拟合点数目； g m a e = 【吒幸】伤，其中k 表示检核误差向量，刀：为检核点数目； m a x a b s ( v ：) 表示检核点最大绝对偏差值； 本文所有实验结果中，高程异常单位均为n l 。对于个别计算实验中，特有符号代表 的含义，将在实验中做具体说明。 2 4 1 平坦地区g p s 高程转换实验 实验区域a 属于地形变化相对平缓的地区