误差理论与测量平差课程设计报告

1、一 设计原始资料水准网严密平差及精度评定例如。如下图水准网，有2个点，3个未知点，7个测段。各数据及观测值见下表（1）（2） 高差观测值(m) 高差观测值(m) 端点号高差观测值测段距离序号1-311-422-332-443-453-565-2-0.59573求各待定点的高程；3-4点的高差中误差；3号点、4号点的高程中误差。提示，本网可采用以测段的高差为平差元素，采用间接平差法编写程序计算。二设计内容及要求正确应用平差模型列出观测值条件方程、误差方程、法方程和解算法方程，得出平差后的平差值及各待定点的高程平差值，评定各平差值的精度和各高程平差值的精度。三水准网间接平差思路根据网型确定水准点数

2、nn，未知水准点数un，总点数tn,总的观测高差段数hn，必要观测数，多余观测数。.确定参数。为平差后能直接求得待定点高程平差值，以3个待定点高程平差值为参数。设3,4,5点的高程平差值分别为， 。列立条件方程.左侧为观测值系数为1，右侧为参数和常数项，并进一步改化成误差方程，最终写成矩阵形式。得到系数矩阵A和常数项L.列立法方程，并解求法方程。由于该水准网间接平差误差方程个数为7个而未知数个数为10个，所列的误差方程是一组相容方程，有无数组解，所以必须在最小二乘原那么VTPV=min的根底上利用拉格朗日乘数法求解.令F= VTPV-2KT(V-A+L),分别对V和求导，并令其导数为零，得到2

3、VTP-2KT=0,ATK=0，将二式合并即得法方程：ATPV=ATPA-ATPL=0。求出Naa= ATPA,W= ATPL,即得到相应的法方程。求解法方程，得到=N-1aaW加上Xi0即可得到待定点的高程平差值，将代入误差方程得到相应的V值，hi+Vi得到各段高差的平差值。精度评定。计算单位权中误差的估值：评定各待定点的高程中误差： 各待定点的精度即为：评定高程平差值的精度：步骤1误差方程的列立。在该水准网中，总观测数n=7，必要观测数t=3,多余观测数r=n-t=7-3=43、4和5点高程平差值为参数分别为：X1、X2、 X3，相应的近似值取为：按图示列出观测方程,化为误差方程如下： 代

4、入写成矩阵形式为：2组成法方程并结算法方程。以1km水准测量的观测高差为单位权观测值，各观测值相互独立，权阵为一主对角阵，由 定权 得权阵为：那么法方程为： 其中 ， 由B，P，,L求出：解得：3)计算V和 1将 代入得：2由得： 3由 得: m4)精度评定。 1计算单位权中误差的估值。2计算3号点、4号点的高程中误差。得 2计算3号点4号点的高差中误差。3至4间高差平差值 的权函数式为 五平差程序设计思路1、数据文件的编辑。 在文本文档输入数据，包括水准网中点数、未知点数以及这些点的点号、高程和高差观测值、距离观测值等。约定文件输入的格式如下：第一行：点个数、未知点个数、观测值个数第二行：点

5、号点在前，未知点在后第三行：高程顺序与上一行的点号对应第四行：高差观测值，按“起点点号、终点点号、高差观测值、距离观测值的顺序输入。本算例数据文件如下：2，3，71，2，3，4，502、界面设计用菜单组织程序，用文本框显示数据的输入、计算和输出情况。菜单设计。本程序的菜单结构如下列图所示：标题名称快捷键文件&FilemnuFile翻开数据mnuOpenCtrl+O保存结果mnuSaveCtrl+S-Aa退出mnuQuitCtrl+Q计算&CalcmnuCalcCtrl+C近似高程mnuHightCtrl+H误差方程mnuEquCtrl+E平差计算mnuAdjCtrl+A精度评定

6、mnuEva窗体、文本框和通用对话框。在主窗体上绘制1个文本框控件和一个通用对话框控件，工具箱中假设无通用对话框那么需在工程部件中找到Microsoft Common Dialog Control 6.0,选中确定，调出通用对话框控件。如下列图。按下表设置各控件的相关属性.对象属性值Text1NametxtShowText1TextText1MultiLineTrueText1ScrollBars3-BothForm1Caption水准网间接平差CommonDialog！NameCDg1界面设计如下列图：3、代码设计1。公共变量声明。 Dim strFileName As StringDim

7、nn%, un%, tn%, hn% '点个数，未知点个数，总点数，观测值个数的定义Dim Pname() As String '点名数组Dim Hknown() As Double '高程数组，存放点高程和高程近似值Dim be%(), en%() '观测值的起点和终点编号数组，储存的是点的序号Dim h#(), s#() '高差观测值数组和距离观测值数组Dim A#(), X#(), P#(), L#() '间接平差的系数阵，解向量，权阵，常数向量2。间接平差通用过程代码设计。 间接平差的通用过程代码在标准模块2中设计，将标准模块加载到过程中

8、，即可在平差计算时直接调用该过程。间接平差通用过程代码设计主要包括以下几步：，定义矩阵转置子过程MatrixTrans(A#(),AT#()。其中A#()中存放原始矩阵，AT#中存放其转置矩阵。，定义矩阵相乘子过程MatrixMulty(a#(),b#(),c#()。其中a#(),b#()为原始矩阵，c#()存放二矩阵相乘所得到的结果。，定义矩阵求逆子过程MatInv(a#().用Gauss-Jordan法求矩阵a#()的逆矩阵，假设a#()是奇异阵，那么提示相关信息，并终止程序运行，假设a#()不是奇异阵，那么a#()中的数据被改写，实矩阵a#()中存放其逆矩阵的计算结果。,调用以上子过程实

9、现间接平差通用过程InAdjust(A#(),P#(),L#(),X#()的编辑。 Sub InAdjust(A#(), P#(), L#(), X#() Dim a1%, a2%, p1%, p2%, l1%, x1% Dim At#(), Atp#(), Naa#(), Na#(), W#() '令Na为Naa的逆矩阵 a1 = UBound(A, 1) - LBound(A, 1) + 1: a2 = UBound(A, 2) - LBound(A, 2) + 1 l1 = UBound(L) - LBound(L) + 1: x1 = UBound(X) - LBound(X)

10、 + 1 p1 = UBound(P, 1) - LBound(P, 1) + 1: p2 = UBound(P, 2) - LBound(P, 2) + 1 ReDim At(1 To a2, 1 To a1), Atp(1 To a2, 1 To a1), Naa(1 To a2, 1 To a2), W(1 To a2) MatrixTrans A, At '调用子过程求出系数矩阵的转置矩阵 MatrixMulty At, P, Atp '调用子过程求出Atp MatrixMulty Atp, A, Naa '调用子过程求出Naa MatrixMulty Atp,

11、 L, W '调用子过程求出W matInv Naa 'Naa假设可逆Naa中的数据被改写，实矩阵Naa中存放Naa逆的结果 Na = Naa MatrixMulty Na, W, XEnd Sub3、数据的输入。 单击“文件翻开文件命令，在mnuOpen_Click()事件中编写代码，并通过给通用对话框设置相应的代码实现翻开过程，弹出翻开对话框。其中CDg1.Filter = "文本文件(*.txt)|*.txt|所有文件(*.*)|*.*"用来指定在通用对话框中列出的文件类型。，待用户选取了文件以后，程序开始读取数据。再对文本框设置相关代码，通过循环过程

12、实现数据的显示，翻开过程和显示结果如下列图：、近似高程的计算。输入数据后，点击“计算近似高程在 mnuHeight_Click()事件中编写代码。用一个数组来存储高程近似值，点的高程放在这个数组的开头，然后按照点号输入顺序依次搜索涉及该点的高差观测值，看该高差涉及的另一点是否，假设未知，那么将检查下一个高差观测值。假设，那么可以计算出当前未知点的高差近似值，并放入高程近似值数组，以此类推，直到所有未知点的高程近似值都被求出为止。最后通过文本框显示出来。显示结果如下列图所示：、列立观测值的误差方程。点击“计算误差方程在mnuEqu_click()事件中编写代码，程序根据输入的各观测值的起止点信息

13、及高差、距离值，根据高差方向，由于每个高差由相应参数表示时，参数前的系数不是1就是-1，所以用循环内嵌If语句可以得出误差方程的系数矩阵A()如下：For i = 1 To hnIf en(i) > nn Then A(i, en(i) - nn) = 1 '假设终点未知，那么给终点对应的系数赋值If be(i) > nn Then A(i, be(i) - nn) = -1 '假设起点未知，那么给起点对应的系数赋值Next i并根据Pi=C/Si定权，令C=1km,得权矩阵P()，并根据从JK段，L=H(K)-H(J)-h(i)得到常数项L()，并将结果显示在文本

14、框中。显示结果如下列图：、平差解算。点击“计算平差计算命令，调用间接平差通用过程实现间接平差计算。并求出高程平差值，并显示在文本框中。、精度评定。点击“计算精度评定命令，在mnuEva_Click()中编写代码，求出高程中误差，并显示在文本框中。、保存和退出。、文件的保存点击“文件保存结果命令，在mnuSave_Click()事件中编码，将文本框中结果保存在指定的文件中。并设置CDg1.Filter = "文本文件(*.txt)|*.txt|所有文件(*.*)|*.*"使其以文本文件的形式保存。、文件的退出点击“文件退出命令，在mnuQuit_Click()事件中编码，退出

15、整个程序。五、平差程序流程图数据文件编辑界面设计代码设计数据输入近似高程计算误差方程列立公共变量声明高程平差值解算精度评定保存,退出菜单文本框通用对话框六、程序源代码及其说明'公共变量的声明Dim strFileName As StringDim nn%, un%, tn%, hn% '点个数，未知点个数，总点数，观测值个数的定义Dim Pname() As String '点名数组Dim Hknown() As Double '高程数组，存放点高程和高程近似值Dim be%(), en%() '观测值的起点和终点编号数组，储存的是点的序号Dim h#(

16、), s#() '高差观测值数组和距离观测值数组Dim A#(), X#(), P#(), L#() '间接平差的系数阵，解向量，权阵，常数向量标准模块1 '点名序号转换函数 Public Function Order(str As String) Dim i% For i = 1 To tn If str = Pname(i) Then Order = i Exit For End If Next iEnd Function标准模块2'实现矩阵的转置Sub MatrixTrans(A() As Double, At() As Double)Dim m%, n%

17、, i%, j%m = UBound(A, 1) '提取A矩阵的行上标n = UBound(A, 2) '提取A矩阵的列上标For i = 1 To m For j = 1 To n At(j, i) = A(i, j) '求转置矩阵 Next jNext iEnd Sub '实现矩阵相乘Sub MatrixMulty(A#(), b#(), c#()Dim m%, y%, n%, g%, i%, j%, k%m = UBound(A, 1) '提取A矩阵的行上标y = UBound(A, 2) '提取A矩阵的列上标n = UBound(b, 1

18、) '提取B矩阵的行上标g = UBound(b, 2) '提取B矩阵的列上标ReDim A(i, k)ReDim b(k, j)If y <> n Then MsgBox "这两个矩阵不能相乘"End IfIf y = n ThenFor i = 1 To m For j = 1 To g For k = 1 To y c(i, j) = A(i, k) * b(k, j) '用循环实现矩阵相乘 Next k Next j Next iEnd IfEnd Sub'实现矩阵的求逆Sub matInv(A() As Double)D

19、im n%, k%, i%, j%, m#, d#Dim ii%(), jj%()n = UBound(A, 1)ReDim ii(n), jj(n)For k = 1 To n d = 0! For i = k To n For j = k To n m = Abs(A(i, j) If m > d Then d = m: ii(k) = i: jj(k) = k '记下主对角线右下角子阵绝对值最大元素的行号i列号j End If Next j Next i If (d + 1!) = 1! Then MsgBox "矩阵不能求逆！", vbCritical

20、, "出错" '奇异阵不可逆 Exit Sub End If If ii(k) <> k Then '将最大元素通过行交换，列交换交换到主元素的位置即主对角线上 For j = 1 To n m = A(k, j): A(k, j) = A(ii(k), j): A(ii(k), j) = m '行交换 Next j End If If jj(k) <> k Then For i = 1 To n m = A(i, k): A(i, k) = A(i, jj(k): A(i, jj(k) = m '列交换 Next i

21、 End If A(k, k) = 1! / A(k, k) '将1/a(k,k)赋值给a(k,k) For j = 1 To n If j <> k Then A(k, j) = A(k, k) * A(k, j) '当j<>k时将a(k,k)*a(k,j)赋值给a(k,j) Next j For i = 1 To n If i <> k Then For j = 1 To n If j <> k Then A(i, j) = A(i, j) - A(i, k) * A(k, j) '当j<>k时将(a(i,

22、j)-a(i,k)*a(k,j)赋值给a(i,j) End If Next j End If Next i For i = 1 To n If i <> k Then A(i, k) = A(i, k) * A(k, k) '当i<>k时将(a(i, k) * a(k, k)赋值给a(i,k) Next iNext kFor k = n To 1 Step -1 If jj(k) <> k Then For j = 1 To n m = A(k, j): A(k, j) = A(jj(k), j): A(jj(k), j) = m Next j En

23、d If If ii(k) <> k Then For i = 1 To n m = A(i, k): A(i, k) = A(i, ii(k): A(ii(k), j) = m Next i End If If ii(k) <> k Then For i = 1 To n m = A(i, k): A(i, k) = A(i, ii(k): A(i, ii(k) = m Next i End If Next kEnd Sub'实现间接平差的通用过程 Sub InAdjust(A#(), P#(), L#(), X#() Dim a1%, a2%, p1%, p

24、2%, l1%, x1% Dim At#(), Atp#(), Naa#(), Na#(), W#() '令Na为Naa的逆矩阵 a1 = UBound(A, 1) - LBound(A, 1) + 1: a2 = UBound(A, 2) - LBound(A, 2) + 1 l1 = UBound(L) - LBound(L) + 1: x1 = UBound(X) - LBound(X) + 1 p1 = UBound(P, 1) - LBound(P, 1) + 1: p2 = UBound(P, 2) - LBound(P, 2) + 1 ReDim At(1 To a2,

25、1 To a1), Atp(1 To a2, 1 To a1), Naa(1 To a2, 1 To a2), W(1 To a2) MatrixTrans A, At '调用子过程求出系数矩阵的转置矩阵 MatrixMulty At, P, Atp '调用子过程求出Atp MatrixMulty Atp, A, Naa '调用子过程求出Naa MatrixMulty Atp, L, W '调用子过程求出W matInv Naa 'Naa假设可逆Naa中的数据被改写，实矩阵Naa中存放Naa逆的结果 Na = Naa MatrixMulty Na, W,

26、 XEnd Sub代码窗口'数据的输入Private Sub mnuOpen_Click()Dim i%, strT1 As String, strT2 As String '定义循环变量和临时变量CDg1.Filter = "文本文件(*.txt)|*.txt|所有文件(*.*)|*.*"CDg1.ShowOpen '翻开对话框strFileName = CDg1.FileName '获得选中的文件名和路径Open strFileName For Input As #1 '翻开文件 Input #1, nn, un, hn '

27、;读入点的个数nn，未知点的个数un，观测值的个数hn tn = nn + unReDim Pname(tn), Hknown(tn), h(hn), s(hn), be(hn), en(hn)For i = 1 To tn '读入点名 Input #1, Pname(i)Next iFor i = 1 To nn '读入高程 Input #1, Hknown(i)Next iFor i = 1 To hn '读入各观测值 Input #1, strT1, strT2, h(i), s(i) be(i) = Order(strT1): en(i) = Order(str

28、T2) '给起终点数组排序Next i'显示读入的数据txtShow.Text = txtShow.Text & "读入水准网数据:" & vbCrLftxtShow = txtShow & "点" & nn & "个，未知点" & un & "个，观测值" & hn & "个。"txtShow.Text = txtShow.Text & vbCrLf & "网中涉及的点名有：&qu

29、ot; For i = 1 To tn txtShow.Text = txtShow.Text & Pname(i) & "，" Next i txtShow.Text = txtShow.Text & vbCrLf txtShow.Text = txtShow.Text & "点高程为：" & vbCrLf For i = 1 To nn txtShow.Text = txtShow.Text & Pname(i) & "的高程为：" & Hknown(i) &

30、 vbCrLf Next i txtShow.Text = txtShow.Text & "各观测值分别为： " & vbCrLf txtShow = txtShow & "起点 终点 高差观测值 距离观测值" & vbCrLf For i = 1 To hn txtShow = txtShow & Pname(be(i) & " " & Pname(en(i) & " " txtShow = txtShow & Format(h(i), &q

31、uot;0.000") & " " & Format(s(i), "0.000") & vbCrLf Next i Close #1 '关闭文件 End Sub'近似高程的计算Private Sub mnuHeight_Click() Dim i%, j% For i = 1 To un For j = 1 To hn If be(j) = nn + i And en(j) < nn + i Then '找到一个起点相同且终点的近似值 Hknown(nn + i) = Hknown(en(j

32、) - h(j) Exit For End If If en(j) = nn + i And be(j) < nn + i Then '找到一个起点相同且终点的近似值 Hknown(nn + i) = Hknown(be(j) + h(j) Exit For End If Next j Next i '显示近似高程计算结果 txtShow.Text = txtShow.Text & "近似高程计算结果：" & vbCrLf For i = 1 To un txtShow = txtShow & Pname(i + nn) &am

33、p; ":" & Format(Hknown(i + nn), "0.000") & vbCrLf Next i End Sub '列误差方程Private Sub mnuEqu_click() Dim i%, j% ReDim A(1 To hn, 1 To un), L(1 To hn), P(1 To hn, 1 To hn) For i = 1 To hn If en(i) > nn Then A(i, en(i) - nn) = 1 '假设终点未知，那么给终点对应的系数赋值 If be(i) > nn

34、 Then A(i, be(i) - nn) = -1 '假设起点未知，那么给起点对应的系数赋值 L(i) = -(Hknown(en(i) - Hknown(be(i) - h(i) '根据起终点计算常数项 P(i, i) = 1 / s(i) '以距离的导数定义权 Next i '显示误差方程 txtShow.Text = txtShow.Text & "列立的误差方程：" & vbCrLf For i = 1 To hn For j = 1 To un txtShow.Text = txtShow.Text &

35、A(i, j) & " " Next j txtShow.Text = txtShow.Text & " " & Format(L(i), "0.0000") & vbCrLf Next i txtShow.Text = txtShow.Text & "权矩阵：" & vbCrLf For i = 1 To hn For j = 1 To hn txtShow.Text = txtShow.Text & P(i, j) & " "

36、Next j txtShow.Text = txtShow.Text & vbCrLf Next iEnd Sub'平差计算并显示结果Private Sub mnuAdj_click() Dim i%, j% ReDim X(1 To un) InAdjust A, P, L, X '调用间接平差间接过程求解 txtShow.Text = txtShow.Text & "平差计算结果：" & vbCrLf '计算并显示高程平差结果 txtShow = txtShow & "点号 初始高程(m) 高程改正值(m

37、) 平差后高程(m)" & vbCrLf For i = 1 To un txtShow.Text = txtShow.Text & Pname(nn + i) & " " & Format(Hknown(nn + i), "0.0000") Hknown(nn + i) = Hknown(nn + i) + X(i) txtShow = txtShow & "" & Format(X(i), "0.000") & "" &

38、; Format(Hknown(nn + i), "0.000") & vbCrLf Next iEnd Sub '精度评定 Private Sub mnuEva_Click()3,4点高程中误差Dim m#,r%,i%Dim V#(),Vt#()MatrixMulty A,X,AXV = AXLr = hn unMatrixTrans V,VtMatrixMulty Vt,p,VtPMatrixMulty VtP,V,VtPVm = sqr ( VtPV / r )Qxx = NaFor i= 1 to r-1 求3,4点高程中误差m(xi) = m * sqr(Na(i)Next iEnd Sub'文件的保存Private Sub m