GPS基线向量网平差VB程序设计(共28页)

1、精选优质文档-倾情为你奉上GPS基线向量网平差程序设计前 言GPS技术以其定位精度高，观测自动化,不需测站间通视及网型与精度关系不大的优势,已成为建立城市及工程控制网的主要技术手段之一。而与常规地面网相比,GPS控制网的数据处理有其自身的特点，由于基线向量是不可独立于坐标系而存在的特殊观测值，所以在平差时或平差后必须转入测区所在的坐标系统。本论文讨论了GPS基线向量的转换和平差问题及工程控制测量实用的方法，并运用VB程序设计语言完成了大地空间直角坐标向大地坐标的转换、大地坐标向高斯平面坐标的转换、二维基线向量网平差的功能。1GPS原理1.1GPS的简介全球定位系统（全局位置系统，GPS）是美国

2、从上世纪70年代开始研制，历时20年，耗资200亿美元，于1994年全面建成的利用导航卫星进行测时和测距，具有在海、陆、空进行全方位实时三维导航与定位能力的新一代卫星导航与定位系统。它是继阿波罗登月计划、航天飞机后的美国第三大航天工程。如今，GPS已经成为当今世界上最实用，也是应用最广泛的全球精密导航、指挥和调度系统。它主要由三大子系统构成：空间卫星系统、地面监控系统、用户接收系统。1.2GPS定位原理 GPS系统采用高轨测距体制，以观测站至GPS卫星之间的距离作为基本观测量。为了获得距离观测量，主要采用两种方法：一是测量GPS卫星发射的测距码信号到达用户接收机的传播时间，即伪距测量；一是测量

3、具有载波多普勒频移的GPS卫星载波信号与接收机产生的参考载波信号之间的相位差，即载波相位测量。采用伪距观测量定位速度最快，而采用载波相位观测量定位精度最高。通过对颗或颗以上的卫星同时进行伪距或相位的测量即可推算出接收机的三维位置。按定位方式，GPS 定位分为单点定位和相对定位（差分定位）。单点定位就是根据一台接收机的观测数据来确定接收机位置的方式，它只能采用伪距观测量。相对定位（差分定位）是根据两台以上接收机的观测数据来确定观测点之间的相对位置的方法，它既可采用伪距观测量也可采用相位观测量。在定位观测时，GPS定位分为动态定位和静态定位。若接收机相对于地球表面运动，则称为动态定位。若接收机相对

4、于地球表面静止，则称为静态定位。GPS定位的实质是根据GPS接收机与其所观测到的卫星之间的距离和所观测卫星的空间位置来求取接收机的空间位置,而这些又是根据GPS卫星发出的导航电文计算出的包括位置、伪距、相位和星历等原始观测量,通过计算来完成的.理论上卫星到接收机的距离为：=c(-)上式中:为第i颗GPS卫星到接收机的距离;c为光速;为接收机相对于统一的时间基准接收到卫星发出信息的时刻;为卫星相对于统一的时间基准发出信息的时刻.事实上,在卫星所发出信息的准确性及其传送还会受到许多诸如发出信息时刻的卫星轨道偏差、电离层与对流层的延迟效应、卫星时钟和接收机时钟与统一的时间基准之间的偏差等因素的影响,

5、造成一定程度的偏差.根据计算GPS卫星到接收机距离的方法,大体可以分为伪距测量定位和相位测量定位两种基本定位方法.1.3GPS发展趋势1991 年的海湾战争中，装在大衣口袋中的GPS 接收机为无地图沙漠作战发挥了巨大作用。在“盟军行动”中，把惯导GPS 集成系统装入导弹和制导导弹，使命中精度达到9 m，而且使机载炸弹具备了在夜间和恶劣天气条件下的精确打击能力。由此可见，GPS早已成为高技术武器平台不可缺少的关键组成部分。在新世纪以及未来军事战争中GPS 将发挥更加巨大的作用。这样的形势下迫使GPS技术必须要有新的突破。经过不懈的努力钻研，如今经取得些成绩。就导航定位卫星技术主要发展趋势如下。（

6、1）采用创新轨道设计欧洲多年来从未中断对导航定位卫星的研究、论证。在第一代中，有“全球导航卫星系统”（GNSS）以及“欧洲静止轨道导航重叠业务系统”（EGNOS）等，它们都是结合利用GPS和静止轨道通信卫星的方案。在第二代中，目前采用创新轨道设计的“伽利略”方案被认为是能够实现最少投入而达到理想应用目的的最佳方案。它既是独立系统，又有开放性特点，可与GPS兼容。这种系统还将在民航选择最佳航线、飞机安全进场着陆等领域有新的应用突破。（2）美国大力开发抗干扰和干扰技术GPS集成到高技术武器平台，使GPS 应用概念发生全新变化。为防止地方干扰，美国在将从2005年发射的第7颗GPS-2F卫星上开始使

7、用新型信号结构。这样，除更加保密外，还可实现6dB 的信号 干扰比的改善。为此，正在研制不受干扰和欺骗的 GPS 接收机应用模块（GRAM）和选择利用抗欺骗模块（SAASM），同时装有这两种模块的接收机被称为“国防部高级GPS 接收机”（DARG）。美国还在开发抗干扰的军事伪系统（Millitary Pseudolites），它可为地域发射GPS差分信号，以改进信号捕获并提高质量。为保护军用飞机使用GPS，美国还在开发微带自适应天线阵列。为使敌方不能使用GPS，美国已开发出GPS干扰机，只有可口可乐瓶大小的干扰机可使敌方无法接收GPS信号。（3）进入21 世纪，GPS在各方面的应用都将加强和发

8、展。一、在综合服务系统中的应用在全球地基GPS连续运行站（约200个）的基础上所组成的IGS（Internationtol GPS），是GPS连续运行站网和综合服务系统的范例。它无偿向全球用户提供GPS 各种信息，如GPS精密星历、快速星历、预报星历、IGS站坐标及其运动速率、IGS站所接收的GPS信号的相位和伪距数据、地球自转速率等。这些信息在大地测量和地球动力学方面支持了无数的科学项目，包括电离层、气象、参考框架、精密时间传递、高分辨的推算地球自转速率及其变化、地壳运动等。二、在电离层监测中的应用GPS在监测电离层方面的应用，也是GPS 空间气象学的开端。太空中充满了等离子体、宇宙射线粒子

9、、各种波段的电磁辐射，由于太阳常在秒钟内抛出百万吨量级的带电物，电离层由此而受到强烈的干扰，这是空间气象学研究的一个对象。通过测定电离层对GPS信号的延迟来确定在单位体积内总自由电子含量（TEC），以建立全球的电离层数字模型。三、在对流层监测中的应用GPS在监测对流层方面的应用，早期主要是由于轨道误差影响定位精度，而且早期的GPS基线相对来说比较短，高差不大，因此对对流层的研究没有给予很大的重视。直到近期由于GPS轨道精度大大提高后，当对流层折射已经成为限制GPS定位精度提高的一个重要障碍时，才开始认真的对对流层的监测研究。我们可以假设在一个高程基本为零的地区，并且如果接收机所接收的GPS信号

10、是从天顶方向传来的话，那么其延迟就可以达到2.22.6 这一量级，而2小时内这一延迟变化可达10 不是少见的，所以IGS 分析中心所提供的对流层参数是采用2小时间隔一次。也正是由于这个实际情况，对流层折射要顾及其随机过程的变化来加以模型化。四、在卫星测高仪中的应用多路径效应是GPS 定位中的一种噪声，至今仍是高精度GPS 定位中一个很不容易解决的“干扰”。过去几年利用大气对GPS 信号延迟的噪声发展了GPS 大气学，目前也正在利用GPS定位中的多路径效应发展GPS测高技术，即利用空载GPS作为测高仪进行测高。它是通过利用海面或冰面所反射的GPS信号，求定海面或冰面地形，测定波浪形态，洋流速度和

11、方向。通常卫星测高或空载测高所测的是一个点，连续测量结果在反向面上是一个截面，而GPS测高则是测量有一定宽度的带，因此可以测定反射表面的起伏（地形）。据报告，试验时空载平面安装2台GPS接收机，1台天线向上用于对载体的定位，1台天线向下，用于接收GPS在反射面上的信号。美国在海上作了测定洋流和波浪的试验。丹麦在格凌兰作了测定冰面地形及其变化的试验。五、在卫星追踪技术中的应用卫星对卫星的追踪（SST）技术的实质是高分辨率 的测定两颗卫星间的距离变化，一般它分为两类，即高低卫星追踪和低低卫星追踪。前一类是高轨卫星（如对地静止卫星，GPS卫星等）追踪低轨（LEO）卫星或空间飞行器，后一类是处于大体为

12、同一低轨道上的两颗卫星之间的追踪，两颗卫星间可以相距数两千米，这两类SST 技术都将LEO 卫星作为地球重力场的传感器，以卫星间单向或双向的微波测距系统测定卫星间的相对速度及其变率。这一速度的不规则变化所反映的信息中， 就包含了地球重力场信息。卫星轨道愈低，这一速度变化受重力场的影响愈明显，所反映重力场的分辨率也愈高。2GPS基线向量网平差模型2.1GPS的坐标转换技术GPS定位成果属于WGS-84坐标系，而实用的成果属于国家大地坐标系，因而GPS成果的坐标系转换是必不可少的。2.1.1工程GPS网和地面网之间的坐标转换模型众所周知,两个空间直角坐标系之间的三维坐标转换一般采用含有7个转换参数

13、的布尔莎、莫洛金斯基或范士模型,对于坐标差转换而言,因其不涉及3个平移转换,可只采用余下的4个转换参数,但这仅适用于点位在两个坐标系中均具有精确的三维坐标的情况。因以往用常规技术测设的地面网属于二维坐标系,由地面点的正常高加上概略的高程异常所得的已知椭球面上的大地高只有数米的精度,由此而进行的三维坐标转换只能达到米级的精度。在工程GPS网转换中,是将GPS基线向量转换为高斯平面上的二维基线向量。因此,适用的GPS网坐标差转换模型就是仅有尺度和方位2个参数(其余2个旋转参数为零)的简化范士模型。2.1.2从投影变换方面保持高斯平面上边长尺度的一致性利用GPS技术建立平面控制网,最后需得到GPS网

14、点的高斯平面直角坐标,这就必须选择一个椭球面作为过渡,并须选定某一经线作为高斯投影的中央子午线,这两者都关系到高斯平面上的边长尺度。GPS网应与地面网从空间边长投影方面取得相一致的边长尺度,与此相应的数据处理方法就是采用与地面网边长归算的高程基准面(常称为投影面)较为接近的区域性椭球面;采用与地面网中央子午线在位置或经度上相同的经线作为GPS网点投影变换的中央子午线。2.1.3按附合网还是按独立网进行平差定位采用GPS技术来加强改善原有地面控制网,如何合理对待、处理地面网的已知点须根据网的用途地面网的实际精度作认真细致的分析、比较而定。若固定地面网一个起始点外,或再固定地面网的一个起始方位角,

15、都不改变平差后GPS网的相对构形,而属于独立网平差,亦可称之为纯GPS网的无约束平差。只要如上所述选定合适的椭球面和高斯投影中央子午线,纯GPS网在高斯平面上的边长尺度就能与地面网坐标系中应有的边长尺度相一致。其优点是,归算到地面网坐标系后仍保持GPS网形的高精度,有利于在已有控制网的坐标基准中鉴定和改善原有的地面网。但由于只是在位置基准点上仍保持原有坐标值,所以其他GPS点上的高斯坐标必与原有坐标值不相一致,而且离位置基准点越远,坐标较差也越大。使近万平方公里的大城市,最大坐标较差可达30。另一种处理方法是对GPS网进行附合网平差,对于重合于多个已有控制点的GPS点,除固定基准点的三维坐标外

16、,其余重合点均固定其二维坐标,在平差的同时就能获得由GPS网转换到地面网坐标系的尺度及方位旋转参数。由于在重合点上保持了原有的二维坐标,非重合点上与原有坐标的差异也会小些,迎合了生产单位不愿改动原有点位坐标的要求。但因难以估计的地面网固有误差会使GPS相对网形遭受扭曲,从而不能保持GPS网应有的高精度,由此建立的GPS网充其量只能达到与原有地面网大致相当的精度。2.2GPS基线向量的转换2.2.1基线向量转换的过程GPS数据处理一般分为基线解算和基线网平差两个阶段。第一阶段的基线处理是从原始相位观测值入手，经过粗加工、预处理、基线解算等步骤，最后获得高精度的基线向量及其协方差阵。第二阶段是利用

17、第一阶段解算出的三维基线向量组成基线向量控制网，进行各种条件检核、平差检验，而后把属于WGS-84坐标系的基线向量转换国家或地方坐标系下，并以国家或地方坐标系的高级控制点的平面坐标为基准对基线网进行约束平差。GPS基线向量网与地面无论在三维参心空间直角坐标系中或三维参心大地坐标系中进行平差，都可得到良好的三维空间位置的平差结果，是严密的解法。特别是在三维大地坐标系中进行联合平差，可将表示平面坐标信息分量同高程位置坐标分量很方便的区分开，并进而简便的转换成工程测量实用的控制成果。因此，这两种联合平差方法，特别是后者在大地测量和工程测量中得到广泛应用。但在这些空间坐标系中进行联合平差时，必须知道满

18、足一定精度要求的地面点的大地高或相应的大地水准面差距N。对于N，目前一般都是采用某种地球重力场模型通过模拟计算的方法得到，但某钟地形地理条件下，在一些地区还很难获得满意的结果。故工程测量中，为避开这个目前尚难以解决的实际问题，现代工程GPS基线向量网与地面网联合平差还常常采用在二维坐标系中进行的办法。二维坐标系可以是参考椭球面也可以是高斯投影（或某种工程施工坐标平面）。由于工程测量的范围一般都比较小，往往都是采用高斯平面直角坐标系（或施工平面直角坐标系），因而联合平差往往是在平面直角坐标系中进行，但也可在参考椭球面二维曲面坐标系中进行。这就是说，在联合平差中，首先应将GPS基线向量观测值及地面

19、网观测值按照一定的数学关系式将它们化算到椭球面或进一步归算到统一的高斯投影平面坐标系中，然后在椭球面或平面直角坐标系中建立观测量的误差方程式，固定量的约束条件方程，并确定它们在该坐标中的随机模型，再进行联合平差，从而得到适宜于工程测量需要的控制测量成果。2.2.2本设计基线向量转换摸型本设计是将GPS基线向量及其随机模型原封不动地按照数学关系式转换到高斯平面直角坐标系中，再结合地面已知数据进行联合平差。首先根据有关GPS观测值和地面观测值在统一坐标系中合并的基本理论，将GPS网中的一个固定点坐标及GPS观测得到的基线向量，求得各点的三维空间直角坐标，并利用公式： (1) (2) (3)将其转换

20、成大地坐标，然后，舍去大地高H，利用，通过高斯坐标正算公式计算高斯平面直角坐标，最后再按取坐标差的办法，得到平面直角坐标系中的GPS基线向量及。在利用它们按数学关系式列立误差方程式及约束条件方程。对GPS基线向量的随机信息，按照：的顺序，也进行方差-协方差的传播并转换到二维平面直角坐标系中。对地面网观测值（如水平方向，斜距等）则按常规地面网平差方法，将其归化到高斯投影平面坐标系上，从而获得属于该坐标系的地面观测值。最后，将归算到高斯平面坐标系中的GPS基线向量及地面观测量这两类观测值在该坐标系中建立联合平差数学模型，按其进行联合平差。2.3二维基线平差的误差方程及求解2.3.1列误差方程GPS

21、与地面数据的二维联合平差中的GPS观测量是经过归算的二维基线向量，。而地面数据有平面边长和平面方向观测值以及作为起始数据的已知地面点坐标，固定的平面边长和坐标方位角。 二维联合平差以网中各点在地面系统的平面坐标为未知参数，记为: (1)由于二维GPS基线向量的归算未顾及其尺度差异，同时虽作了方位归算，但还会存在残余的方位差异，因此还应设定尺度参数和方位参数。也可以说二维GPS基线向量与地面系统的平面坐标属于两个不同的平面坐标系。平差后应有 : (2)上式中也就是二维GPS基线向量的观测方程，其误差方程为: (3)写成矩阵形式为： (4) 上式中，矩阵为二维基线向量（坐标差）改正数组成的列向量，

22、 为单位阵，为转换参数 的系数阵，为常数列向量。对于一端为固定点的基线向量 其误差方程式为： (4)a (4)b将（2）式进一步简化为： (5)上式中 ， 。 2.3.2法方程式的组成及解算由于各基线向量观测值之间互相独立 因而可分别对每个基线向量观测值的误差方程式组成法方程 将单个法方程的系数阵及常数项加到总法方程的对应系数项和常数项上去。对应于（3）式和的法方程式分别为： (6) (7)设网中有条基线向量 则总法方程的系数阵和常数项为： ； (8)总的法方程可写为: NdX-U=0 (9) 在程序设计中 法方程组成可有下述两种算法:一 、按点累加法设法方程系数阵按照点分块 则分块阵为： ；

23、 在法方程系数阵中， 对角块是第i点与周围其它控制点连测的基线权阵之和, 而非对角元（ij）是第i点到第j点连测基线权阵之和乘以-；而常数项是第点连测基线权阵与（4）式中的对应乘积之和。二、 按基线累加法首先对基线向量的误差方程进行改化，使各误差方程互相独立，使权阵变为单位阵，即用 乘以相应的误差方程，得到改化后的误差方程，然后累加法方程，从而使算法简单化。解算法方程后得到未知数为： (10)各待定点坐标平差值为： 2.3.3精度评定单位权方差估值为:上式中,为基线向量个数，n平差未知数个数，为s为网中点数,平差未知数的方差估值为: 。 主要的程序介绍3.1VB程序设计语言简介Visual B

24、asic(简称VB)，它是Microsoft公司推出的一种Windows应用程序开发工具。由于它是在Windows平台下设计应用程序的最迅速、最简捷的工具之一，而且具有简单易学、操作方便、功能强大等特点，已成为普通用户首选的程序设计语言。“Visual”指的是采用可视化的开发图形用户界面（GUI）的方法，一般不需要编写大量代码去描述界面元素的外观和位置，而只要把需要的控件拖放到屏幕上的相应位置即可方便设计图形用户界面；“Basic”指的是BASIC语言，因为VB是在原有的BASIC语言的基础上发展起来的。VB包含Microsoft Excel、Microsoft Access等众多Window

25、s应用软件中的VBA都使用VB语言，以供用户进行二次开发；目前制作网页使用较多的VBScript脚本语言也是VB的子集。利用VB的数据访问特性，用户可以对Microsoft SQL Server和其他企业数据库在内的大部分数据库格式创建数据库忽然前端应用程序，以及可调整的服务器端部件。利用ActiveX(TM)技术，VB可使用如Microsoft Word字处理器、Microsoft Excel电子数据表及其他Windows应用程序提供的功能，甚至可直接使用有VB专业版或企业版创建的应用程序和对象。用户最终创建的程序是一个真正的.EXE文件，可以自由发布。3.2程序流程图建立三维基线向量数据文

26、本 读入GPS三维基线向量数据 根据三维基线向量观测值和固定点坐标推算待定点在WGS-84坐标系下的近似坐标 将待定点坐标转换到高斯投影平面直角坐标系中，并求出它们的坐标差，作为二维基线向量 根据地面已知点坐标及二维基线向量，将所有点坐标转到国家坐标系下 根据误差方程组法方程 解法方程，输出最终平差结果3.主要程序代码3.3.1读GPS基线网数据其数据结构如下 2，105说明 边长比例误差（单位为十分之一毫米），中央子午线经度（度）y,B068,-.925,.847,.320 说明 识别符号，固定点点名，固定点的空间坐标B064 测站点点名B068,-421.318,-15.016,-202.

27、943 基线目标点点名 , 三维基线向量end主要程序代码如下： CommonDialog1.ShowOpenfname = CommonDialog1.FileName 将用户在打开对话框中选择的文件名对变量fname赋值。If fname Then 若无此判断当对话框中选择取消时、下面赋值语句将出错。Set ts = fso.OpenTextFile(fname) 将fname作为文本文件打开,并设置句柄。 k = 0: p = 0: h = 0: ds(j) = 0: j = 0: d = 0 j是测站数累计变量,k是已知点累计变量,ds()为基线累积变量 Do 后测型循环，结束循环的条

28、件是读到文件结束符end q = ts.ReadLine 读一行，置入q。 q = Trim(q): i = 1: 删除q可能有的前导和尾随空格，i是工作变量。 m(i) = InStr(q, ,) 查行中第一个逗号的左数位置,并保存在整形数组变量m(i)中 Do While m(i) 0 前测型Do. Loop循环，成立条件是该行字符串中有逗号。 tr(i) = Mid(q, m(i - 1) + 1, m(i) - m(i - 1) - 1) 提取指定位置开始的指定数目字符。 i = i + 1 ：m(i) = InStr(m(i - 1) + 1, q, ,) 从上一个找到的逗号位置起，

29、查找下一个逗号的位置。LoopIf m(i) = 0 And i 1 Then tr(i) = Right(q, Len(q) - m(i - 1) 处理一行中最后一个逗号后的字符串。以下部分是将存储在数组变量 m(i)中的字符分类存放到基线向量、已知坐标、网型信息等数组中If p = 0 Then 读到的是文件第一行。mk = tr(1): L0 = tr(2): p = 1 提取边长比例误差mk（用于边长定权）和中央子午线经度L0 ,并使该句以后不能再执行。 ElseIf m(1) = 0 And q end Then j = j + 1: dm(j) = q: ds(j) = ds(j

30、- 1) 该行中没有逗号,但又不是结束符，则一定是测站点名。将读出的字符串赋值到点名数组变量dm(j)。If p 0 And m(1) 0 Then 不是第一行，并且该行中有逗号分割的多个字串。If tr(1) = y Then 有识别符号y则这一行不是基线向量，而是固定点坐标值。k = k + 1: ym(k) = tr(2): xo(k) = Val(tr(3): yo(k) = Val(tr(4): zo(k) = Val(tr(5) 存储固定点点名及坐标值。 ElseIf tr(1) y Then 累计基线目标点号到点名数组变量dn(j)中ds(j) = ds(j) + 1:dn(ds

31、(j) = tr(1): dx(ds(j) = Val(tr(2): dy(ds(j) = Val(tr(3): dz(ds(j) = Val(tr(4) 提取三维基线向量.Loop Until Trim(q) = end 数据读取完毕。3.3.2空间近似坐标推算（1）累加所有点名 For i = 1 To ds(cds) 依次访问所有测站（cds）的目标方向. p = 0 设识别变量. For j = 1 To d 依次访问所有测站（d）.If dm(j) = dn(i) Then p = 1 查看目标点是否设过测站，是则对识别变量p赋值1。If p = 0 Then d = d + 1:

32、dm(d) = dn(i) 如 p=0，表明目标点未设过测站,将该点点名加入点名数组 zds = d 将总点数存入模块级变量zds （2）计算所有测点的空间近似坐标 x(1)=100:y(1)=100:z(1)=100 先给一个测站假定一个坐标值。For i = 1 To cds 按测站循环。For j =ds(i-1)+1 to ds(i) 遍访i测站上所有的目标点. h=seqn(dn(j) 调用函数查目标点对应的序号。If Abs(x(i)0 Then 如果测站点已解出，则计算目标点的假定坐标值。x(h)=x(i)+dx(j): y(h)=y(i)+dy(j): z(h)=z(i)+dz

33、(j)计算公式。end ifif Abs(x(h)0 then 如果目标点已解出，则计算测站点的假定坐标值。 x(i)=x(h)-dx(j): y(i)=y(h)-dy(j): z(i)=z(h)-dz(j) 计算公式。 End if Next jNext iDo 后测型循环，当所有点的近似坐标算出则退出循环。s = s + 1 s是循环计数变量，控制循环次数，避免假定坐标计算不出时，进入死循环。 For i = 1 To cds 按测站再次进入循环。 If Abs(x(i)0 Then 在该测站假定坐标已计算出的前提下，计算这条基线另一点的假定坐标。 For j = ds(i - 1) +

34、1 To ds(i) 遍访i测站上所有基线。 h = seqn(dn(j) 调用函数查目标点对应的序号。 If Abs(x(h)0 Then Exit For 如目标点坐标已解出，退出ForNext循环。End ifNext jFor j = ds(i - 1) + 1 To ds(i) 再次遍访i测站上的所有基线。 h = seqn(dn(j) 调用函数查目标点对应的序号。If Abs(x(h)0 Then 如果该目标点坐标已解出，反算k点的假定坐标值。 x(k) = x(h) - dx(j): y(k) = y(h) - dy(j): z(k) = z(h) - dz(j) If Abs(

35、x(h) zds 坐标值已全部算出或虽还有未算出的，但根据循环次数已不能算出时结束循环。For i = 1 To yds 查找固定点坐标For j = 1 To zds 查点名数组中那个点是固定点,用m()数组存其序号If ym(i) = dm(j) Then m(i) = j X1=xo(i)-x(j): Y1=yo(i)-y(j): Z1=yo(i)-y(j) 计算固定点真实坐标与假定的差值。 End If Next jNext iFor i = 1 To zds 按站点循环，将所有点的坐标值都加上坐标改正数。 x(i) = x(i) + X1: y(i) = y(i) + Y1: z(i

36、) = z(i) + Z1Next iFor i = 1 To yds 置入固定坐标值 x(m(i) = xo(i): y(m(i) = yo(i): z(m(i) = zo(i)Next i3.3.3二维基线向量的转换（1）大地坐标B,L计算的程序代码 For i = 1 To zds 按测站循环 B(i) = bth(x(i), y(i), z(i) 调用纬度反算的函数公式，计算i点的纬度值赋给数组变量B() L(i) = ath(x(i), y(i) 调用经度反算的函数公式计算i测站的大地经度,并赋值给数组变量L() Next i一、 Private Function ath(x As

37、Double, y As Double) 经度坐标反算通用公式的函数程序代码。ath = Atn(y / x) If x 0 Then 坐标象限的转换，pi为圆周率。 ath = ath + piElse If y 0 Then ath = ath + 2 * piEnd IfEnd Function二、 Private Function bth(x As Double, y As Double, z As Double) 纬度坐标反算公式通用函数的程序代码。t=0Do 后测型循环。 t=t0 n=a/Sqr(1-e2*sin(t0) e为克拉索夫斯基椭球体的第一偏心率,a为克拉索夫斯基椭球长

38、半轴。 t = Atn(z+n*e2*sin(t0)/sqr(x2+y2)t与t0进行叠代，直到满足精度要求。Loop Until Abs(t-t0)0.1 叠代前后值之差的绝对值小于千分之一秒时退出。bth = t 将叠代值赋给函数。End Function 二、高斯投影正算的程序代码For i = 1 To zds 按站点循环,利用高斯投影正算公式实现坐标转换，并用数组变量gx(),gy()储存平面近似坐标 b1 = Sin(B(i): b2 = Cos(B(i): t = Tan(B(i)B(i)为i点的纬度。 n = c / (Sqr(1 + e2 * b2 2) n为该站点所在卯酉圈

39、的曲率半径，c 为克拉索夫斯基椭球体极点处的子午线曲率半径,e为克拉索夫斯基椭球体的第二偏心率。 x0 = .861 * (B(i) * 180 / pi) - 32005.78 * b1 * b2 - 133.929 * b2 * b1 3 - 0.697 * b2 * b1 5 x0为i点的子午线弧长。 dl = L(i) - L0 * pi / 180 L0（整度）为中央子午线，L(i)为i点的经度。 r = b2 * dl v = e1 * b2 2 A1 = (5 - t 2 + (9 + 4 * v) * v) / 24 A2 = (61 + (t 2 - 58) * t 2) /

40、 720 A3 = (1 - t 2 + v) / 6 A4 = (5 + (t 2 - 18 - 58 * v) * t 2 + 14 * v) / 120 gx(i) = x0 + n * t * r 2 * (0.5 + (A1 + A2 * r 2) * r 2) gy(i) = n * r * (1 + (A3 + A4 * r 2) * r 2)（3）计算二维基线向量For i = 1 To cds 按测站循环，计算二维基线向量并用数组变量kx(),ky()储存。 For j = ds(i - 1) + 1 To ds(i) k = seqn(dn(j) 调用函数查目标点对应的序号

41、。kx(j) = gx(k) - gx(i): ky(j) = gy(k) - gy(i) 计算公式。Next jPrivate Function seqn(str As String) As Integer 由点名查计算序号函数For i = 1 To zdsIf str = dm(i) Then seqn = i 将查到的序号赋给函数名,返回调用处Next iEnd Function3.3.4平面近似坐标计算的程序代码这部分的程序与前面空间近似坐标的计算相同，首先读入平面已知点坐标，格式为：点名，X坐标，Y坐标 end 接着进行平面近似坐标的推算，推算程序和前面相同。3.3.5组法方程的程

42、序代码Dim n1 As Integer, l1 As Double, pp As Double 定义过程级变量 l1 = 0 常数项清零 n1 = 2 * (zds - yds) + 2 未知数数目其中有两个为转换参数未知数。 n2 = n1 * (n1 + 1) / 2 一维存储法方程系数数组上限 ReDim NX(n2), UX(n1) 重新定义法方程系数、常数数组 For i = 1 To cds 按测站循环 For j = ds(i - 1) + 1 To ds(i) 在i测站上按方向循环，求误差方程系数、常数。 ReDim nb(n1) 重新定义误差方程系数数组,并且每循环到一新边

43、长前清零。 h = seqn(dn(j) 查目标点在点名数组中的序号 ss = Sqr(kx(j) 2 + ky(j) 2) 反算边长 pp = 1/(ss + mk * ss * 10 -4) 边长观测值定权,mk为边长比例误差。 cha = charact(i, k) 自定义函数，查测站点i是否已知点，如不是,用k返回i前面有几个已知点 If cha = n Then 测站点i不是已知点d = 2 * (i - k) 1 计算测站i点x未知数在未知数点集中的序号。chu = charact(h, s) 自定义函数，查目标点j是否已知点，如不是,用s返回序号h前面有几个已知点。 If chu

44、 =n then 如果目标点不是已知点，则组法方程。 t=2*(h-s)-1 先计算目标站点j的x未知数在未知数点集中的序号。 nb(t) = 1: nb(d) = -1 nb(n1 - 1) = kx(j): nb(n1) = -ky(j): l1 = x(h)-x(i) - kx(j) 根据误差方程给x坐标系数数组赋值并计算常数项l1，将转换参数固定为未知数点集最后两个位置。 Call equation(nb(), pp, l1) 调用函数组法方程，参数分别是x坐标误差方程系数数组和转换参数数组、权、误差方程常数项。 nb(t + 1) = 1: nb(d + 1) = -1: nb(n1

45、 - 1) = ky(j): nb(n1) = kx(j): l1 = y(h)-y(i) - ky(j) 根据误差方程给y坐标系数数组赋值并计算常数项l1，将转换参数固定为未知数点集最后两个位置。 Call equation(nb(), pp, l1) 组法方程，参数分别是y坐标误差方程系数数组和转换参数数组、权、误差方程常数项。End if 这个条件判别完，进入下一个条件判断。 If chu=y then 如果目标点为已知点则以nb(h) 和nb(h+1)都为零组法方程。 nb(d) = -1: nb(n1 - 1) = kx(j): nb(n1) = -ky(j): l1 = x(h)-

46、x(i) - kx(j) 系数数组赋值并计算常数项。 Call equation(nb(), pp, l1) 调用函数组法方程，同上。 nb(d + 1) = -1: nb(n1 - 1) = ky(j): nb(n1) = kx(j): l1 = y(h)-y(i) - ky(j) 根据误差方程给系数数组赋值并计算常数项。 Call equation(nb(), pp, l1) 调用函数组法方程，同上。 End if End if 当测站不是已知点时，判断并组法方程完毕。 If cha= =y Then 如果i测站为已知点时组则以nb(d)和nb(d+1)为零组法方程。chu = chara

47、ct(h, s) 自定义函数，先查目标点j是否已知点.If chu=n Then 如不是,用s返回序号h前面有几个已知点。 nb(h) = 1: nb(n1 - 1) = kx(j): nb(n1) = -ky(j): l1 = x(h)-x(i) - kx(j) 对误差方程系数数组赋值同时计算常数项。 Call equation(nb(), pp, l1) 调用函数组法方程，同上。 nb(h + 1) = 1: nb(n1 - 1) = ky(j): nb(n1) = kx(j): l1 = y(h) y(i)- ky(j) 对误差方程 y系数数组赋值. Call equation(nb()

48、, pp, l1) 调用函数组法方程，同上。 End If 这种情况判断完毕，进行下一种判断。chu = charact(h, s) 自定义函数，先查目标点j是否已知点，如不是,用s返回序号h前面有几个已知点。If chu=y then 如果目标点也是已知点，则退出循环。 Exit for End if End if Next j 一个目标点的误差方程处理完毕，进入下一个目标方向。 ll = ll + pp * l1 2 累积和方程常数项。 Next i 一个测站的误差方程处理完毕，进入下一测站4实例平差对比与结果分析4.1实例数据准备先将南方GPS adj4.0处理所得的在WGS84-坐标系

49、下经过自由网平差的三维基线向量减去相应的改正数，作为原始数据。做法如下：基 线 名 X Y Z X改正 Y改正 Z改正B064，B068 -421.317 -15.016 -202.943 0.987 0.133 0.068B077，B078 -118.733 -153.881 185.225 -10.399 8.121 4.074改后为：基 线 名 X Y Z B064，B068 -421.318 -15.016 -202.943 B077，B078 -118.723 -153.889 185.221 然后将改过的三维基线向量制成适应于本设计运行的exe文本。下面为本实例的原始数据：2,10

50、5y,B064,-.609,.863,.263 B077B078,-118.723,-153.889,185.221 B085,-306.783,-221.570,192.144 B127,665.204,-76.724,491.886 B129,605.475,2.523,318.949 B085B078,188.040,67.722,-6.907B093,-472.221,-175.392,76.949 B109,-440.175,43.702,-256.151 B089B109,-267.982,156.320,-359.022 B078,360.235,180.330,-109.788

51、 B085,172.189,112.614,-102.875 B093,-300.032,-62.787,-25.939 B093B068,571.449,476.216,-500.935 B109,32.050,219.101,-333.087 B106N050,1039.566,126.420,317.926 B068,246.987,120.199,-73.433 B093,-324.466,-356.008,427.509 B109,-292.406,-136.900,94.431 B109B068,539.390,257.101,-167.862 B078B129,724.238,156.416,133.731 B127,783.968