后方交会实验报告
- 格式:doc
- 大小:414.50 KB
- 文档页数:17
摄影测量学实验报告实验名称:空间后方交会班级:测绘工程12-1一.实验目的掌握运用空间后方交会求解外方元素的过程。
学会运用空间后方交会的原理,根据所给控制点的地面摄影测量坐标系坐标以及相应的像平面坐标系中的坐标,利用已知公式和计算机编程语言实现空间后方交会的过程,完成外方位元素的求解。
二.实验工具个人电脑,visual basic 6.0三.实验源代码Private Sub Command1_Click()Dim c(), kx(), ky(), l(), aa(), a11(), a12(), a13(), a14(), a15(), a16(), a21(), a22(), a23(), a24(), a25()Dim a26(), v(), mm(5)nn = Val(Text2.Text)nnt = nn - 1ReDim c(nnt, 5), kx(nnt), ky(nnt), l(2 * nnt + 1, 0), aa(2 * nnt + 1, 5), a11(nnt), a12(nnt), a13(nnt), a14(nnt), a15(nnt), a16(nnt), a21(nnt), a22(nnt), a23(nnt), a24(nnt), a25(nnt)ReDim a26(nnt)f = Val(Text1.Text)fei = 0#w = 0#k = 0#If Label2.Caption = "" ThenMsgBox "尚未读取文本", , "提示信息"Elsea = Split(Label2.Caption, vbCrLf)For i = 0 To nnthh = Split(a(i), " ")For j = 0 To 5c(i, j) = Val(hh(j))Next jNext im1 = (c(0, 1) - c(1, 1)) ^ 2 + (c(0, 2) - c(1, 2)) ^ 2 m2 = (c(0, 3) - c(1, 3)) ^ 2 + (c(0, 4) - c(1, 4)) ^ 2 m = Sqr(m2 / m1)zs = m * fxs = 0ys = 0For i = 0 To nntxs = xs + c(i, 3)ys = xs + c(i, 4)Next ixs = xs / nnys = ys / nnDoa1 = Cos(fei) * Cos(k) - Sin(fei) * Sin(w) * Sin(k)a2 = -Cos(fei) * Sin(k) - Sin(fei) * Sin(w) * Cos(k)a3 = -Sin(fei) * Cos(w)b1 = Cos(w) * Sin(k)b2 = Cos(w) * Cos(k)b3 = -Sin(w)c1 = Sin(fei) * Cos(k) + Cos(fei) * Sin(w) * Sin(k)c2 = -Sin(fei) * Sin(k) + Cos(fei) * Sin(w) * Cos(k)c3 = Cos(fei) * Cos(w)For i = 0 To nntkx(i) = -f * (a1 * (c(i, 3) - xs) + b1 * (c(i, 4) - ys) + c1 * (c(i, 5) - zs)) / (a3 * (c(i, 3) - xs) + b3 * (c(i, 4) - ys) + c3 * (c(i, 5) - zs))ky(i) = -f * (a2 * (c(i, 3) - xs) + b2 * (c(i, 4) - ys) + c2 * (c(i, 5) - zs)) / (a3 * (c(i, 3) - xs) + b3 * (c(i, 4) - ys) + c3 * (c(i, 5) - zs))Next iFor i = 0 To nntl(i * 2, 0) = c(i, 1) - kx(i)l(i * 2 + 1, 0) = c(i, 2) - ky(i)Next iFor i = 0 To nntzb = (a3 * (c(i, 3) - xs) + b3 * (c(i, 4) - ys) + c3 * (c(i, 5) - zs))a11(i) = (a1 * f + a3 * c(i, 1)) / zba12(i) = (b1 * f + b3 * c(i, 1)) / zba13(i) = (c1 * f + c3 * c(i, 1)) / zba21(i) = (a2 * f + a3 * c(i, 2)) / zba22(i) = (b2 * f + b3 * c(i, 2)) / zba23(i) = (c2 * f + c3 * c(i, 2)) / zba14(i) = c(i, 2) * Sin(w) - (c(i, 1) * (c(i, 1) * Cos(k) - c(i, 2) * Sin(k)) / f + f * Cos(k)) * Cos(w)a15(i) = -f * Sin(k) - c(i, 1) / f * (c(i, 1) * Sin(k) + c(i, 2) * Cos(k))a16(i) = c(i, 2)a24(i) = -c(i, 1) * Sin(w) - (c(i, 2) * (c(i, 1) * Cos(k) - c(i, 2) * Sin(k)) / f - f * Sin(k)) * Cos(w)a25(i) = -f * Cos(k) - c(i, 2) / f * (c(i, 1) * Sin(k) + c(i, 2) * Cos(k))a26(i) = -c(i, 1)Next iFor i = 0 To nntaa(2 * i, 0) = a11(i)aa(2 * i, 1) = a12(i)aa(2 * i, 2) = a13(i)aa(2 * i, 3) = a14(i)aa(2 * i, 4) = a15(i)aa(2 * i, 5) = a16(i)aa(2 * i + 1, 0) = a21(i)aa(2 * i + 1, 1) = a22(i)aa(2 * i + 1, 2) = a23(i)aa(2 * i + 1, 3) = a24(i)aa(2 * i + 1, 4) = a25(i)aa(2 * i + 1, 5) = a26(i)Next iReDim at(5, 2 * nnt + 1), ata(5, 5), nata(5, 5), nataat(5, 2 * nnt + 1), detx(5, 0)at = zhuanzhi(aa)ata = cheng(at, aa)nata = qiuni(ata)nataat = cheng(nata, at)detx = cheng(nataat, l)xs = xs + detx(0, 0)ys = ys + detx(1, 0)zs = zs + detx(2, 0)fei = fei + detx(3, 0)w = w + detx(4, 0)k = k + detx(5, 0)Loop Until Abs(detx(0, 0)) < 0.000001 And Abs(detx(1, 0)) < 0.000001 And Abs(detx(2, 0)) < 0.000001 And Abs(detx(3, 0)) <0.000001 And Abs(detx(4, 0)) < 0.000001 And Abs(detx(5, 0)) < 0.000001ReDim v(2 * nnt + 1, 0), ax(2 * nnt + 1, 0)ax = cheng(aa, detx)v = jian(ax, l)vv = 0For i = 0 To 2 * nnt + 1vv = v(i, 0) ^ 2 + vvNext imo = Sqr(vv / (2 * nn - 6))For i = 0 To 5mm(i) = mo * Sqr(nata(i, i))Next itxtshow.Text = txtshow.Text & "外方元素为中误差为" & vbCrLf & "Xs " & xs & " " & mm(0) & vbCrLf & "Ys " & ys & " " & mm(1) & vbCrLf & "Zs " & zs & " " & mm(2) & vbCrLf & "Φ" & fei & " " & mm(3) & vbCrLf & "W " & w & " " & mm(4) & vbCrLf & "K " & k & " " & mm(5) & vbCrLf & vbCrLftxtshow.Text = txtshow.Text & " 报告时间:" & Year(Date) & Format(Month(Date), "00") & Format(Day(Date), "00") & vbCrLftxtshow.Text = txtshow.Text & " 审批人签字:"End IfEnd SubPrivate Sub Command2_Click()CommonDialog1.Filter = "文档|*.txt"CommonDialog1.ShowOpenLabel1.Caption = CommonDialog1.FileNameEnd SubPrivate Sub Command3_Click()Label2.Caption = ""Dim mylineIf Label1.Caption <> "" ThenOpen Label1.Caption For Input As #1Do While Not EOF(1)Line Input #1, mylineLabel2.Caption = Label2.Caption + myline + vbCrLfLoopClose #1End SubPrivate Sub Command4_Click()CommonDialog1.DialogTitle = "保存文件" CommonDialog1.Filter = "文档|*.doc" CommonDialog1.Action = 2If CommonDialog1.FileName <> "" ThenOpen CommonDialog1.FileName For Output As #1Print #1, txtshow.TextClose #1End IfEnd SubPrivate Sub Command5_Click()txtshow.Text = ""End SubPrivate Sub Command6_Click()If MsgBox("确认退出程序?", 68, "提示信息") = 6 Then EndEnd IfEnd SubPublic Function cheng(a(), b())ReDim s(UBound(a, 1), UBound(b, 2))If UBound(a, 2) <> UBound(b, 1) Then MsgBox "不符合乘法要求"EndElseFor i = 0 To UBound(s, 1)For j = 0 To UBound(s, 2)s(i, j) = 0For Q = 0 To UBound(a, 2)s(i, j) = Val(s(i, j)) + Val(a(i, Q)) * Val(b(Q, j)) Next QNext jNext icheng = sEnd IfEnd FunctionPublic Function zhuanzhi(a())Dim sReDim s(UBound(a, 2), UBound(a, 1))For i = 0 To UBound(a, 1)For j = 0 To UBound(a, 2)s(j, i) = a(i, j)Next jNext izhuanzhi = sEnd FunctionPublic Function qiuni(m())Dim s()If UBound(m, 1) <> UBound(m, 2) Then MsgBox ("求逆错误")EndElsen = UBound(m, 1)Dim e(), m1, m2, m3ReDim e(n, 2 * n + 1)ReDim s(n, 2 * n + 1)For i = 0 To nFor j = 0 To ne(i, j) = m(i, j)Next jNext iFor i = 0 To nFor j = n + 1 To 2 * n + 1e(i, j) = 0Next jNext iFor i = 0 To ne(i, n + i + 1) = 1 Next iFor t = 0 To nIf e(t, t) = 0 ThenFor i = t + 1 To nFor j = 0 To 2 * n + 1 s(i, j) = e(i, j)e(i, j) = e(t, j)e(t, j) = s(t, j)Next jIf e(t, t) <> 0 Then GoTo daima1End IfNext iIf e(t, t) = 0 Then MsgBox ("求逆错误") EndGoTo lastlineEnd Ifdaima1:m1 = e(t, t)For j = 0 To 2 * n + 1e(t, j) = e(t, j) / m1Next jFor i = t + 1 To nm2 = e(i, t)For j = 0 To 2 * n + 1e(i, j) = e(i, j) - m2 * e(t, j) Next jNext iNext tFor t = 0 To n - 1For i = t + 1 To nm3 = e(t, i)For j = 0 To 2 * n + 1e(t, j) = e(t, j) - m3 * e(i, j) Next jNext iNext tReDim c(n, n)For i = 0 To nFor j = 0 To nc(i, j) = Round(e(i, j + n + 1), 2)Next jNext iqiuni = cEnd Iflastline:End FunctionPublic Function jian(a(), b())Dim s, tReDim s(UBound(a, 1), UBound(a, 2))If UBound(a, 1) <> UBound(b, 1) Or UBound(a, 2) <> UBound(b, 2) ThenMsgBox "不符合减法要求"EndElseFor i = 0 To UBound(a, 1)For j = 0 To UBound(a, 2)s(i, j) = Val(a(i, j)) - Val(b(i, j))Next jNext ijian = sEnd IfEnd FunctionPrivate Sub Form_Load()Form1.Image1.Top = 0Form1.Image1.Left = 0Form1.Image1.Width = Form1.Width Form1.Image1.Height = Form1.Height End SubPrivate Sub Form_Resize()Form1.Image1.Top = 0Form1.Image1.Left = 0Form1.Image1.Width = Form1.Width Form1.Image1.Height = Form1.Height End Sub四.实验框图获取已知数据量测控制点的像点坐标确定未知数初始值五.实验数据1 -0.08615 -0.06899 36589.41 25273.32 2195.172 -0.05340 0.08221 37631.08 31324.51 728.693 -0.01478 -0.07663 39100.97 24934.98 2386.54 0.01046 0.06443 40426.54 30319.81 757.31 六.实验截图七.实验心得此次实验让我更加了解空间后方交会的计算过程,加强了对空间后方交会基本公式和误差方程式,法线方程式的记忆。