GPS基线向量网平差VB程序设计

GPS基线向量网平差程序设计

前 言

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

1

1 GPS原理

1.1 GPS的简介

全球定位系统（全局位置系统，GPS）是美国从上世纪70年代开始研制，历时20年，耗资200亿美元，于1994年全面建成的利用导航卫星进行测时和测距，具有在海、陆、空进行全方位实时三维导航与定位能力的新一代卫星导航与定位系统。它是继阿波罗登月计划、航天飞机后的美国第三大航天工程。如今，GPS已经成为当今世界上最实用，也是应用最广泛的全球精密导航、指挥和调度系统。它主要由三大子系统构成：空间卫星系统、地面监控系统、用户接收系统。

1.2 GPS定位原理

GPS系统采用高轨测距体制，以观测站至GPS卫星之间的距离作为基本观测量。为了获得距离观测量，主要采用两种方法：一是测量GPS卫星发射的测距码信号到达用户接收机的传播时间，即伪距测量；一是测量具有载波多普勒频移的GPS卫星载波信号与接收机产生的参考载波信号之间的相位差，即载波相位测量。采用伪距观测量定位速度最快，而采用载波相位观测量定位精度最高。通过对４颗或４颗以上的卫星同时进行伪距或相位的测量即可推算出接收机的三维位置。按定位方式，GPS 定位分为单点定位和相对定位（差分定位）。单点定位就是根据一台接收机的观测数据来确定接收机位置的方式，它只能采用伪距观测量。相对定位（差分定位）是根据两台以上接收机的观测数据来确定观测点之间的相对位置的方法，它既可采用伪距观测量也可采用相位观测量。在定位观测时，GPS定位分为动态定位和静态定位。若接收机相对于地球表面运动，则称为动态定位。若接收机相对于地球表面静止，则称为静态定位。

GPS定位的实质是根据GPS接收机与其所观测到的卫星之间的距离和所观测卫星的空间位置来求取接收机的空间位置,而这些又是根据GPS卫星发出的导航电文计算出的包括位置、伪距、相位和星历等原始观测量,通过计算来完成的.理论上卫星到接收机的距离为：Li=c(tp-ts) ．

2

上式中:Li为第i颗GPS卫星到接收机的距离;c为光速;tp为接收机相对于统一的时间基准接收到卫星发出信息的时刻;ts为卫星相对于统一的时间基准发出信息的时刻.事实上,在卫星所发出信息的准确性及其传送还会受到许多诸如发出信息时刻的卫星轨道偏差、电离层与对流层的延迟效应、卫星时钟和接收机时钟与统一的时间基准之间的偏差等因素的影响,造成一定程度的偏差.根据计算GPS卫星到接收机距离的方法,大体可以分为伪距测量定位和相位测量定位两种基本定位方法.

1.3 GPS发展趋势

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

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

GPS集成到高技术武器平台，使GPS 应用概念发生全新变化。为防止地方干扰，美国在将从2005年发射的第7颗GPS-2F卫星上开始使用新型信号结构。这样，除更加保密外，还可实现6dB 的信号／ 干扰比的改善。为此，正在研制不受干扰和欺骗的 GPS 接收机应用模块（GRAM）和选择利用抗欺骗模块（SAASM），

3

同时装有这两种模块的接收机被称为“国防部高级GPS 接收机”（DARG）。美国还在开发抗干扰的军事伪系统（Millitary Pseudolites），它可为地域发射GPS差分信号，以改进信号捕获并提高质量。为保护军用飞机使用GPS，美国还在开发微带自适应天线阵列。为使敌方不能使用GPS，美国已开发出GPS干扰机，只有可口可乐瓶大小的干扰机可使敌方无法接收GPS信号。

（3） 进入21 世纪，GPS在各方面的应用都将加强和发展。 一、在综合服务系统中的应用

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

二、在电离层监测中的应用

GPS在监测电离层方面的应用，也是GPS 空间气象学的开端。太空中充满了等离子体、宇宙射线粒子、各种波段的电磁辐射，由于太阳常在１秒钟内抛出百万吨量级的带电物，电离层由此而受到强烈的干扰，这是空间气象学研究的一个对象。通过测定电离层对GPS信号的延迟来确定在单位体积内总自由电子含量（TEC），以建立全球的电离层数字模型。

三、在对流层监测中的应用

GPS在监测对流层方面的应用，早期主要是由于轨道误差影响定位精度，而且早期的GPS基线相对来说比较短，高差不大，因此对对流层的研究没有给予很大的重视。直到近期由于GPS轨道精度大大提高后，当对流层折射已经成为限制GPS定位精度提高的一个重要障碍时，才开始认真的对对流层的监测研究。我们可以假设在一个高程基本为零的地区，并且如果接收机所接收的GPS信号是从天顶方向传来的话，那么其延迟就可以达到2.2～2.6ｍ 这一量级，而2小时内这一延迟变化可达10ｃｍ 不是少见的，所以IGS 分析中心所提供的对流层参数是采用2小时间隔一次。也正是由于这个实际情况，对流层折射要顾及其随机过程

4

的变化来加以模型化。

四、在卫星测高仪中的应用

多路径效应是GPS 定位中的一种噪声，至今仍是高精度GPS 定位中一个很不容易解决的“干扰”。过去几年利用大气对GPS 信号延迟的噪声发展了GPS 大气学，目前也正在利用GPS定位中的多路径效应发展GPS测高技术，即利用空载GPS作为测高仪进行测高。它是通过利用海面或冰面所反射的GPS信号，求定海面或冰面地形，测定波浪形态，洋流速度和方向。通常卫星测高或空载测高所测的是一个点，连续测量结果在反向面上是一个截面，而GPS测高则是测量有一定宽度的带，因此可以测定反射表面的起伏（地形）。据报告，试验时空载平面安装2台GPS接收机，1台天线向上用于对载体的定位，1台天线向下，用于接收GPS在反射面上的信号。美国在海上作了测定洋流和波浪的试验。丹麦在格凌兰作了测定冰面地形及其变化的试验。

五、在卫星追踪技术中的应用

卫星对卫星的追踪（SST）技术的实质是高分辨率 的测定两颗卫星间的距离变化，一般它分为两类，即高低卫星追踪和低低卫星追踪。前一类是高轨卫星（如对地静止卫星，GPS卫星等）追踪低轨（LEO）卫星或空间飞行器，后一类是处于大体为同一低轨道上的两颗卫星之间的追踪，两颗卫星间可以相距数两千米，这两类SST 技术都将LEO 卫星作为地球重力场的传感器，以卫星间单向或双向的微波测距系统测定卫星间的相对速度及其变率。这一速度的不规则变化所反映的信息中， 就包含了地球重力场信息。卫星轨道愈低，这一速度变化受重力场的影响愈明显，所反映重力场的分辨率也愈高。

5

设法方程系数阵按照点分块 则分块阵为：

?N11?N21 N???????Nn1?U1?N1n????N22?N2n??U2? ； U????????????Nn2?Nnn???Un??N12?在法方程系数阵中， 对角块Nii是第i点与周围其它控制点连测的基线权阵之和, 而非对角元Nij（i?j）是第i点到第j点连测基线权阵之和乘以-１；而常数项Ui是第ｉ点连测基线权阵与（4）式中的Li对应乘积之和。

二、 按基线累加法

首先对基线向量的误差方程进行改化，使各误差方程互相独立，使权阵变为单位阵，即用

P 乘以相应的误差方程，得到改化后的误差方程，然后累加i法方程，从而使算法简单化。

?为：X???N?1U (10) 解算法方程后得到未知数X~??Xi??Xi0??Xi?各待定点坐标平差值为： ????0???~? Y?i??Yi???Yi?2.3.3 精度评定

单位权方差估值为: ?2?oVTPV/??m?n?s?? ?上式中,ｍ为基线向量个数，n平差未知数个数，为s为网中点数,平差未知数X2的方差估值为: Di??0N?1 。

11

３ 主要的程序介绍

3.1 VB程序设计语言简介

Visual Basic(简称VB)，它是Microsoft公司推出的一种Windows应用程序开发工具。由于它是在Windows平台下设计应用程序的最迅速、最简捷的工具之一，而且具有简单易学、操作方便、功能强大等特点，已成为普通用户首选的程序设计语言。

“Visual”指的是采用可视化的开发图形用户界面（GUI）的方法，一般不需要编写大量代码去描述界面元素的外观和位置，而只要把需要的控件拖放到屏幕上的相应位置即可方便设计图形用户界面；“Basic”指的是BASIC语言，因为VB是在原有的BASIC语言的基础上发展起来的。

VB包含Microsoft Excel、Microsoft Access等众多Windows应用软件中的VBA都使用VB语言，以供用户进行二次开发；目前制作网页使用较多的VBScript脚本语言也是VB的子集。

利用VB的数据访问特性，用户可以对Microsoft SQL Server和其他企业数据库在内的大部分数据库格式创建数据库忽然前端应用程序，以及可调整的服务器端部件。利用ActiveX(TM)技术，VB可使用如Microsoft Word字处理器、Microsoft Excel电子数据表及其他Windows应用程序提供的功能，甚至可直接使用有VB专业版或企业版创建的应用程序和对象。

用户最终创建的程序是一个真正的.EXE文件，可以自由发布。

12

3.2 程序流程图

建立三维基线向量数据文本 读入GPS三维基线向量数据 根据三维基线向量观测值和固定点坐标推算待定点在WGS-84坐标系下的近似坐标 将待定点坐标转换到高斯投影平面直角坐标系中，并求出它们的坐标差，作为二维基线向量 根据地面已知点坐标及二维基线向量，将所有点坐标转到国家坐标系下

根据误差方程组法方程 解法方程，输出最终平差结果

3.３ 主要程序代码

3.3.1 读GPS基线网数据

其数据结构如下 2，105

说明 边长比例误差（单位为十分之一毫米），中央子午线经度（度） 13

y,B068,-1380207.925,5263626.847,3316898.320 说明 识别符号，固定点点名，固定点的空间坐标 B064 测站点点名

B068,-421.318,-15.016,-202.943 基线目标点点名 , 三维基线向量 ?? end

主要程序代码如下： CommonDialog1.ShowOpen

fname = CommonDialog1.FileName '将用户在\打开\对话框中选择的文件名对变量fname赋值。

If fname <> \若无此判断当对话框中选择取消时、下面赋值语句将出错。

Set ts = fso.OpenTextFile(fname) '将fname作为文本文件打开,并设置句柄。

k = 0: p = 0: h = 0: ds(j) = 0: j = 0: d = 0 'j是测站数累计变量,k是已知点累计变量,ds()为基线累积变量

Do '后测型循环，结束循环的条件是读到文件结束符\ 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, \从上一个找到的逗号位置起，查找下一个逗号的位置。

14

Loop

If 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 ,并使该句以后不能再执行。

Else

If m(1) = 0 And q <> \- 1) '该行中没有逗号,但又不是结束符，则一定是测站点名。将读出的字符串赋值到点名数组变量dm(j)。

If p > 0 And m(1) <> 0 Then '不是第一行，并且该行中有逗号分割的多个字串。

If tr(1) = \有识别符号y则这一行不是基线向量，而是固定点坐标值。

k = k + 1: ym(k) = tr(2): xo(k) = Val(tr(3)): yo(k) = Val(tr(4)): zo(k) = Val(tr(5)) '存储固定点点名及坐标值。

Else

If tr(1) <> \累计基线目标点号到点名数组变量dn(j)中 ds(j) = ds(j) + 1:dn(ds(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) = \’数据读取完毕。

3.3.2 空间近似坐标推算

（1） 累加所有点名

For i = 1 To ds(cds) '依次访问所有测站（cds）的目标方向. p = 0 '设识别变量.

For j = 1 To d '依次访问所有测站（d）.

15

If dm(j) = dn(i) Then p = 1 ' 查看目标点是否设过测站，是则对识别变量p赋值1。

If p = 0 Then d = d + 1: 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(j)’计算公式。 end if

if 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 j

Next i

Do ’后测型循环，当所有点的近似坐标算出则退出循环。

s = s + 1 's是循环计数变量，控制循环次数，避免假定坐标计算不出时，进入死循环。

For i = 1 To cds '按测站再次进入循环。

If Abs(x(i))>0 Then '在该测站假定坐标已计算出的前提下，计算这条基线另一点的假定坐标。

For j = ds(i - 1) + 1 To ds(i) '遍访i测站上所有基线。 h = seqn(dn(j)) '调用函数查目标点对应的序号。 If Abs(x(h))>0 Then

Exit For '如目标点坐标已解出，退出For－Next循环。 End if

16

Next j

For j = ds(i - 1) + 1 To ds(i) '再次遍访i测站上的所有基线。 h = seqn(dn(j)) '调用函数查目标点对应的序号。

If Abs(x(h))<0 Then '如该基线目标点坐标未求出，则用下面公式计算其假定坐标值。

x(h) = x(i) + dx(j): y(h) = y(i) + dy(j): z(h) = z(i) + dz(j) End if

Next j p = 0 ’p为识别变量。

For k = 1 To zds ’查找所有测站点。

If Abs(x(k))〈0 Then '查看是否还有没解算出坐标的点，有则进入循环再次搜索计算。

For j = ds(k - 1) + 1 To ds(k) '再次遍访k测站上所有基线 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(x(h))<0 Then p = 1 '如果目标点也未解出，则对P赋值再进入DO循环。 End if Next j Next k

Loop Until p = 0 Or s > 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) ’计算固定点真实坐标与假定的差值。

17

End If Next j Next i

For i = 1 To zds ’按站点循环，将所有点的坐标值都加上坐标改正数。 x(i) = x(i) + X1: y(i) = y(i) + Y1: z(i) = z(i) + Z1 Next i

For i = 1 To yds '置入固定坐标值

x(m(i)) = xo(i): y(m(i)) = yo(i): z(m(i)) = zo(i) Next i

3.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 Double, y As Double) '经度坐标反算通用公式的函数程序代码。

ath = Atn(y / x)

If x < 0 Then ’坐标象限的转换，pi为圆周率。 ath = ath + pi Else

If y < 0 Then ath = ath + 2 * pi End If End Function

二、 Private Function bth(x As Double, y As Double, z As Double) '纬度坐标反算公式通用函数的程序代码。

18

t=0

Do ’ 后测型循环。 t=t0

n=a/Sqr(1-e^2*sin(t0)) ’e为克拉索夫斯基椭球体的第一偏心率,a为克拉索夫斯基椭球长半轴。

t = Atn((z+n*e^2*sin(t0))/sqr(x^2+y^2))’t与t0进行叠代，直到满足精度要求。

Loop Until Abs(t-t0)<0.000000000001 ’ 叠代前后值之差的绝对值小于千分之一秒时退出。

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 + e^2 * b2 ^ 2)) 'n为该站点所在卯酉圈的曲率半径，c 为克拉索夫斯基椭球体极点处的子午线曲率半径,e为克拉索夫斯基椭球体的第二偏心率。

x0 = 111134.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) / 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)

19

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)) '调用函数查目标点对应的序号。

kx(j) = gx(k) - gx(i): ky(j) = gy(k) - gy(i) ’计算公式。 Next j

Private Function seqn(str As String) As Integer '由点名查计算序号函数

For i = 1 To zds

If str = dm(i) Then seqn = i '将查到的序号赋给函数名,返回调用处 Next i End Function

3.3.4 平面近似坐标计算的程序代码

这部分的程序与前面空间近似坐标的计算相同，首先读入平面已知点坐标，格式为：点名，X坐标，Y坐标

?? end

接着进行平面近似坐标的推算，推算程序和前面相同。

3.3.5 组法方程的程序代码

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 '按测站循环

20

下表2是用本程序平差的结果 测站点名 x坐标 y坐标 点位误差mm 0.0 7.0 13.6 10.4 12.7 10.1 0.0 6.1 13.2 10.1 11.7 9.2 9.0 Vx(mm) 0.00 11.76 10.27 11.54 11.60 0.68 0.00 1.57 6.85 7.71 8.34 6.19 6.33 Vy(mm) 0.00 -11.00 -14.44 12.22 -13.95 1.15 0.00 -0.89 -4.95 -9.80 -11.57 -6.22 -7.71 B064 26060.725 16974.147 B068 25838.305 17384.832 B093 26428.482 18060.018 B106 25928.200 17654.482 B109 26038.260 17972.342 B115 26313.755 16431.216 B127 26652.060 16582.485 B129 26456.555 16619.608 N050 26275.994 16617.868 B085 26338.220 17558.521 B089 26459.481 17753.973 B077 26104.054 17204.929 B078 26323.572 17359.419 对两表的点位误差进行比较，B093点的差值11.4mm和N050的差值11.7mm最大，这两点为最弱点。平差时，所取已知点的数目不同也会影响到平差的结果。

26

结 论

本人在设计工作中，深入理解了GPS数据后处理、测量平差基础的理论与方法，提高了运用VB语言设计复杂数值计算程序的能力，在此过程中，深刻体会到：

（1）在基线转换过程中应仔细检查及认真调试高斯投影正算公式，直至计算结果准确无误。

（2）组法方程时，应考虑所有可能情况，避免少组或漏组法方程，未知数序号的计算也要准确无误，否则都有可能导致程序无法运行或平差结果有误。

（3）对于程序设计过程中遇到的各种问题，要细心检查，逐步调试，直到找出问题的根源，才能解决问题完善程序。

由于时间仓促及个人能力有限，本设计所完成的GPS数据处理程序还不够完善，没有做高程拟合、与地面观测值联合平差等内容，但通过设计，我已经掌握了相关的理论和方法，初步具备了数值计算从设计能力，为今后进一步的摄入学习、研究奠定了基础。

27

参考文献

[1] 孔祥元 梅是义编. 控制测量学[M]. 武汉： 武汉大学出版社，2002 [2] 刘基余 李征航 王跃虎 桑吉章编. 全球定位系统原理及其应用[M]. 北京：

测绘出版社，1999

[3] 郑阿奇编. Visual Basic 实用教程[M]. 北京：电子工业出版社. 2000 [4] 李玉宝编. 测量平差程序讲义. 绵阳：西南科技大学测绘教研室. 1996 [5] 马明栋 王之勋 沈 蔚编. ＧＰＳ基线向量网平差模型与程序设计[J]. 辽

宁： 辽宁工程技术大学学报第3期，2001

[6] 施一民编. 工程ＧＰＳ控制网平差转换的要点与模型[J]. 上海:同济大学

出版社. 2003

[7] 王晓海编. ＧＰＳ及应用新发展[J]. 西安: 航天科技集团总公司五院五零

四研究所.2004

[8] 赵长胜 金继读 王忠义编. GPS基线网数据处理系统的设计和实现[J]. 辽

宁： 辽宁工程技术大学学报第6期，2001

[9] 武汉大学测量平差教研室编. 测量平差基础[M]. 北京：测绘出版社，1996 [10] T. P. Yunck. The Anatomy of Global Positioning System[J]. GPS World,

Vol.3, No.5, 1992

[11] A. Leick, GPS Satellite Surveying[J]. John Wiley and Sons, 1990

28

For j = ds(i - 1) + 1 To ds(i) '在i测站上按方向循环，求误差方程系数、常数。

ReDim nb(n1) ' 重新定义误差方程系数数组,并且每循环到一新边长前清零。

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 = \测站点i不是已知点

d = 2 * (i - k) – 1 '计算测站i点x未知数在未知数点集中的序号。 chu = charact(h, s) '自定义函数，查目标点j是否已知点，如不是,用s返回序号h前面有几个已知点。

If chu =\’ 如果目标点不是已知点，则组法方程。

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 - 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=\then ’如果目标点为已知点则以nb(h) 和nb(h+1)都为零组法方程。

nb(d) = -1: nb(n1 - 1) = kx(j): nb(n1) = -ky(j): l1 = x(h)-x(i) -

21

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= =\’如果i测站为已知点时组则以nb(d)和nb(d+1)为零组法方程。

chu = charact(h, s) '自定义函数，先查目标点j是否已知点. If chu=\’如不是,用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(), pp, l1) '调用函数组法方程，同上。 End If ’这种情况判断完毕，进行下一种判断。

chu = charact(h, s) '自定义函数，先查目标点j是否已知点，如不是,用s返回序号h前面有几个已知点。

If chu=\’如果目标点也是已知点，则退出循环。 Exit for End if End if

Next j '一个目标点的误差方程处理完毕，进入下一个目标方向。 ll = ll + pp * l1 ^ 2 '累积和方程常数项。

Next i '一个测站的误差方程处理完毕，进入下一测站

22

4 实例平差对比与结果分析

4.1 实例数据准备

先将南方GPS adj4.0处理所得的在WGS84-坐标系下经过自由网平差的三维基线向量减去相应的改正数，作为原始数据。做法如下：

基 线 名 △X △Y △Z △X改正 △Y改正 △Z改正 B064，B068 -421.317 -15.016 -202.943 0.987 0.133 0.068 B077，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,105

y,B064,-1379786.609,5263641.863,3317101.263 B077

B078,-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 B085

B078,188.040,67.722,-6.907 B093,-472.221,-175.392,76.949 B109,-440.175,43.702,-256.151 B089

23

B109,-267.982,156.320,-359.022 B078,360.235,180.330,-109.788 B085,172.189,112.614,-102.875 B093,-300.032,-62.787,-25.939 B093

B068,571.449,476.216,-500.935 B109,32.050,219.101,-333.087 B106

N050,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 B109

B068,539.390,257.101,-167.862 B078

B129,724.238,156.416,133.731 B127,783.968,77.167,306.665 B115

B127,-105.592,-189.944,301.135 B129

B127,59.729,-79.248,172.937 N050

B068,-792.580,-6.202,-391.356 B064

B068,-421.318,-15.016,-202.943 B106,-668.305,-135.230,-129.509 B115,557.066,18.996,219.687 B127,451.477,-170.947,520.816 N050,371.264,-8.787,188.415

24

end （“end” 为结束识符）

地面已知点数据如下：

B064, 26060.725, 16974.147 B127, 26652.060, 16582.485 end

4.2 平差成果对比及分析

下表1是用南方测绘GPS数据处理软件平差的结果 点名 B064 X坐标 Y坐标 rms(mm) 0.000 1.582 2.220 1.478 2.162 2.672 0.000 2.835 1.543 2.435 2.494 2.607 2.587 Vx(mm) 0.000 1.181 1.579 1.129 1.528 2.048 0.000 2.110 1.162 1.771 1.830 1.935 1.908 Vy(mm) 0.000 1.053 1.561 0.953 1.530 1.716 0.000 1.893 1.014 1.672 1.693 1.747 1.748 26060.72 16974.147 B068 25838.323 17384.844 B093 26428.474 18060.000 B106 25928.201 17654.479 B109 26038.258 17972.332 B115 26313.759 16431.227 B127 26652.060 16582.485 B129 26456.558 16619.611 N050 26275.998 16617.882 B085 26338.220 17558.501 B089 26459.479 17753.952 B077 26104.053 17204.910 B078 26323.577 17359.403

25

本文来源：https://www.bwwdw.com/article/z6df.html