daoxianwangpingcha.txt

这是精密导线平差的VB源代码

'***********************************************************************************
'    本程序用于单一附（闭）合导线严密平差计算，采用按角度条件平差法。以方向观测值中
'误差的先验值作为单位权中误差。计算结果可求得各待定点的坐标平差值及其点位精度Mx，My，
'及M，并计算出各待定点误差椭圆元素E，F，Z
'
'参考文献： 郭久训．《控制网平差程序设计》　北京：原子能出版社，2004．8
'
'平差数据来源：潘正风等．《数字测图原理与方法》　武汉大学出版社　　186页表6-5
'
'等级：所用平差数据为首级图根导线数据（精度很不高）。本程序中写出了方位角和导线全长相对闭合差的判
'别，但考虑到程序的通用性，将这些限差判断当作了注释处理，而不实际运行。
'
'程序不足：没有导线网的图形表达。
'***********************************************************************************

Private i%, j%, n%, tc#, tb#, B_x!, B_y!, e1!, e2!, m!, m0#, z#, aa As Boolean, bb As Boolean, cc As Boolean ',dd As Boolean
Private Naa#(2, 2), Naa逆#(2, 2), W#(2), K#(2), qq#(2), fx#(2), fy#(2)
Private A#(), Q#(), V#(), C#(), mx#(), my#(), mk#(), e#(), f#(), zz#()
Private Po() As Point





'文件格式说明：
'文件格式详见文件"平差数据.txt"
'
Private Sub 打开文件_Click()                  '打开文件
    Dim ff$, temp$, A_name$, A_x!, A_y!, A_l#, A_s!, B_name$, B_l#, C_name$, C_x!, C_y!, D_name$, D_x!, D_y!
    Form1.Cls  '清屏
    CommonDialog1.DialogTitle = "打开数据文件"
    CommonDialog1.FileName = ""
    CommonDialog1.ShowOpen
    '出错处理
    On Error GoTo FileErr
        ff = CommonDialog1.FileName           'ff是文件路径名
        Open ff For Input As #1               '以顺序文件方式打开文件，使用input
        
        Line Input #1, temp                   '读取文件中的说明语句
        Line Input #1, temp                   '读取文件中的说明语句
        Line Input #1, temp                   '读取文件中的说明语句
        Input #1, n                           '读取n，n为 （测站数-1）
        ReDim Po(n) As Point                  '定义Po(n)，其中Po(0)存A点数据，Po(n)存B点数据，Po(1)到Po(n-1)存n-1个未知点数据。
        Line Input #1, temp                   '读取文件中的说明语句
        Line Input #1, temp                   '读取文件中的说明语句
        Input #1, m, e1, e2                   '读取先验方向观测值中误差m，测距仪固定误差e1，比例误差e2
        Line Input #1, temp                   '读取文件中的说明语句
        Line Input #1, temp                   '读取文件中的说明语句
        Input #1, C_name, C_x, C_y            '读取已知点C
        Input #1, A_name, A_x, A_y, A_l, A_s  '读取已知点A
        
        For i = 1 To n - 1                    '读取n-1个未知点
            Input #1, Po(i).name, Po(i).l, Po(i).s
        Next i
        
        Input #1, B_name, B_x, B_y, B_l       '读取已知点B
        Input #1, D_name, D_x, D_y            '读取已知点D
                      
        Close #1
        
        tc = ZBiaoFSuan(C_x, C_y, A_x, A_y)    '坐标反算,求点C到A的坐标方位角，并记作tc，单位是度
        tb = ZBiaoFSuan(B_x, B_y, D_x, D_y)    '坐标反算,求点B到D的坐标方位角，并记作tb，单位是度
        
        Po(0).name = A_name: Po(0).x = A_x: Po(0).y = A_y: Po(0).l = A_l: Po(0).s = A_s  'Po(0)存A点数据
        Po(n).name = B_name: Po(n).x = B_x: Po(n).y = B_y: Po(n).l = B_l                 'Po(n)存B点数据
        For i = 0 To n
            Po(i).l = deg(Po(i).l)             '将观测方向左角的单位度分秒化作度
        Next i
        
        ReDim Q#(2 * n), A#(2, 2 * n), V#(2 * n), C#(1, 2 * n - 2), mx#(n - 1), my#(n - 1), _
        mk#(n - 1), e#(n - 1), f#(n - 1), zz#(n - 1)        '变量重新定义
        bb = True
        MsgBox "文件已成功打开", , "提示"
        显示平差数据 ff
        Exit Sub
FileErr:
        MsgBox "您的文件未打开或打开的文件格式有误！注意：请重新运行本程序！！", , "提示"
End Sub

'
'
'
Private Sub 开始计算_Click()
    
    If bb = True Then
        推算方位角和坐标         '求近似坐标
        求条件式
        'If dd = True Then        '若平差数据没有超限，则进行下面的计算
        组法方程式
        求逆求K
        求改正数和平差值
        推算方位角和坐标         '求平差后坐标
        精度评定
        'Else
            'MsgBox "限差超限！", , "提示"
            'Exit Sub
        'End If
    Else
        MsgBox "数据文件未打开！！", , "提示"
        Exit Sub
    End If
    aa = True
    MsgBox "计算完毕！", , "提示"
    显示平差结果
End Sub

'
'
'
Private Sub 保存_Click()
    Dim ff$
    
    If aa = True Then
        CommonDialog1.DialogTitle = "保存平差结果"
        CommonDialog1.Filter = "(*.txt)|*.txt"
        CommonDialog1.FileName = "平差结果.txt"
        CommonDialog1.ShowSave
        ff = CommonDialog1.FileName
        If ff = "" Then
            MsgBox "文件名不能为空", , "警告"
            Exit Sub
        End If
          
        Open ff For Output As #2                    '以顺序文件方式保存文件，使用output
            
            Print #2, "平差结果："
            Print #2, "---------------------------------------------------------------------------------------------"
            Print #2, " 单位权中误差 m0=" & Format(m0, "####.##") & " 秒"
            Print #2, "---------------------------------------------------------------------------------------------"
            Print #2, "点名", "x(m)", "y(m)", "点位误差(cm)", "E(cm)", "F(cm)", "Z(度分秒)"
            For i = 1 To n - 1
                Print #2, " " & Po(i).name, Format(Po(i).x, "########.####"), Format(Po(i).y, "########.####"), "     " & Format(mk(i), "#####.##"), _
                        "      " & Format(e(i), "#####.##"), "      " & Format(f(i), "#####.##"), "      " & Format(zz(i), "###.#####")
            Next i
            Print #2, "---------------------------------------------------------------------------------------------"
        Close #2
          
        cc = True
        MsgBox "保存完毕！", , "提示"
    Else
        MsgBox "没有需要保存平差结果！", , "提示"
        Exit Sub
    End If
End Sub

'
'
'
Private Sub 退出_Click()
    If aa = True And cc = False Then
            提示.Show                         '若平差结果未保存则进行提示
    Else
            End                               '若没有平差结果或平差结果已保存则退出
    End If
End Sub

'
'
'
Private Sub 推算方位角和坐标()
    Po(0).t = tc + Po(0).l - 180
    If Po(0).t > 360 Then Po(0).t = Po(0).t - 360
    If Po(0).t < 0 Then Po(0).t = Po(0).t + 360
    
    For i = 1 To n
        '求方位角
        Po(i).t = Po(i - 1).t + Po(i).l - 180
        If Po(i).t > 360 Then Po(i).t = Po(i).t - 360
        If Po(i).t < 0 Then Po(i).t = Po(i).t + 360
        '求坐标
        Po(i).x = Po(i - 1).x + Po(i - 1).s * Cos(Po(i - 1).t / p0)   't要以度为单位
        Po(i).y = Po(i - 1).y + Po(i - 1).s * Sin(Po(i - 1).t / p0)
    Next i
End Sub

'
'角度改正数的单位是秒，边长改正数的单位是厘米
'方向观测值中误差为单位权
'方向观测值中误差的先验值为 m(秒)
'角度观测值的权为 0.5
'距离观测值中误差的先验值为 ms=e1+e2*s*0.0001(厘米)
'
Private Sub 求条件式()
    Dim sums!
    
    '求A
    For i = 0 To n       '求角度改正数系数
        A(0, n + i) = 1
        A(1, n + i) = (Po(i).y - Po(n).y) / 2062.65
        A(2, n + i) = (Po(n).x - Po(i).x) / 2062.65
        Q(n + i) = 2     '求角度观测值的权倒数
    Next i
    For i = 0 To n - 1   '求边长改正数系数
        A(1, i) = Cos(Po(i).t / p0)
        A(2, i) = Sin(Po(i).t / p0)
        Q(i) = (e1 + e2 * Po(i).s * 0.0001) ^ 2 / m ^ 2 '求距离观测值的权倒数
    Next i
    '求W
    W(0) = (Po(n).t - tb) * 3600  '单位是s
    W(1) = (Po(n).x - B_x) * 100  '单位是cm
    W(2) = (Po(n).y - B_y) * 100  '单位是cm
    For i = 0 To n - 1
        sums = sums + Po(i).s
    Next i
   
    '限差判断
    '这是首级图根导线的精度要求，其中 n+1 为测站数
    'If Abs(W(0)) > 40 * Sqr(n + 1) Or Sqr((W(1) / 100) ^ 2 + (W(2) / 100) ^ 2) / sums > 0.00025 Then '方位角闭合差和导线全长相对闭合差
        'Exit Sub
    'End If
    'dd = True
End Sub

'
'
'
Private Sub 组法方程式()
    Dim g%
    For i = 0 To 2
        For j = i To 2        '由于Naa是对称的，只求其上三角即可。
            For g = 0 To 2 * n
                Naa(i, j) = Naa(i, j) + A(i, g) * A(j, g) * Q(g)
            Next g
        Next j
    Next i
End Sub

'
'
'对于3*3的矩阵来说，利用公式　Naa逆=Naa*/|Naa| 　去求Naa的逆即可，其中Naa*是伴随矩阵。
Private Sub 求逆求K()
    Dim ee#, ii#, jj#, kk#, det#       'det是|Naa|的值
    ee = 1
    For i = 0 To 2
        ee = ee * Naa(i, i)
    Next i
    ee = ee + 2 * Naa(0, 1) * Naa(1, 2) * Naa(0, 2)
    ii = Naa(0, 1) ^ 2
    jj = Naa(0, 2) ^ 2
    kk = Naa(1, 2) ^ 2
    det = ee - (ii * Naa(2, 2) + jj * Naa(1, 1) + kk * Naa(0, 0))
    '求伴随矩阵Naa*，用Naa逆存储
    Naa逆(0, 0) = Naa(1, 1) * Naa(2, 2) - kk
    Naa逆(0, 1) = Naa(0, 2) * Naa(1, 2) - Naa(0, 1) * Naa(2, 2)
    Naa逆(1, 0) = Naa逆(0, 1)
    Naa逆(1, 1) = Naa(0, 0) * Naa(2, 2) - jj
    Naa逆(0, 2) = Naa(0, 1) * Naa(1, 2) - Naa(0, 2) * Naa(1, 1)
    Naa逆(2, 0) = Naa逆(0, 2)
    Naa逆(2, 2) = Naa(0, 0) * Naa(1, 1) - ii
    Naa逆(1, 2) = Naa(0, 2) * Naa(0, 1) - Naa(0, 0) * Naa(1, 2)
    Naa逆(2, 1) = Naa逆(1, 2)
    '求Naa逆
    For i = 0 To 2
        For j = 0 To 2
            Naa逆(i, j) = Naa逆(i, j) / det
        Next j
    Next i
    '求K
    For i = 0 To 2
        For j = 0 To 2
            K(i) = K(i) - Naa逆(i, j) * W(j)
        Next j
    Next i
End Sub

'
'
'
'
Private Sub 求改正数和平差值()
    Dim pvv#
    pvv = 0
    For i = 0 To 2 * n
        For j = 0 To 2
            V(i) = V(i) + K(j) * A(j, i) * Q(i)     '求改正数V
        Next j
        pvv = pvv + V(i) * V(i) / Q(i)              '求 pvv
    Next i
    m0 = Sqr(pvv / 3)                               '求单位权中误差 m0
    For i = 0 To n - 1
        Po(i).s = Po(i).s + V(i) / 100              '求导线边长的平差值，其中边长改正数的单位是厘米，边长的单位是米
    Next i
    For i = 0 To n
        Po(i).l = Po(i).l + V(n + i) / 3600         '求导线左角的平差值，其中角度改正数的单位是秒，角度的单位是度
    Next i
End Sub

'
'列出待定点坐标的权函数式，据此进行精度评定。
'
Private Sub 精度评定()
    Dim uu#, vv#, ww#, ss#, r#, kk!
    For j = 1 To n - 1
        uu = 0: vv = 0: ww = 0
        For i = 0 To j - 1
            C(0, i) = Cos(Po(i).t / p0)
            C(1, i) = Sin(Po(i).t / p0)
            C(0, n + i) = (Po(i).y - Po(j).y) / 2062.65
            C(1, n + i) = (Po(j).x - Po(i).x) / 2062.65
            uu = uu + C(0, i) ^ 2 * Q(i) + C(0, n + i) ^ 2 * Q(n + i)
            vv = vv + C(1, i) ^ 2 * Q(i) + C(1, n + i) ^ 2 * Q(n + i)
            ww = ww + C(0, i) * C(1, i) * Q(i) + C(0, n + i) * C(1, n + i) * Q(n + i)
        Next i
        For kk = 0 To 2
            ss = 0: r = 0
            For i = 0 To j - 1
                ss = ss + A(kk, i) * C(0, i) * Q(i) + A(kk, n + i) * C(0, n + i) * Q(n + i)
                r = r + A(kk, i) * C(1, i) * Q(i) + A(kk, i + n) * C(1, n + i) * Q(n + i)
            Next i
            fx(kk) = ss: fy(kk) = r
        Next kk
        For kk = 0 To 2
            ss = 0: r = 0
            For i = 0 To 2
                ss = ss - Naa逆(kk, i) * fx(i)
                r = r - Naa逆(kk, i) * fy(i)
            Next i
            K(kk) = ss: qq(kk) = r
        Next kk
        For i = 0 To 2
            uu = uu + fx(i) * K(i)
            vv = vv + fy(i) * qq(i)
            ww = ww + fx(i) * qq(i)
        Next i
        mx(j) = m0 * Sqr(Abs(uu))
        my(j) = m0 * Sqr(Abs(vv))
        mk(j) = Sqr((mx(j)) ^ 2 + (my(j)) ^ 2)        '每个点的点位中误差，单位cm
        r = uu + vv
        ss = Sqr((uu - vv) ^ 2 + 4 * ww ^ 2)
        e(j) = m0 * Sqr((r + ss) / 2)                 '点位误差椭圆元素E，单位cm
        f(j) = m0 * Sqr(Abs(r - ss) / 2)              '点位误差椭圆元素F，单位cm
        z = Atn(2 * ww / (uu - vv))
        z = z / 2 * p0
        If z < 0 Then z = z + 90
        If ww < 0 Then z = z + 180
        zz(j) = dms(z)                                '点位误差椭圆元素Z
    Next j
End Sub

'
'
'
Private Function 显示平差数据(ByVal fname$)
    Dim linesfromFile$, Nextline$
    
    Open fname For Input As #1
        Do Until EOF(1)
            Line Input #1, Nextline
            linesfromFile = linesfromFile + Nextline + Chr(13)
        Loop
    Close #1
    Print linesfromFile
End Function

'
'
'
Private Sub 显示平差结果()
    Print "平差结果："
    Print "---------------------------------------------------------------------------------------------"
    Print " 单位权中误差 m0=" & Format(m0, "####.##") & " 秒"
    Print "---------------------------------------------------------------------------------------------"
    Print "点名", "x(m)", "y(m)", "点位误差(cm)", "E(cm)", "F(cm)", "Z(度分秒)"
    For i = 1 To n - 1
        Print " " & Po(i).name, Format(Po(i).x, "########.####"), Format(Po(i).y, "########.####"), "  " & Format(mk(i), "#####.##"), _
                  " " & Format(e(i), "#####.##"), Format(f(i), "#####.##"), Format(zz(i), "###.#####")
    Next i
    Print "---------------------------------------------------------------------------------------------"
    '
    '
    '
    '可以将平差所得的B点坐标和B点的实际坐标输出进行比较，以验证平差结果是否正确
    'Print Format(Po(n).x, "########.######"), Format(Po(n).y, "########.######")
    'Print B_x, B_y
End Sub