近景摄影测量实验报告

课间实验报告2015年——2016年第2学期实验课程：近景摄影测量实验班级：学生姓名：学号：指导教师：重庆交通大学测量与空间数据处理实验室目录实验一近景单张像片空间后方交会.......................................3实验二近景立体像对空间前方交会.....................................8实验三直接线性变换求L系数精确解.................................13实验四直接线性变换求空间坐标精确值.............................17

实验一近景单张像片空间后方交会一、实验目的通过单张像片空间后方交会算法编程，掌握空间后方交会的基本原理和基本步骤，为近景摄影测量解析处理方法打下理论基础。二、实验内容利用VB编程平台，通过给定的单片像点的像点坐标、内方位元素和控制点物方空间坐标，计算出像片的外方位元素。三、实验步骤：1、编程流程图开始获取已知数据：比例尺，地面点坐标，像点坐标和内方位元素定义旋转矩阵R组成矩阵A和L计算内方位元素改正值判断内方位元素改正值是否收敛否是输出外方位元素结束2、程序代码'定义内方位元素变量Dimx0,yo,fAsDouble'定义外方位元素变量DimXs,Ys,ZsAsDoubleDimQ,w,kAsDouble'定义旋转变量元素Dima1,a2,a3AsDoubleDimb1,b2,b3AsDoubleDimc1,c2,c3AsDouble'定义A矩阵Dima11,a12,a13,a14,a15,a16AsDoubleDima21,a22,a23,a24,a25,a26AsDouble'定义像片共线方程变量DimX8AsDouble,Y8AsDouble,Z8AsDouble'定义像点坐标变量Dimx,y,zAsDouble'定义像片L变量DimL1AsDouble,L2AsDouble'定义外方位角元素三角函数值Dimsk,sw,ck,cwAsDoublePrivateSubCommand1_Click()x0=48.6401Txtx0=x0y0=11.4765Txty0=y0f=4548.5949Txtf=fXs=3391.658TxtXs=XsYs=-135.645

TxtYs=YsZs=97.25TxtZs=ZsTxtQ="-0°35′9.00″"Txtw="3°33′13.00″"Txtk="-5°4′59.00″"EndSubPrivateSubCommand3_Click()'计算外方位角元素的三角函数值sw=0.061982457cw=0.998077239sk=-0.088599727ck=0.996067311'输入像点坐标x=599.331y=1894.198z=-4929.805'输入旋转矩阵a1=0.995959:a2=0.089231:a3=0.010205b1=-0.088434:b2=0.994152:b3=-0.061983c1=-0.015676:c2=0.06083:c3=0.998025'计算像片共线方程系数X8=a1*(x-Xs)+b1*(y-Ys)+c1*(z-Zs)Y8=a2*(x-Xs)+b2*(y-Ys)+c2*(z-Zs)Z8=a3*(x-Xs)+b3*(y-Ys)+c3*(z-Zs)'计算矩阵Aa11=(a1*f+a3*(x-x0))/Z8a12=(b1*f+b3*(x-x0))/Z8a13=(c1*f+c3*(x-x0))/Z8a14=(y-y0)*sw-((x-x0)*((x-x0)*ck-(y-y0)*sk)/f+f*ck)*cwa15=-f*sk-(x-x0)*((x-x0)*sk+(y-y0)*ck)/fa16=y-y0a21=(a2*f+a3*(y-y0))/Z8a22=(b2*f+b3*(y-y0))/Z8a23=(c2*f+c3*(y-y0))/Z8a24=-(x-x0)*sw-((y-y0)*((x-x0)*ck-(y-y0)*sk)/f-f*sk)*

cwa25=-f*ck-(y-y0)*((x-x0)*sk+(y-y0)*ck)/fa26=-(x-x0)'计算LL1=x-(x0-f*X8/Z8)L2=y-(y0-f*Y8/Z8)Printa11,a12,a13,a14,a15,a16Printa21,a22,a23,a24,a25,a26PrintL1,L2Xs=3391.658-0.768Ys=-135.645-0.016Zs=97.25-0.819Q="-0°35′8.00″"w="3°33′9.00″"Txtk="-5°4′59.00″"MsgBox"改正后的Xs为："&Xs&Chr(13)&"改正后的Ys为："&Ys&Chr(13)&"改正后的Zs为："&Zs&Chr(13)&"改正后的φ为："&Q&Chr(13)&"改正后的ω为："&w&Chr(13)&"改正后的κ为："&Txtk,vbOKOnly,"计算结果"EndSub四、实验结果1、实验数据物方坐标：像点坐标：以毫米为单位点名，x,第一列是点号；单位：毫米第二列是X坐标，朝南为正向；第三列是Y坐标，朝上为正向；y第四列为Z坐标，朝前(西)为正向。43011617.7121883.870-4917.30843021615.5721687.514-4913.75043031614.6401315.066-4916.72443041615.7841001.358-4916.69143051617.197690.040-4917.07043061615.933387.492-4917.23443071618.946145.491-4917.80743081616.753-186.640-4918.53143091617.714-482.264-4918.56143101617.859-680.767-4918.7994301,1105.312313,3367.9792424302,1111.483438,3203.9980094303,1128.936009,2886.1463744304,1145.729428,2615.3853014305,1163.579521,2343.6035844306,1179.483899,2076.6727384307,1196.664025,1861.9718854308,1216.048311,1565.0460874309,1236.196347,1299.2754134310,1249.998721,1120.52425x0=48.6401mm，y0=11.4765mm，f=4548.5949mm

五、心得体会在做此次试验之前，先对上学期所学的该知识点进行了系统的复习，从整体上把握了单相空间后方交会原理及解算过程。也理解了将控制点的地面点坐标视为真值，把相应的像点坐标视为观测值，加入相应的改正数，最后列出每个点的误差方程，然后用矩阵表示，最后构成总误差方程，从而求出外方位元素的改正数。实验二近景立体像对空间前方交会一、实验目的通过立体像对空间前方交会算法编程，掌握空间前方交会的基本原理和基本步骤，为近景摄影测量解析处理方法打下理论基础。二、实验内容利用VB编程平台，通过给定的立体像对像点的像点坐标以及对应的内外方位元素，实现待测控制点的物方空间坐标的求解。三、实验步骤：1、编程流程图开始获取数据将获取数据强制转换为double型计算旋转矩阵中的元素计算u1，v1,w1和u2，v2，w2计算投影系数N1和N2解求地面点坐标将结果导出至指定窗体中结束2、程序代码Dimf#,x1#,y1#,x2#,y2#,i%,j%,u1#,u2#,v1#,v2#,w1#,w2#,fai1#,kab1#,omg1#,fai2#,kab2#,omg2#Dima12#,a13#,b11#,b12#,b13#,c11#,c12#,c13#,a21#,a22#,a23#,b21#,b22#,b23#,c21#,c22#,c23#Dimn1#,n2#,bu#,bv#,bw#Dimxs1#,xs2#,ys1#,ys2#,zs1#,zs2#PrivateSubCommand1_Click()fai1=Val(Text1.Text):kab1=Val(Text3.Text):omg1=Val(Text2.Text)fai2=Val(Text11.Text):kab2=Val(Text10.Text):omg2=Val(Text7.Text)f=Val(Text13.Text)xs1=Val(Text4.Text):xs2=Val(Text9.Text):ys1=Val(Text5.Text):ys2=Val(Text8.Text):zs1=Val(Text6.Text):zs2=Val(Text12.Text)EndSubPrivateSubCommand3_Click()x1=Val(InputBox("输入第一片坐标x的值"))y1=Val(InputBox("输入第一片坐标y的值"))x2=Val(InputBox("输入第二片坐标x的值"))y2=Val(InputBox("输入第二片坐标y的值"))EndSubPrivateSubCommand2_Click()Text1.Visible=FalseText2.Visible=FalseText3.Visible=FalseText4.Visible=FalseText5.Visible=FalseText6.Visible=FalseText7.Visible=FalseText8.Visible=FalseText9.Visible=FalseText10.Visible=FalseText11.Visible=FalseText12.Visible=FalseText13.Visible=FalseLabel1.Visible=FalseLabel2.Visible=FalseLabel3.Visible=FalseLabel4.Visible=FalseLabel5.Visible=FalseLabel6.Visible=FalseLabel7.Visible=FalseLabel8.Visible=FalseLabel9.Visible=FalseLabel10.Visible=False

Label11.Visible=FalseLabel12.Visible=FalseLabel13.Visible=FalseLabel14.Visible=Falsea11=Cos(fai1)*Cos(kab1)-Sin(fai1)*Sin(omg1)*Sin(kab1)a12=-Cos(fai1)*Sin(kab1)-Sin(fai1)*Sin(omg1)*Cos(kab1)a13=-Sin(fai1)*Cos(omg1)b11=Cos(omg1)*Sin(kab1)b12=Cos(omg1)*Cos(kab1)b13=-Sin(omg1)c11=Sin(fai1)*Cos(kab1)+Cos(fai1)*Sin(omg1)*Sin(kab1)c12=-Sin(fai1)*Sin(kab1)+Cos(fai1)*Sin(omg1)*Cos(kab1)c13=Cos(fai1)*Cos(omg1)a21=Cos(fai2)*Cos(kab2)-Sin(fai2)*Sin(omg2)*Sin(kab2)a22=-Cos(fai2)*Sin(kab2)-Sin(fai2)*Sin(omg2)*Cos(kab2)a23=-Sin(fai2)*Cos(omg2)b21=Cos(omg2)*Sin(kab2)b22=Cos(omg2)*Cos(kab2)b23=-Sin(omg2)c21=Sin(fai2)*Cos(kab2)+Cos(fai2)*Sin(omg2)*Sin(kab2)c22=-Sin(fai2)*Sin(kab2)+Cos(fai2)*Sin(omg2)*Cos(kab2)c23=Cos(fai2)*Cos(omg2)u1=a11*x1+a12*y1-a13*fv1=b11*x1+b12*y1-b13*fw1=c11*x1+c12*y1-c13*fu2=a21*x2+a22*y2-a23*fv2=b21*x2+b22*y2-b23*fw2=c21*x2+c22*y2-c23*fbu=xs2-xs1bv=ys2-ys1bw=zs2-zs1n1=(bu*w2-bw*u2)/(u1*w2-u2*w1)n2=(bu*w1-bw*u1)/(u1*w2-u2*w1)Dimxx#,yy#,yy1#,yy2#,zz#xx=Fix((xs1+u1*n1)*1000+0.5)/1000yy1=Val(Text5.Text)+v1*n1yy2=Val(Text8.Text)+v2*n2yy=Fix(((yy1+yy2)/2)*1000+0.5)/1000zz=Fix((Val(Text6.Text)+w1*n1)*1000+0.5)/1000

Form1.Print"地面坐标为:"Form1.Print"xx=";xx,"yy=";yy;"zz=";zzEndSub四、实验结果1、实验数据相片1坐标：相片2坐标：点名，x,y点名，x,y4201,896.942591,3200.1324684202,903.373104,3012.7248834203,918.450557,2802.550964204,937.004622,2449.9082384205,957.037515,2086.3138984206,970.889881,1858.4535284207,979.044342,1714.8198454208,989.950247,1549.5977224209,1001.556951,1384.7020554210,1015.53176,1189.9886894201,595.600431,3258.0235544202,601.34579,3061.8024614203,615.818572,2841.9570754204,633.003988,2472.9369064205,651.907837,2092.444594206,664.95803,1854.1427144207,672.855606,1703.7227394208,683.251761,1530.7373154209,694.71749,1357.8616454210,708.184825,1154.038955第1张像片：Xs=4310.470mmYs=-145.342mmZs=-295.236mmR=0.9496930.1011730.296391-0.0865340.994309-0.062137-0.3009910.0333630.953043φ=ω=κ=-4°58′25.00″第2张像片：Xs=4079.418mmYs=-145.179mmZs=-298.330mmR=0.9730900.0961040.209429-0.0851000.994506-0.060955-0.2141360.0414920.975922φ=-12°6′42.00″ω=3°29′40.00″κ=-4°53′27.00″x0=48.6401mmy0=11.4765mm，f=4548.5949mm

五、心得体会通过这次编程实验，对于共线方程和立体像对前方交会公式的理解更为透彻了，对于方程中Bu矩阵和常数L矩阵的求解也更为熟悉。各个等式中每个字母所代表的含义也更为清楚。前方交会的编程其实不算难，只是在计算中需要定义的变量很多，很容易混淆，因此在定义变量时要加上注释并定义得与变量所代表的含义相近。由于在VB中计算矩阵比较困难，因此可以将矩阵中的每个元素定义为一个相应的变量，在VB计算出后再到Matlab中计算矩阵求出改正数，将所求出的改正数在VB中输入并与之前的近似坐标计算再输出最后物方坐标。这样求解比在VB中计算矩阵要省时省力。所以在实验中，不妨多借助几样软件来达到高效运算的效果。实验三直接线性变换求L系数精确解一、实验目的通过直接线性变换基本原理编程解算L系数精确值，掌握直接线性变换的基本原理和基本步骤，为近景摄影测量解析处理方法打下理论基础。二、实验内容利用MATLAB编程平台，通过给定的单片像点的像点坐标和控制点物方空间坐标，计算出像片的L系数精确值。三、实验步骤：1、编程流程图开始输入点位的像点坐标和物方坐标按照所给点位的坐输入像点矩标计算L近似值公阵B式中的矩阵A按照公式计算L的近似值计算基本关系式分计算像点像径r母C计算矩阵M计算矩阵W计算L的精确值结束2、编程代码①矩阵AA=1.0e+007*0.00020.0002-0.00050.000000000.17880.2082-0.543500000.00020.0002-0.00050.00000.54480.6345-1.65610.00020.0002-0.00050.000000000.17960.1876-0.546200000.00020.0002-0.00050.00000.51760.5407-1.57440.00020.0001-0.00050.000000000.18230.1485-0.555100000.00020.0001-0.00050.00000.46600.3795-1.41900.00020.0001-0.00050.000000000.18510.1147-0.563300000.00020.0001-0.00050.00000.42260.2619-1.28590.00020.0001-0.00050.000000000.18820.0803-0.572100000.00020.0001-0.00050.00000.37900.1617-1.15240.00020.0000-0.00050.000000000.19060.0457-0.5800②矩阵BB=1.0e+003*1.10533.36801.11153.20401.12892.88611.14572.61541.16362.34361.1795③L的近似值L=1.0e+003*-0.00010.0000-0.0002-1.1128-0.00010.0000-0.0007-3.22460.0000-0.00000.0000

④矩阵MM=1.0e+011*-0.0000-0.00000.0001-0.00000000-0.0201-0.02350.06120.13080000-0.0000-0.00000.0001-0.0000-0.0614-0.00000.18650.4156-0.0000-0.00010.0004-0.00000000-0.0152-0.01590.04630.12030000-0.0001-0.00010.0004-0.0000-0.43840.45791.33340.3615-0.0000-0.00000.0001-0.00000000-0.0242-0.01970.07370.10190000-0.0000-0.00000.0001-0.0000-0.0619-0.05040.26510.2711-0.0000-0.00000.0001-0.00000000-0.0301-0.01870.09170.08760000-0.0000-0.00000.0001-0.0000-0.0688-0.04260.20940.2079-0.0000-0.00000.0001-0.00000000-0.0340-0.01450.10330.07450000-0.0000-0.00000.0001-0.0000-0.0684-0.02920.20810.1558-0.0000-0.00000.0001-0.00000000-0.0335-0.00800.10190.06270000-0.0000-0.00000.0001-0.0000-0.0590-0.01410.17940.1145⑤矩阵WW=1.0e+007*0.12450.37930.94142.71370.14990.38330.18650.42580.21010.42320.2073

五、实验结果分析1、实验数据像片编号2-40-g11-k1-0012201,1358.55679,2940.5965532202,1375.104882,2714.0732622203,1396.439124,2359.6673272206,1454.452586,1516.0217162207,1467.812313,1322.0091122209,1493.982413,993.473865控制点空间坐标2201597.7471785.157-5930.7482202600.3771446.575-5934.1462203599.126922.699-5934.2242206599.171-299.231-5933.7772207597.520-577.132-5933.3442209598.152-1046.647-5934.4452、实验结果LG=1.0e+003*0.00070.0000-0.0008-5.16810.00550.00000.0004-6.6576-0.00000.0000-0.00000.0000五、实验心得通过本次实验我的体会很多，就比如在求解L的近似值，过程比较繁琐，且求解物方空间坐标的精确值需要和法方程联立才能得到结果。直接线性变换虽然是直接建立了像点坐标和物方空间坐标的关系，但其编程求解过程很是繁琐。矩阵的运算量非常的大。相比之下倒是共线方程更为快捷。

实验四直接线性变换求空间坐标精确值一、实验目的通过直接线性变换算法编程，掌握直接线性变换的基本原理和基本步骤，为近景摄影测量解析处理方法打下理论基础。二、实验内容利用VB编程平台，通过给定的立体像对的同名像点坐标和L系数，解算待测点的物方空间坐标。三、实验步骤：1、编程流程图利用L矩阵求内方位元素：开始输入物方空间坐标原点在像空间坐标系中的坐标输入的精确值Lx0=-(L1*L9+L2*L10+L3*L11)/(L9^2+L10^2+L11^2)按照以上公式计算x0y0=-(L5*L9+L6*L10+L7*L11)/(L9^2+L10^2+L11^2)按照以上公式计算y0计算ds=E计算dβ=arcsinDfx=cosdβ*[(L1^2*L2^2*L3^2)/(L9^2+L10^2+L11^2)-x0^2]^1/2fy=fx/(1+ds)输出内方位元素结束利用L矩阵求外方位元素：开始输入物方空间坐标原点在像空间输入L的精确值坐标系中的坐标按公式计算外方位元素的线元素按公式计算外方位元素的角元素输出内方位元素结束利用立体像对像点坐标和内外方位元素求物方坐标：开始获取数据将获取数据强制转换为double型计算旋转矩阵中的元素计算u1，v1,w1和u2，v2，w2计算投影系数N1和N2解求地面点坐标将结果导出至指定窗体中结束2、编程代码(1)用VB编程计算内方位元素和外方位元素角元素：'定义L的精确值变量DimL1,L2,L3,L4,L5,L6,L7,L8,L9,L10,L11AsDouble'定义r3Dimr3AsDouble'定义计算变量DimA,B,CAsDouble'定义内方位元素Dimx0,y0AsDoubleDimfx,fyAsDoubleDimdSAsDoubleDimdBAsString'定义外方位角元素DimQ,W,KAsDouble'定义计算变量Dima3,c3,b3AsDoubleDimb1,b2AsDoublePrivateSubCommand1_Click()'输入L的精确值L1=0.691694694954L2=0.000608383183L3=-0.825372309722L4=-5168.07051610071L5=5.470559066537L6=0.008254963815L7=0.441376855385L8=-6657.62884269957L9=-0.000695566593L10=0.000000547457L11=-0.000023826626'计算r3r3=Sqr(L9^2+L10^2+L11^2)'计算x0,y0x0=-(L1*L9+L2*L10+L3*L11)/(L9^2+L10^2+L11^2)

A=r3^2*(L1^2+L2^2+L3^2)-x0^2B=r3^2*(L5^2+L6^2+L7^2)-y0^2C=r3^2*(L1*L5+L2*L6+L3*L7)-x0*y0vbOKOnly,"内方位元素计算结果"外方位角元素计算结果"EndSub(3)用VB编程计算物方空间坐标：'定义两个像点坐标Dimx1,y1AsDoubleDimx2,y2AsDouble'定义物方空间坐标DimX,Y,ZAsDouble'定义两张像片的L系数DimL1,L2,L3,L4,L5,L6,L7,L8,L9,L10,L11AsDoubleDimLL1,LL2,LL3,LL4,LL5,LL6,LL7,LL8,LL9,LL10,LL11AsDouble'定义基本关系式分母DimA1,A2,A3AsDouble'定义内方位