【VB.net】测量平差程序设计（条件平差和间接平差）

先说一下，如果是陶老师的学生，建议不要抄！

1.窗口

依然是见仁见智，打码的部分是我的名字学号。

2.Form代码

Imports System.IO
Public Class Form1
    '//***********************************文件夹读取**************************************************//
    Private Sub Button1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button1.Click
        OpenFileDialog1.InitialDirectory = "C:\"
        OpenFileDialog1.Filter = "文本文件(*.txt)|*.txt|所有文件(*.*)|*.*"
        OpenFileDialog1.FilterIndex = 1
        OpenFileDialog1.RestoreDirectory = True
        If OpenFileDialog1.ShowDialog = DialogResult.OK Then
            TextFileName.Text = OpenFileDialog1.FileName
        End If
        If TextFileName.Text <> "" Then
            Dim fs As FileStream
            fs = New FileStream(TextFileName.Text, FileMode.OpenOrCreate, FileAccess.Read)
            Dim sr As StreamReader
            sr = New StreamReader(fs)
            Text1.Text = sr.ReadToEnd
            sr.Close()
            fs.Close()
        End If
    End Sub
    '//************************************条件平差*******************************************//
    Private Sub Button2_Click(sender As Object, e As EventArgs) Handles Button2.Click
        str1 = Text1.Text
        M1 = Split(Trim(str1), vbCrLf)
        '/******************************************/
        M2 = M1(0).Trim().Split(New Char() {" ", ",", "，"})            '结果为（n=4,t=2,y=2）
        M3 = M2(0).Split("=")                                           '结果为（n,4）
        n = CInt(M3(1))                                 '读取观测个数n
        M3 = M2(1).Split("=")
        t = CInt(M3(1))                       '读取必要观测个数t，即待定点个数
        r = n - t                               '计算多余观测量=n-t
        M3 = M2(2).Split("=")
        y = CInt(M3(1))                         '读取已知点个数y
        Dim H1(y - 1), S(n - 1), X(t - 1), h(n - 1), L(n - 1) As Double
        Dim A(t - 1, n - 1) As Double                 'A矩阵t行n列
        Dim A_1(n - 1, t - 1) As Double              'A的逆矩阵
        Dim Q(n - 1, n - 1), P(n - 1, n - 1), VP(0, n - 1), V(n - 1, 0), V_1(n - 1, t - 1), V_Z(0, n - 1) As Double
        Dim NAA(t - 1, t - 1), NAA_1(t - 1, n - 1), NAA_n(t - 1, t - 1), K1(t - 1, 0), W(t - 1, 0) As Double
        For i = 0 To y - 1
            M2 = Split(Trim(M1(i + 2)), "=")      '取已知点数据,即H1数组中的数值
            H1(i) = CDbl(M2(1))
        Next
        '系数矩阵A的读取
        For i = 0 To t - 1
            For j = 0 To n - 1
                M4 = M1(y + 3 + i).Trim().Split(New Char() {" ", ",", "，"})
                A(i, j) = CDbl(M4(j))
            Next
        Next
        A_1 = ZZ(A, A_1)                             '转置
        For i = 0 To n - 1
            M5 = M1(y + t + 5 + i).Trim().Split(New Char() {" ", ",", "，"})
            h(i) = CDbl(M5(1))
            S(i) = CDbl(M5(2))                           'S，用于定权
        Next
        For k = 0 To n - 1
            For j = 0 To n - 1
                If k = j Then
                    Q(j, k) = S(k)
                Else Q(j, k) = 0
                End If
            Next
        Next                                            'Q矩阵
        P = QN(Q)
        'W矩阵的读取
        For i = 0 To t - 1
            M6 = M1(y + t + 3).Trim().Split(New Char() {" ", ",", "，"})
            W(i, 0) = CDbl(M6(i))
        Next
        NAA_1 = Mmul(A, Q, NAA_1)
        NAA = Mmul(NAA_1, A_1, NAA)
        NAA_n = QN(NAA)                               'NAA的逆矩阵
        For i = 0 To t - 1
            For j = 0 To t - 1
                NAA_n(i, j) = -NAA_n(i, j)              '负的NAA
            Next
        Next
        K1 = Mmul(NAA_n, W, K1)                      '计算K矩阵
        V_1 = Mmul(Q, A_1, V_1)
        V = Mmul(V_1, K1, V)                        '计算V矩阵
        For i = 0 To n - 1
            L(i) = h(i) + V(i, 0) / 1000
        Next                                       '计算观测量的平差值L
        M2 = Split(Trim(M1(2)), "=")
        For k = 0 To t - 1
            M7 = M1(y + t + 6 + n + k).Trim().Split(New Char() {" ", ",", "，"})
            X(k) = M2(1)
            For i = 0 To n - 1
                xx = M7(i + 1) * L(i)       '这一步是将与未知点计算有关的平差值带入，有关为0，无关为1
                X(k) = X(k) + xx
            Next
        Next
        '计算单位权中误差
        V_Z = ZZ(V, V_Z)
        VP = Mmul(V_Z, P, VP)
        sgm = Mmul(VP, V, sgm)
        sigm = Math.Sqrt(sgm(0, 0) / r)
        Text2.Text = "" + "未知点平差后高程值：（单位，米）" + vbCrLf
        For j = 0 To t - 1
            M7 = M1(y + t + 6 + n + j).Trim().Split(New Char() {" ", ",", "，"})
            Text2.Text = Text2.Text + M7(0) + "：" + CStr(X(j)) + vbCrLf
        Next
        Text2.Text = Text2.Text + "单位权中误差为：（单位：毫米）" + vbCrLf
        Text2.Text = Text2.Text + CStr（sigm） + vbCrLf

    End Sub
    Private Sub Button4_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button4.Click
        Dim sf As New SaveFileDialog
        sf.Filter = "文本文件(*.txt)|*.txt|所有文件(*.*)|*.*"
        sf.ShowDialog()
        My.Computer.FileSystem.WriteAllText（sf.FileName， Text2.Text， False）
        MsgBox("已保存文件到:" & sf.FileName)
    End Sub
    Private Sub Button5_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button5.Click
        Dim sf As New SaveFileDialog
        sf.Filter = "文本文件(*.txt)|*.txt|所有文件(*.*)|*.*"
        sf.ShowDialog()
        My.Computer.FileSystem.WriteAllText（sf.FileName， Text1.Text， False）
        MsgBox("已保存文件到:" & sf.FileName)
    End Sub

    Private Sub Button3_Click(sender As Object, e As EventArgs) Handles Button3.Click
        str2 = Text1.Text
        M1 = Split(Trim(str2), vbCrLf)
        M2 = M1(0).Trim().Split(New Char() {" ", ",", "，"})            '结果为（n=4,t=1,y=4）
        M3 = M2(0).Split("=")                                           '结果为（n,4）
        n = CInt(M3(1))                                 '读取观测个数n,即方程数
        M3 = M2(1).Split("=")
        t = CInt(M3(1))                         '读取必要观测个数t，即参数个数
        r = n - t                               '计算多余观测量=n-t
        Dim B(n - 1, t - 1), B_Z(t - 1, n - 1), BZP(t - 1, n - 1), V(n - 1, 0), VZ(0, n - 1), xx(t - 1, 0), BX(n - 1, 0) As Double
        Dim P(n - 1, n - 1), NBB(t - 1, t - 1), NBB_N(t - 1, t - 1), W(t - 1, 0), L(n - 1, 0), VZP(0, n - 1) As Double
        Dim h(n - 1), S(n - 1), X(t - 1)， L_1(n - 1) As Double
        '矩阵L的读取
        M4 = Split(M1(n + 3), ";")
        For j = 0 To n - 1
            L(j, 0) = CDbl(M4(j))
        Next
        '系数矩阵B的读取
        For i = 0 To n - 1
            For j = 0 To t - 1
                M5 = M1(n + 5 + i).Trim().Split(New Char() {" ", ",", "，"})
                B(i, j) = CDbl(M5(j))
            Next
        Next
        B_Z = ZZ(B, B_Z)                             'B的转置
        For i = 0 To n - 1
            M6 = M1(2 + i).Trim().Split(New Char() {" ", ",", "，"})
            h(i) = CDbl(M6(1))
            S(i) = CDbl(M6(2))                           'S，用于定权
        Next
        'P矩阵的读取
        For k = 0 To n - 1
            For j = 0 To n - 1
                If k = j Then
                    P(j, k) = 10 / S(k)
                Else P(j, k) = 0
                End If
            Next
        Next
        BZP = Mmul(B_Z, P, BZP)
        NBB = Mmul(BZP, B, NBB)
        W = Mmul(BZP, L, W)
        NBB_N = QN(NBB)
        xx = Mmul(NBB_N, W, xx)
        BX = Mmul(B, xx, BX)
        For k = 0 To n - 1
            V(k, 0) = BX(k, 0) - L(k, 0)
        Next
        For i = 0 To t - 1
            M7 = M1(2 * n + i + 6).Trim().Split("=")
            X(i) = CDbl(M7(1)) + xx(i, 0) / 1000
        Next
        For j = 0 To n - 1
            L_1(j) = h(j) + V(j, 0) / 1000
        Next
        VZ = ZZ(V, VZ)
        VZP = Mmul(VZ, P, VZP)
        sgm = Mmul(VZP, V, sgm)
        sigm = Math.Sqrt(sgm(0, 0) / r)
        Text2.Text = "" + "未知点平差后高程值：（单位，米）" + vbCrLf
        For j = 0 To t - 1
            M8 = M1(2 * n + 6 + j).Trim().Split("=")
            Text2.Text = Text2.Text + M8(0) + "：" + CStr(X(j)) + vbCrLf
        Next
        Text2.Text = Text2.Text + "观测值的平差值：（单位，米）" + vbCrLf
        For k = 0 To n - 1
            M9 = M1(2 + k).Trim().Split(New Char() {" ", ",", "，"})
            Text2.Text = Text2.Text + M9(0) + "：" + CStr(L_1(k)) + vbCrLf
        Next
        Text2.Text = Text2.Text + "单位权中误差为：（单位：毫米）" + vbCrLf
        Text2.Text = Text2.Text + CStr（sigm） + vbCrLf
    End Sub
End Class

3.Module代码

Module Module1
    Public str1, str2 As String
    Public y, n, t, i, j, k, r As Integer
    Public xx, sgm(0, 0), sigm As Double
    Public M1(), M2(), M3(), M4(), M5(), M6(), M7(), M8(), M9() As String
    '//********************************下面是函数部分*********************************//
    '矩阵相乘函数
    Function Mmul(mtxA(,) As Double, mtxB(,) As Double, ByRef mtxC(,) As Double)
        Dim m As Integer
        Dim n As Integer
        Dim l As Integer
        Dim i As Integer, j As Integer, K As Integer
        m = UBound(mtxA, 1) - LBound(mtxA, 1)
        n = UBound(mtxA, 2) - LBound(mtxA, 2)
        l = UBound(mtxB, 2) - LBound(mtxB, 2)
        For i = 0 To m
            For j = 0 To l
                mtxC(i, j) = 0#
                For K = 0 To n
                    mtxC(i, j) = mtxC(i, j) + mtxA(i, K) * mtxB(K, j)
                Next K
            Next j
        Next i
        Return mtxC
    End Function
    '矩阵转置函数
    Function ZZ(mtxA(,) As Double, mtxB(,) As Double)
        Dim m As Integer
        Dim n As Integer
        Dim i As Integer, j As Integer
        m = UBound(mtxA, 1) - LBound(mtxA, 1)       '原来的行数
        n = UBound(mtxA, 2) - LBound(mtxA, 2)      '原来的列数
        For i = 0 To m
            For j = 0 To n
                mtxB(j, i) = mtxA(i, j)
            Next j
        Next i
        Return mtxB
    End Function
    '矩阵求逆函数
    Function QN(Mtx(,) As Double)
        Dim N As Double
        N = UBound(Mtx, 1) - LBound(Mtx, 1)               '矩阵的行数/列数
        Dim row(0 To N) As Integer, col(0 To N) As Integer
        Dim d As Double, p As Double
        For k = 0 To N
            d = 0.0#
            For i = k To N
                For j = k To N
                    p = Math.Abs(Mtx(i, j))
                    If (p > d) Then
                        d = p
                        row(k) = i
                        col(k) = j    '记录变换的位置
                    End If
                Next j
            Next i
            If (d + 1.0# = 1.0#) Then '找到主元，并存储到行列号
            End If
            If (row(k) <> k) Then
                For j = 0 To N
                    p = Mtx(k, j)
                    Mtx(k, j) = Mtx(row(k), j)
                    Mtx(row(k), j) = p '行变换
                Next j
            End If
            If (col(k) <> k) Then
                For i = 0 To N
                    p = Mtx(i, k)
                    Mtx(i, k) = Mtx(i, col(k))
                    Mtx(i, col(k)) = p '列变换
                Next i
            End If
            Mtx(k, k) = 1.0# / Mtx(k, k)
            For j = 0 To N
                If (j <> k) Then
                    Mtx(k, j) = Mtx(k, j) * Mtx(k, k)
                End If
            Next j
            For i = 0 To N
                If (i <> k) Then
                    For j = 0 To N
                        If (j <> k) Then
                            Mtx(i, j) = Mtx(i, j) - Mtx(i, k) * Mtx(k, j)
                        End If
                    Next j
                End If
            Next i
            For i = 0 To N
                If (i <> k) Then
                    Mtx(i, k) = -Mtx(i, k) * Mtx(k, k)
                End If
            Next i
        Next k
        For k = N To 0 Step -1
            If (col(k) <> k) Then
                For j = 0 To N
                    p = Mtx(k, j)
                    Mtx(k, j) = Mtx(col(k), j)
                    Mtx(col(k), j) = p
                Next j
            End If
            If (row(k) <> k) Then
                For i = 0 To N
                    p = Mtx(i, k)
                    Mtx(i, k) = Mtx(i, row(k))
                    Mtx(i, row(k)) = p
                Next i
            End If
        Next k
        Return Mtx
    End Function
End Module

当时费劲切片分数据的做法真的笨的可以ε＝ε＝ε＝(#>д<)ﾉ