1、VB测量平差程序的设计讲稿Case 0 读入观测值文件Text1.Visible = FalseCommonDialog1.ShowOpenfname = CommonDialog1.FileName 将用户在打开对话框中选择的文件名对变量fname赋值If fname Then 若无此判断当对话框中选择取消时、下面赋值语句将出错 Set ts = fso.OpenTextFile(fname) 将fname作为文本文件打开,并设置句柄 j = 0: k = 0: p = 0: h = 0 j是测站数累计变量,k是已知点累计变量,l（j）、ns（j）分别是方向值、边长累积计数 Do While
2、 ts.AtEndOfLine True 前测型循环，进入循环的条件是没有读到文件结束尾 B = ts.ReadLine 读一行，置入b B = Trim(B): i = 1: 删除B可能有的前导和尾随空格，i是工作变量， m(i) = InStr(B, ,) 查行中第一个逗号的左数位置,并保存在整形数组变量m(i) Do While m(i) 0 前测型Do. Loop循环，成立条件是该行字符串中有逗号 tr(i) = Mid(B, m(i - 1) + 1, m(i) - m(i - 1) - 1) 提取指定位置开始的指定数目字符。 i = i + 1 m(i) = InStr(m(i -
3、 1) + 1, B, ,) 从上一个找到的逗号位置起，查找下一个逗号的位置 Loop If m(i) = 0 And i 1 Then tr(i) = Right(B, Len(B) - m(i - 1) 处理一行中最后一个逗号后的字符串 以下部分是将存储在数组变量 m(i)中的字符分类存放到方向、边长、已知坐标、网型信息等数组中 If p = 0 Then 读到的是文件第一行。 ma = tr(1): ms = tr(2): mk = tr(3): p = 1 提取观测方向,边先验精度值，并使该句以后不能再执行。 Else If m(1) = 0 Then j = j + 1: ReDim
4、 Preserve dm(j): ReDim Preserve nl(j): ReDim Preserve ns(j): dm(j) = B: nl(j) = nl(j - 1): ns(j) = ns(j - 1) 该行中没有逗号,但又不是结束符，则一定是测站点名。将读出的字符串赋值到点名数组变量dm(j)， Next j If p = 0 Then d = d + 1: ReDim Preserve ns(d): ns(d) = ns(d - 1): ReDim Preserve dm(d): dm(d) = lb(i) 如 p=0，表明目标点未设过测站,将该点点名加入点名数组 Next
5、i zds = d 将总点数存入模块级变量zds ReDim x(zds), y(zds) 重新定义坐标数组 x(1) = 10000: y(1) = 10000 为推算近似坐标,对第一个点赋假设坐标值 k = 1 For i = 1 To yds If lb(1) = ym(i) Then k = k + 1 Next i ss = sid(1, k) 调出第一点到未知点方向的边长，参数是测站点序号，照准方向号 h = seqn(lb(k) 查k方向照准点的计算序号 x(h) = x(1) + ss * Cos(0): y(h) = y(1) + ss * Sin(0) 计算第一点上第k方向
6、值的目标点假设坐标 For i = 1 To nl(cds) 遍访所有方向值，将其由角度值转换为弧度值. If l(i) 0.001 Then l(i) = radian(l(i) 零方向值不参加转换 Next i n = 0 Do n = n + 1 n是循环计数变量，控制循环次数，避免假定坐标计算不出时，进入死循环。 For i = 1 To cds 按测站循环 If x(i) 1 Then 在该测站假设坐标已计算出的前提下，求照准点假设坐标 p1 = 0 For j = nl(i - 1) + 1 To nl(i) 遍访i测站上所有方向值 h = seqn(lb(j) 查目标点对应的序号
7、 If x(h) 1 Then 目标点坐标已解出 t2 = azimuth(xo(1), yo(1), xo(2), yo(2) dt = t2 - t1: x1 = x(m(1): y1 = y(m(1) For i = 1 To zds 将假设坐标转换到实际坐标 dx = x(i) - x1: dy = y(i) - y1 x(i) = xo(1) + dx * Cos(dt) - dy * Sin(dt) y(i) = yo(1) + dx * Sin(dt) + dy * Cos(dt) Next i For i = 1 To yds 置入已知点坐标 x(m(i) = xo(i): y
8、(m(i) = yo(i) Next iCase 2 组法方程Text1.Visible = FalseDim l1 As Double, pp As Double, n2 As Long 定义过程级变量q = 206265: ll = 0n1 = 2 * (zds - yds) 未知数数目n2 = n1 * (n1 + 1) / 2 一维存储法方程系数数组上限ReDim NX(n2), UX(n1) 重新定义法方程系数、常数数组Call order(m(), yds) 对保存已知点序号的m()数组排序For i = 1 To cds 按测站循环 z = 0 将按测站累积的变量清零 下面开始处
9、理一个测站的方向观测值 k1 = nl(i - 1) + 1: k2 = nl(i) 一测站上最小和最大方向号 For j = k1 To k2 在i测站上按方向循环，求定向角未知数 h = seqn(lb(j) t = azimuth(x(i), y(i), x(h), y(h) f = t - l(j): If f = 0 Then i测站有观测边 For j = ns(i - 1) + 1 To ns(i) 依次遍访i测站上各观测边 ReDim nb(n1) 重新定义误差方程系数数组,并且每循环到一新边长前清零 h = seqn(sb(j) t = azimuth(x(i), y(i),
10、 x(h), y(h) A = Cos(t): B = Sin(t) cha = charact(i, k) 自定义函数，查测站点i是否已知点，如不是,用k返回i前面有几个已知点 If cha = n Then 测站点i不是已知点 d = 2 * (i - k - 1) + 1 计算测站i点x坐标未知数在未知数点集中的序号 nb(d) = -A: nb(d + 1) = -B End If cha = charact(h, k) If cha = n Then 照准方向点h不是已知点 d = 2 * (h - k - 1) + 1 计算照准方向h点x坐标未知数在未知数点集中的序号 nb(d)
11、= A: nb(d + 1) = B End If dx = x(h) - x(i): dy = y(h) - y(i) ss = Sqr(dx 2 + dy 2) 反算边长，置于判断式外是因为两已知点之间不会测边 pp = (ma / (0.1 * ms + mk * ss * 10 -4) 2 边长观测值定权 l1 = (ss - s(j) * 100 求边误差方程常数(单位是厘米) 并对该测站最大方向号，最大边长号数组变量nl(j)、ns(j)赋值累计起始值 If m(1) 0 Then 不是第一行，并且该行中有逗号分割的多个字串 tr(2) = LCase(tr(2) If tr(2)
12、 l And tr(2) s Then 这一行不是方向或边长观测值，而是已知点坐标值 k = k + 1: ReDim Preserve ym(k): ReDim Preserve xo(k): ReDim Preserve yo(k): ym(k) = tr(1): xo(k) = Val(tr(2): yo(k) = Val(tr(3) 输入已知点坐标 Else If tr(2) = l Then nl(j) = nl(j) + 1: ReDim Preserve lb(nl(j): ReDim Preserve l(nl(j): lb(nl(j) = tr(1): l(nl(j) = V
13、al(tr(3) 累计测站方向数、提取照准点、方向值 If tr(2) = s Then ns(j) = ns(j) + 1: ReDim Preserve sb(ns(j): ReDim Preserve s(ns(j): sb(ns(j) = tr(1): s(ns(j) = Val(tr(3) 提取边观测数、提取照准点、观测边 End If End If End If Loop ts.Close cds = j: yds = k 用模块级变量cds 、 yds保存测站点总数、已知点总数 MsgBox 数据已成功读入, 0 + 64 + 0, 数据输入End IfCase 1 解算近似坐标
14、 d = cds d是测站数 For i = 1 To nl(cds) 依次访问所有的方向值 p = 0 设识别变量 For j = 1 To d 依次访问所有测站 If dm(j) = lb(i) Then p = 1 查看目标点是否设过测站，是则对识别变量p赋值1。 t = azimuth(x(i), y(i), x(h), y(h): p1 = j 反算方位角并记下方向号 Exit For 退出ForNext循环 End If Next j For j = nl(i - 1) + 1 To nl(i) 再次遍访i测站上的所有方向值 h = seqn(lb(j) 查目标点对应的序号 If
15、x(h) 1 Then 该照准方向点坐标未求出 tt = t + l(j) - l(p1) 计算j方向的方位角 ss = sid(i, j) 根据测站号、方向号去查找边长 If ss = 0 Then ss = side(i, h) x(h) = x(i) + ss * Cos(tt): y(h) = y(i) + ss * Sin(tt) 计算近似坐标 End If Next j End If Next i p = 0 For k = 1 To zds If x(k) zds 坐标已全部算出或虽还有未算出的，但根据循环次数已不能算出时结束循环。 For i = 1 To yds 查点名数组中
16、那些点是已知点,用m()数组存其序号 For j = 1 To zds If ym(i) = dm(j) Then m(i) = j Next j Next i t1 = azimuth(x(m(1), y(m(1), x(m(2), y(m(2) 利用一对公共点求旋转角 z = z + f Next j zo = z / (k2 - k1 + 1) zo是定向角未知数（零方向的坐标方位角） ReDim nc(n1): ln = 0 重新定义和方程系数数组，并且每循环到一新测站前清零，ln是和方程常数项 For j = k1 To k2 在i测站上按方向循环，求误差方程系数、常数 ReDim
17、nb(n1) 重新定义误差方程系数数组,并且每循环到一新方向前清零 h = seqn(lb(j) dx = x(h) - x(i): dy = y(h) - y(i) ss = Sqr(dx 2 + dy 2) 反算边长 t = azimuth(x(i), y(i), x(h), y(h) 反算坐标方位角 A = q * Sin(t) * 10 -2 / ss: B = -q * Cos(t) * 10 -2 / ss 求方向误差方程系数 cha = charact(i, k) 自定义函数，查测站点i是否已知点，如不是,用k返回i前面有几个已知点 If cha = n Then 序号为i的测站
18、点不是已知点 d = 2 * (i - k - 1) + 1 计算测站i点x坐标未知数在未知数点集中的序号 nb(d) = A: nb(d + 1) = B 对误差方程系数数组赋值 nc(d) = nc(d) + nb(d): nc(d + 1) = nc(d + 1) + nb(d + 1) 累积到和方程数组中去 End If cha = charact(h, k) 查照准方向点h是否已知点，如不是,用k返回h前面有几个已知点 If cha = n Then 照准方向点不是已知点 d = 2 * (h - k - 1) + 1 计算照准方向点x坐标未知数在未知数点集中的序号 nb(d) =
19、-A: nb(d + 1) = -B 对误差方程系数数组赋值 nc(d) = nc(d) + nb(d): nc(d + 1) = nc(d + 1) + nb(d + 1) 累积到和方程数组中去 End If If t - zo -0.01 Then t = t + 2 * pi 当t是零方向方位角时，可能比平均值zo略小，所以不用0的比较方法 l1 = q * (t - zo - l(j) 误差方程常数项,以秒为单位 ln = ln + l1: ll = ll + l1 2 累积和方程常数项及求pvv的Pll值Call equation(nb(), pp, l1) ll = ll + pp
20、 * l1 2 累积用于求pvv的Pll值 Next j End If一个测站的误差方程处理完毕，进入下一测站Next iText1.Visible = FalseMsgBox (法方程组成完毕)Case 3 高斯约化法解法方程Dim m2 As Longn = 2 * (zds - yds)ReDim nb(n), nc(n) 数组清零For i = 1 To n nb(i) = -UX(i) 将法方程常数项赋予一工作数组,保留原值用于求pvvNext iFor k = 1 To n - 1 从法方程的第一行到倒数第二行, k实际上控制约化次数 m1 = k * (k + 1) / 2 m1
21、是k行的自乘元素在法方程系数阵中序号（一维按列上三角阵） For j = k + 1 To n 在j循环内完成k次约化 m2 = (j - 1) * j / 2) + k m2在j循环内依次是k行自乘系数右边各元素序号 For i = j To n i循环内完成法方程系数阵j行元素的k次约化 d = (i - 1) / 2 * i n1 = d + k 与约化元素同列的 k行元素 n2 = d + j 约化元素在法方程系数阵中编号 NX(n2) = NX(n2) - NX(n1) * NX(m2) / NX(m1) Next i Next j For i = k + 1 To n 在i循环内完
22、成法方程常数项的k次约化 nub = (i - 1) * i / 2 + k nb(i) = nb(i) - nb(k) * NX(nub) / NX(m1) Next iNext k 法方程约化完毕pvv1 = ll ll是pllFor k = 1 To n 在循环内按pvv=pll+(l)*(l)求pvv nub = k * (k + 1) / 2pvv1 = pvv1 - nb(k) * nb(k) / NX(nub)Next kFor k = n To 1 Step -1 在k循环内回代求解未知数 u = nb(k) If k n Then For i = k + 1 To n nub
23、 = (i - 1) * i / 2 + k u = u - NX(nub) * nb(i) Next i nb(k) = u End If nub = k * (k + 1) / 2 nb(k) = nb(k) / NX(nub)Next kpvv2 = llFor k = 1 To n 按pvv=pll+W*X求pvv pvv2 = pvv2 + UX(k) * nb(k)MsgBox (法方程解算完毕)End SelectEnd SubPrivate Function azimuth(x1 As Double, y1 As Double, x2 As Double, y2 As Doub
24、le) 反算坐标方位角dx = x2 - x1: dy = y2 - y1azimuth = Atn(dy / dx)If dx 0 Then azimuth = azimuth + piElse If dy = 0 Then 该测站有观测边，确定m方向的边长 For k = ns(n - 1) + 1 To ns(n) If sb(k) = lb(m) Then ss = s(k) 依次将边的照准点名与m方向照准点名对比，相等则赋予ss。 Next kEnd IfIf ns(n) - ns(n - 1) - 1 0 Or ss 1 Then 如测站n未找到m方向的边，则转到lb(m)站去找
25、h = seqn(lb(m) 由点名查计算编号 NX(h) = NX(h) + B(j) * B(i) * p 组法方程系数阵 End If Next j End IfNext iEnd SubPrivate Sub order(m() As Integer, n As Integer) 排序通用过程For i = 1 To n - 1 For j = i + 1 To n If m(j) m(i) Then A = m(j): m(j) = m(i): m(i) = A End If Next jNext iEnd SubPrivate Function side(n As Long, m
26、As Long) 反算已知点间边长函数过程 For i = 1 To yds If ym(i) = dm(n) Then k1 = i If ym(i) = dm(m) Then k2 = i Next i dx = xo(k1) - xo(k2): dy = yo(k1) - yo(k2) side = Sqr(dx 2 + dy 2)End FunctionPrivate Sub inversion(A() As Double, W() As Double) 法方程系数阵求逆子程序n = 2 * (zds - yds) 形参数组A()是约化后法方程系数数组,W()是工作数组。For i =
27、 n To 1 Step -1 i控制所求逆阵元素的行数，由下而上。 End IfNext iEnd Sub 求逆完毕,此时约化法方程系数阵中元素已是逆阵元素8888888777770000999900000000009999Next kFor i = 1 To zds 求坐标平差值 cha = charact(i, k) 自定义函数，查点i是否已知点，如不是,用k返回i前面有几个已知点 If cha = n Then 序号i的点不是已知点 h = i - k d = 2 * (h - 1) + 1 计算i点x坐标未知数在未知数点集中的序号 x(i) = x(i) + nb(d) / 100:
28、 y(i) = y(i) + nb(d + 1) / 100 Else End IfNext ima = Sqr(pvv1 / (nl(cds) + ns(cds) - n - cds) 求单位权中误差Call inversion(NX(), nb() 将约化后的法方程系数阵送入通用过程inversion()求逆ReDim nb(n) 工作数组清零For i = 1 To zds 按点号循环,求点位中误差 cha = charact(i, k) 自定义函数，查点i是否已知点，如不是,用k返回i前面有几个已知点 If cha = n Then 不是已知点 n1 = 2 * (i - k - 1)
