摄影测量实验报告.doc

《摄影测量实验报告.doc》由会员分享，可在线阅读，更多相关《摄影测量实验报告.doc（21页珍藏版）》请在三一办公上搜索。

1、精选优质文档-倾情为你奉上 目录 实验一 数字内定向编程实验一、目的要求： 掌握数字内定向的原理，利用计算机编程语言实现数字内定向的过程，完成一立体像对两幅像片影像的数字内定向。 二、仪器用具： 计算机，编程软件（VB）三、实验数据： 立体像对两像片（1504、1505）上四个框标点的理论位置以及四个框标点的像素坐标，具体如表3-1。 15041505ijijxyxy1417.229190.964412.029198.964-110.000109.999-110.000109.99928277.836335.6148271.936367.414109.999110.006109.999110.

2、00638130.9078195.1648101.3078226.264109.987-109.997109.987-109.9974269.8648050.136241.0648058.036-109.999-109.999-109.999-109.999表3-1 四、实验内容： 4.1 内定向参数求取 方法：利用航摄像片上的四个框标点的理论位置以及四个框标点的像素坐标为依据。 根据条件：框标点的像素坐标通过公式转换后应等于其理论坐标。 利用式（3-1）列出方程，通过最小二乘法计算内定向参数: h0、 h1、 h2、 k0、k1、 k2。 计算公式： x=h0 + h1*i + h2*j y

3、=k0 + k1*i + k2*j （3-1） 求解过程： 列误差方程 法化 求参数 4.2 将像点扫描坐标转化为框标坐标 通过公式： x=h0 + h1*i + h2*j y=k0 + k1*i + k2*j 将测量的同名点的扫描坐标转换位像框标坐标系中坐标。 五、实验程序：Private Sub read_Click() 读入数据 Dim i%, j%, k% CommonDialog1.ShowOpen Open CommonDialog1.FileName For Input As #1 For i = 1 To 4 Line Input #1, s s = Split(s, ) Fo

4、r j = 1 To 2 A(i, j + 1) = s(j - 1) Next j Line Input #1, s s = Split(s, ) Lx(i, 1) = s(0) Ly(i, 1) = s(1) Next i For i = 1 To 4 A(i, 1) = 1 Next iEnd Sub内定向求解Private Sub jisuan_Click() Dim i%, j%, k% For i = 1 To 3 For j = 1 To 4 at(i, j) = A(j, i) Next j Next i矩阵求解 JZC at(), A(), 3, 4, 3, B() Rect

5、 B(), 3, BdN() JZC BdN(), at(), 3, 3, 4, bb() JZC bb(), Lx(), 3, 4, 1, HH() JZC bb(), Ly(), 3, 4, 1, KK() For i = 1 To 3 Text1 = Text1 & Format(HH(i, 1), 0.) & vbCrLf Next i Text1 = Text1 & vbCrLf For i = 1 To 3 Text1 = Text1 & Format(KK(i, 1), 0.) & vbCrLf Next i Close #1End Sub坐标转换Private Sub Comm

6、and1_Click()xx = Val(Text2)yy = Val(Text3)Text4 = HH(1, 1) + xx * HH(2, 1) + yy * HH(3, 1)Text5 = KK(1, 1) + xx * KK(2, 1) + yy * KK(3, 1)Text4.Text = Format(Text4, 0.000)Text5.Text = Format(Text5, 0.000)End Sub六、实验结果：1504定向参数：x方向：h0=-121.,h1=0.,h2=0.,精度=0.y方向：k0=115.,k1=0.,k2=-0.,精度=0.1505定向参数：x方向：

7、h0=-121.,h1=0.,h2=0.,精度=0.y方向：k0=115.,k1=0.,k2=-0.,精度=0.实验二 空间后方交会编程实验一、目的要求： 掌握空间后方交会的原理，根据所给控制点的地面物方坐标以及相应的像点在像平面坐标系中的坐标，利用计算机编程语言实现空间后方交会的过程，完成所给立体像对（1504、1505）中两张像片各自的6个外方位元素的解求。 二、仪器用具： 计算机，编程软件（VB）三、实验数据 实验数据为一立体像对（1504、1505）的像片控制点数据，包括像点坐标及相应地面摄测坐标。具体如表4-1所示： 左片（1504）像平面坐标地面摄测坐标xyXtpYtpZtp-2.

8、81673.965.070.38014.250-6.459-87.783.140.3305.580-76.415-46.343.380.9805.430-30.041-40.404.290.8008.810-65.89075.337.750.2305.760右片（1505）像平面坐标地面摄测坐标xyXtpYtpZtp72.50774.412.070.38014.25073.384-85.905.140.3305.5803.767-47.117.380.9805.43049.064-39.709.290.8008.81010.24473.307.750.2305.760表4-1两张像片的内方位元

9、素为：x0=y0=0mm,f=210.681mm。 四、实验内容： 利用已知地面控制点数据以及相应像点坐标，编写空间后方交会程序，求解所给像对（1504、1505）中两张像片各自的6个外方位元素。 4.1 数学模型：共线条件方程式 4.2 求解过程： 1）获取已知数据： a. 像片比例尺1/m： 由像片上两点间距离与相应地面点间距之比求得。 b. 平均航高H： H=m*f+Zi/n 其中，n为已知控制点数。 c. 内方位元素x0,y0,f： d. 地面摄影测量坐标Xtp,Ytp,Ztp e. 控制点的像点坐标（表4-1） 2） 确定未知数初始值： 在竖直摄影情况下，角元素的初始值为0，即=0；

10、线元素中，Zs0=H,Xs0,Ys0的取值可用五个已知控制点坐标的平均值，即：Xs0=Xtpi/5，Ys0=Ytpi/5。 3） 计算旋转矩阵R：利用角元素值的近似值按下式计算方向余弦值，组成R阵： a1 = Cos() * Cos() - Sin() * Sin() * Sin() a2 = -Cos() * Sin() - Sin() * Sin() * Cos() a3 = -Sin() * Cos() b1 = Cos() * Sin() b2 = Cos() * Cos() b3 = -Sin() c1 = Sin() * Cos() + Cos() * Sin() * Sin()

11、c2 = -Sin() * Sin() + Cos() * Sin() * Cos() c3 = Cos() * Cos() 4） 逐点计算像点坐标的近似值： 利用未知数的近似值按共线方程式计算控制点像点坐标的近似值（x），（y） 5） 组成误差方程式： 按下式逐点计算误差方程式的系数和常数项： 系数： zba = (a3 * (Xa - Xs) + b3 * (Ya - Ys) + c3 * (Za - zs) a11 = (a1 * f + a3 * x) / zba a12 = (b1 * f + b3 * x) / zba a13 = (c1 * f + c3 * x) / zba a

12、21 = (a2 * f + a3 * y) / zba a22 = (b2 * f + b3 * y) / zba a23 = (c2 * f + c3 * y) / zba a14= y * Sin() - (x * (x * Cos() - y * Sin() / f + f * Cos() * Cos() a15= -f * Sin() - x * (x * Sin() + y * Cos() / f a16= y a24 = -x * Sin() - (y * (x * Cos() - y * Sin() / f - f * Sin() * Cos() a25 = -f * Cos(

13、) - y * (x * Sin() + y * Cos() / f a26 = -x 常数项： lx=x-（x） ly=y-（y） 用矩阵形式表示为： V=AX-l 式中： V=Vx Vy A= a11 a12 a13 a14 a15 a16 a21 a22 a23 a24 a25 a26 X=dXs dYs dZs d d d L=lx ly 所有控制点，列出误差方程式，构成总误差方程式： V=AX-l 式中： V=V1 V2 V3 V4 V5 A=A1 A2 A3 A4 A5 L=l1 l2 l3 l4 l5 6） 组成法方程式： 计算法方程的系数矩阵Na=AA与常数项U=Al。 7）

14、解求外方位元素： 根据法方程，按下式解求外方位元素改正数： X=Nani*U 并与相应的近似值求和，得到外方位元素新的近似值。 8） 检查计算是否收敛： 将求得的外方位元素的改正数与规定的限差比较，小于限差则计算终止；否则用新的近似值重复4-8的计算，直到满足要求为止。 迭代限差为： dXs1m,dYs1m,dZs1m, d0.00001,d0.00001,d0.00001 最后得出的六个外方位元素的解为： Xs=Xs0+dXs1+dXs2+ Ys=Ys0+dYs1+dYs2+ Zs=Zs0+dZs1+dZs2+ =0+d1+d2+ =0+d1+d2+ =0+d1+d2+ 五、实验程序：读取数

15、据Private Sub read_Click()Dim i%CommonDialog1.ShowOpenOpen CommonDialog1.FileName For Input As #1For i = 1 To 6 Line Input #1, s If i 1 Then s = Split(s, ) M(i - 1) = s(0) N(i - 1) = s(1) Xtp(i - 1) = s(2) Ytp(i - 1) = s(3) Ztp(i - 1) = s(4) End IfNext iEnd Sub计算Private Sub jisuan_Click()Dim i%, j%Di

16、m s1 As Single, s2 As Single, w As Single, EZ As Single, H As Single, Xs0 As Single, Ys0 As Single, Zs0 As Single, fai As Single, omg As Single, kab As SingleDim a1 As Single, a2 As Single, a3 As Single, b1 As Single, b2 As Single, b3 As Single, c1 As Single, c2 As Single, c3 As SingleDim A(10, 6) A

17、s Single, L(10, 1) As Single, P(5, 2, 6) As SingleDim at(6, 10) As Single, BB(6, 6) As Single, CC(6, 6) As Single, DD(6, 10) As Single, EE(6, 1) As SingleDim x0(5) As Single, y0(5) As Single, lx(5) As Single, ly(5) As Single, zba(5) As Singlef = 0.像片比例尺1/ms1 = Sqr(M(1) - M(2) 2 + (N(1) - N(2) 2)s2 =

18、 Sqr(Xtp(1) - Xtp(2) 2 + (Ytp(1) - Ytp(2) 2)w = s2 / s1平均航高HEZ = (Ztp(1) + Ztp(2) + Ztp(3) + Ztp(4) + Ztp(5) / 5H = EZ控制点坐标平均值Xs0 = (Xtp(1) + Xtp(2) + Xtp(3) + Xtp(4) + Xtp(5) / 5Ys0 = (Ytp(1) + Ytp(2) + Ytp(3) + Ytp(4) + Ytp(5) / 5Zs0 = EZ + w * ffai = 0omg = 0kab = 0lj = 1Do While lj 6 R矩阵系数 a1 =

19、Cos(fai) * Cos(kab) - Sin(fai) * Sin(omg) * Sin(kab) a2 = -Cos(fai) * Sin(kab) - Sin(fai) * Sin(omg) * Cos(kab) a3 = -Sin(fai) * Cos(omg) b1 = Cos(omg) * Sin(kab) b2 = Cos(omg) * Cos(kab) b3 = -Sin(omg) c1 = Sin(fai) * Cos(kab) + Cos(fai) * Sin(omg) * Sin(kab) c2 = -Sin(fai) * Sin(kab) + Cos(fai) *

20、Sin(omg) * Cos(kab) c3 = Cos(fai) * Cos(omg) 像点坐标的近似值 For i = 1 To 5 x0(i) = -f * (a1 * (Xtp(i) - Xs0) + b1 * (Ytp(i) - Ys0) + c1 * (Ztp(i) - Zs0) / (a3 * (Xtp(i) - Xs0) + b3 * (Ytp(i) - Ys0) + c3 * (Ztp(i) - Zs0) y0(i) = -f * (a2 * (Xtp(i) - Xs0) + b2 * (Ytp(i) - Ys0) + c2 * (Ztp(i) - Zs0) / (a3 *

21、(Xtp(i) - Xs0) + b3 * (Ytp(i) - Ys0) + c3 * (Ztp(i) - Zs0) lx(i) = M(i) - x0(i) ly(i) = N(i) - y0(i) Next i 误差方程的系数 For i = 1 To 5 zba(i) = (a3 * (Xtp(i) - Xs0) + b3 * (Ytp(i) - Ys0) + c3 * (Ztp(i) - Zs0) P(i, 1, 1) = (a1 * f + a3 * M(i) / zba(i) P(i, 1, 2) = (b1 * f + b3 * M(i) / zba(i) P(i, 1, 3)

22、= (c1 * f + c3 * M(i) / zba(i) P(i, 1, 4) = N(i) * Sin(omg) - (M(i) * (M(i) * Cos(kab) - N(i) * Sin(kab) / f + f * Cos(kab) * Cos(omg) P(i, 1, 5) = -f * Sin(kab) - M(i) * (M(i) * Sin(kab) + N(i) * Cos(kab) / f P(i, 1, 6) = N(i) P(i, 2, 1) = (a2 * f + a3 * N(i) / zba(i) P(i, 2, 2) = (b2 * f + b3 * N(

23、i) / zba(i) P(i, 2, 3) = (c2 * f + c3 * N(i) / zba(i) P(i, 2, 4) = -M(i) * Sin(omg) - (N(i) * (M(i) * Cos(kab) - N(i) * Sin(kab) / f - f * Sin(kab) * Cos(omg) P(i, 2, 5) = -f * Cos(kab) - N(i) * (M(i) * Sin(kab) + N(i) * Cos(kab) / f P(i, 2, 6) = -M(i) Next i 误差方程常数项 For i = 1 To 6 A(1, i) = P(1, 1,

24、 i) A(2, i) = P(1, 2, i) A(3, i) = P(2, 1, i) A(4, i) = P(2, 2, i) A(5, i) = P(3, 1, i) A(6, i) = P(3, 2, i) A(7, i) = P(4, 1, i) A(8, i) = P(4, 2, i) A(9, i) = P(5, 1, i) A(10, i) = P(5, 2, i) Next i For i = 1 To 2 L(2 * i - 1, 1) = lx(i) L(2 * i, 1) = ly(i) Next i For i = 1 To 6 For j = 1 To 10 at

25、(i, j) = A(j, i) Next j Next i 求解法方程的解 JZC at(), A(), 6, 10, 6, BB() Rect BB(), 6, CC() JZC CC(), at(), 6, 6, 10, DD() JZC DD(), L(), 6, 10, 1, EE() If EE(1, 1) 1 And EE(2, 1) 1 And EE(3, 1) 1 And EE(4, 1) 0.00001 And EE(5, 1) 0.00001 And EE(6, 1) 0.00001 Then Exit Do Else Xs0 = Xs0 + EE(1, 1) Ys0 =

26、 Ys0 + EE(2, 1) Zs0 = Zs0 + EE(3, 1) fai = fai + EE(4, 1) omg = omg + EE(5, 1) kab = kab + EE(6, 1) lj = lj + 1 End IfLoopZs0 = w * f * 1000 + EZText1.Text = 外方位元素： & vbCrLf & Xs= & Format(Xs0, 0.) & vbCrLf & Ys= & Format(Ys0, 0.) & vbCrLf & Zs= & Format(Zs0, 0.) & vbCrLf & fai= & Format(fai, 0.) &

27、vbCrLf & omg= & Format(omg, 0.) & vbCrLf & kab= & Format(kab, 0.)Text1.Text = Text1 & vbCrLf & f=210. & vbCrLf & x0=0. & vbCrLf & y0=0. & vbCrLf & H= & Format(H, 0.)End Sub六、实验结果1504外方位元素： 1505外方位元素：Xs=. Xs=.Ys=. Ys=.Zs=911. Zs=910.=0. =0.=-0. =-0.=-0. =-0.f=210. f=210.x0=0. x0=0.y0=0. y0=0.H=7. H=7

28、.实验三 空间前方交会编程实验一、目的要求： 掌握空间前方交会的原理，根据所给立体像对两张像片的内、外方位元素，利用同名像点在左右像片上的坐标，解求其对应的地面点的物方坐标，利用计算机编程语言实现空间前方交会的过程，完成所给立体像对（1504、1505）上若干对同名点对应的地面物方点的坐标计算。 二、仪器用具： 计算机，编程软件（VB），通用图像处理软件photoshop。 三、实验数据 实验数据为一立体像对（1504、1505）两像片影像数据，以及它们的内、外方位元素，具体见表5-1。 15041505内方位元素h0-121.-121.h10.0.h20.0.k0115.115.k10.0.

29、k2-0.-0.外方位元素Xs.Ys.Zs910.911.0.0.-0.-0.-0.-0.表5-1 四、实验内容： 利用所给立体像对两张像片的内、外方位元素，编写空间前方交会程序，根据所给像对（1504、1505）中若干同名像点在左右像片上的坐标，解求其对应的地面点的物方坐标，实现空间前方交会的过程。 4.1数学模型： 基于点投影系数的空间前方交会式： 4.2 同名点量测： 利用photoshop图像处理软件将所给立体像对的两张像片影像打开，目视判断地物同名点，并记录下同名点在左右影像上的像素坐标（il，jl），（ir，jr），具体形式如表5-2。 点号 1505 1504 i（横） j（综）

30、 i（横） j（综） 1 7682 4862 4941 4857 2 7649 4895 4906 4890 3 7595 4947 4850 4941 4 7562 4981 4814 4974 5 7509 5034 4758 5026 6 7474 5068 4723 5058 7 7387 5154 4631 5143 表5-2 4.3 求解过程： 利用上面4.2中量测的同名像点的像素坐标以及对应像片的内定向参数，计算出同名点在左右像片的像框标坐标系中坐标。再根据同名点在左右像片的像平面坐标，和已知的左右像片外方位元素、以及4.1中的前方交会计算式，计算该像点对应地面点的坐标（X,Y,

31、Z），并计算地面点在Y方向上的残余上下视差Y。五、实验程序：读取同名点数据Private Sub read_Click() Dim i%CommonDialog1.ShowOpenOpen CommonDialog1.FileName For Input As #1For i = 1 To 10 Line Input #1, S If i 1 Then S = Split(S, ) P2(i - 1, 1) = S(0) Q2(i - 1, 1) = S(1) P1(i - 1, 1) = S(2) Q1(i - 1, 1) = S(3) End IfNext iText1.Text = 15

32、04 & & 1505 & vbCrLfFor i = 1 To 9 Text1.Text = Text1 & P1(i, 1) & & Q1(i, 1) & & P2(i, 1) & & Q2(i, 1) & vbCrLfNext iEnd Sub计算空间坐标Private Sub jisuan_Click()Dim Bx As Single, By As Single, Bz As SingleDim R1(3, 3) As Single, R2(3, 3) As SingleDim x1(9) As Single, y1(9) As Single, x2(9) As Single, y2

33、(9) As SingleDim a(9, 3, 1) As Single, b(9, 3, 1) As Single, c(9, 3, 1) As Single, d(9, 3, 1) As SingleDim M(9) As Single, N(9) As Single, R(9) As Single, S(9) As Single, T(9) As SingleDim k1(3, 1) As Single, k2(3, 1) As Single, k3(3, 1) As Single, k4(3, 1) As Single, k5(3, 1) As Single, k6(3, 1) As

34、 Single, k7(3, 1) As Single, k8(3, 1) As Single, k9(3, 1) As SingleDim l1(3, 1) As Single, l2(3, 1) As Single, l3(3, 1) As Single, l4(3, 1) As Single, l5(3, 1) As Single, l6(3, 1) As Single, l7(3, 1) As Single, l8(3, 1) As Single, l9(3, 1) As SingleDim g1(3, 1) As Single, g2(3, 1) As Single, g3(3, 1

35、) As Single, g4(3, 1) As Single, g5(3, 1) As Single, g6(3, 1) As Single, g7(3, 1) As Single, g8(3, 1) As Single, g9(3, 1) As SingleDim h1(3, 1) As Single, h2(3, 1) As Single, h3(3, 1) As Single, h4(3, 1) As Single, h5(3, 1) As Single, h6(3, 1) As Single, h7(3, 1) As Single, h8(3, 1) As Single, h9(3,

36、 1) As Single已知内外方位元素h01 = -121.h11 = 0.h21 = 0.k01 = 115.k11 = 0.k21 = -0.h02 = -121.h12 = 0.h22 = 0.k02 = 115.k12 = 0.0006k22 = -0.Xs1 = .Ys1 = .Zs1 = 911.5227fai1 = 0.omg1 = -0.01454kab1 = -0.Xs2 = .Ys2 = .Zs2 = 910.fai2 = 0.omg2 = -0.kab2 = -0.F = 210.681Bx = Xs2 - Xs1By = Ys2 - Ys1Bz = Zs2 - Zs

37、1正交矩阵R1，R2R1(1, 1) = Cos(fai1) * Cos(kab1) - Sin(fai1) * Sin(omg1) * Sin(kab1)R1(1, 2) = -Cos(fai1) * Sin(kab1) - Sin(fai1) * Sin(omg1) * Cos(kab1)R1(1, 3) = -Sin(fai1) * Cos(omg1)R1(2, 1) = Cos(omg1) * Sin(kab1)R1(2, 2) = Cos(omg1) * Cos(kab1)R1(2, 3) = -Sin(omg1)R1(3, 1) = Sin(fai1) * Cos(kab1) +

38、Cos(fai1) * Sin(omg1) * Sin(kab1)R1(3, 2) = -Sin(fai1) * Sin(kab1) + Cos(fai1) * Sin(omg1) * Cos(kab1)R1(3, 3) = Cos(fai1) * Cos(omg1)R2(1, 1) = Cos(fai2) * Cos(kab2) - Sin(fai2) * Sin(omg2) * Sin(kab2)R2(1, 2) = -Cos(fai2) * Sin(kab2) - Sin(fai2) * Sin(omg2) * Cos(kab2)R2(1, 3) = -Sin(fai2) * Cos(o

39、mg2)R2(2, 1) = Cos(omg2) * Sin(kab2)R2(2, 2) = Cos(omg2) * Cos(kab2)R2(2, 3) = -Sin(omg2)R2(3, 1) = Sin(fai2) * Cos(kab2) + Cos(fai2) * Sin(omg2) * Sin(kab2)R2(3, 2) = -Sin(fai2) * Sin(kab2) + Cos(fai2) * Sin(omg2) * Cos(kab2)R2(3, 3) = Cos(fai2) * Cos(omg2)像素坐标转换为像平面坐标For i = 1 To 9 x1(i) = h01 + P

40、1(i, 1) * h11 + Q1(i, 1) * h21 y1(i) = k01 + P1(i, 1) * k11 + Q1(i, 1) * k21 x2(i) = h02 + P2(i, 1) * h12 + Q2(i, 1) * h22 y2(i) = k02 + P2(i, 1) * k12 + Q2(i, 1) * k22Next iFor i = 1 To 9 a(i, 1, 1) = x1(i) a(i, 2, 1) = y1(i) a(i, 3, 1) = -F b(i, 1, 1) = x2(i) b(i, 2, 1) = y2(i) b(i, 3, 1) = -FNext iFor i = 1 To 3 k1(i, 1) = a(1, i, 1) k2(i, 1) = a(2, i, 1) k3(i, 1) = a(3, i, 1) k4(i, 1) = a(4, i, 1) k5(i, 1) = a(5, i, 1) k6(i, 1) = a(6, i, 1) k7(i, 1) = a(7, i, 1) k8(i, 1) = a(8, i, 1) k9(i, 1)