摄影测量实验报告

名目试验一数字内定向编程试验 2试验二空间前方交会编程试验 6试验三空间前方交会编程试验 16CAD成图试验 25六、试验成果29试验一数字内定向编程试验一、目的要求：像片影像的数字内定向。二、仪器用具：计算机，编程软件〔VB〕三、试验数据：1504150515041505ixjyixjy1417.229-110.000190.964109.999412.029-110.000198.964109.99928277.836109.999335.614110.0068271.936109.999367.414110.00638130.907109.9878195.164-109.9978101.307109.9878226.264-109.9974269.8648050.136241.0648058.036-109.999 -109.999 -109.999 -109.999-109.999 -109.999 -109.999 -109.999表3-1四、试验内容：内定向参数求取方法：利用航摄像片上的四个框标点的理论位置以及四个框标点的像素坐标为依据。依据条件：框标点的像素坐标通过公式转换后应等于其理论坐标。利用式〔3-1〕h0、h1、h2、k0、k1、计算公式：x=h0+h1*i+h2*jy=k0+k1*i+k2*j〔3-1〕求解过程：列误差方程法化求参数将像点扫描坐标转化为框标坐标通过公式：x=h0+h1*i+h2*jy=k0+k1*i+k2*j将测量的同名点的扫描坐标转换位像框标坐标系中坐标。五、试验程序：PrivateSubread_Click ”读入数据Dimi%,j%,k%CommonDialog1.ShowOpenOpenCommonDialog1.FileNameForInputAs#1Fori=1To4

k2。LineInput#1,ss=Split(s,““)Forj=1To2A(i,j+1)=s(j-1)NextjLineInput#1,ss=Split(s,““)Lx(i,1)=s(0)Ly(i,1)=s(1)NextiFori=1To4A(i,1)=1NextiEndSub”内定向求解PrivateSubjisuan_ClickDimi%,j%,k%Fori=1To3Forj=1To4at(i,j)=A(j,i)NextjNexti”矩阵求解JZCat,A,3,4,3,BRectB,3,BdNJZCBdN,at,3,3,4,bbJZCbb,Lx,3,4,1,HHJZCbb,Ly,3,4,1,KKFori=1To3Text1=Text1&Format(HH(i,1),“0.000000“)&vbCrLfNextiText1=Text1&vbCrLfFori=1To3Text1=Text1&Format(KK(i,1),“0.000000“)&vbCrLfNextiClose#1EndSub”坐标转换PrivateSubCommand1_Clickxx=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“)EndSub六、试验结果：1504定向参数：=0.006349y方向：k0=115.125954,k1=0.000516,k2=-0.027982,精度=0.0013381505定向参数：=0.005566y方向：k0=115.321828,k1=0.000600,k2=-0.027981,精度=0.002853试验二空间前方交会编程试验一、目的要求：〔15041505〕中两张像片各自的6个外方位元素的解求。二、仪器用具：计算机，编程软件〔VB〕三、试验数据左片〔1504〕试验数据为一立体像对〔1504、1505〕的像片掌握点数据，包括像点坐标及相应地面摄测坐标。具体如表4-1所示：左片〔1504〕像平面坐标像平面坐标地面摄测坐标xyXtpYtpZtp-2.81673.965501286.070543471.38014.250-6.459-87.783501261.140542778.3305.580-76.415-46.343500966.380542964.9805.430-30.041-40.404501163.290542986.8008.810-65.89075.337501019.750543480.2305.760xyXtpYtpZtp72.50774.412501286.070543471.38014.25073.384-85.905501261.140542778.3305.5803.767-47.117500966.380542964.9805.43049.064-39.709501163.290542986.8008.81010.24473.307501019.750543480.2305.760右片〔1505〕像平面坐标右片〔1505〕像平面坐标地面摄测坐标两张像片的内方位元素为：x0=y0=0mm,f=210.681mm。四、试验内容：利用地面掌握点数据以及相应像点坐标〔1504、1505〕中两张像片各自的6个外方位元素。数学模型：共线条件方程式xfyf

a3(XXs)b3(YYs)c3(ZZs)a2(XXs)b2(YYs)c2(ZZs)a3(XXs)b3(YYs)c3(ZZs)求解过程：猎取数据：像片比例尺1/m：由像片上两点间距离与相应地面点间距之比求得。平均航高H：H=m*f+ΣZi/n其中，n为掌握点数。内方位元素x0,y0,f：地面摄影测量坐标Xtp,Ytp,Ztp掌握点的像点坐标〔表4-1〕确定未知数初始值：在竖直摄影状况下，角元素的初始值为0，即φ=ω=κ=0；线元素中，Zs0=H,Xs0,Ys0的取值可用五个掌握点坐标的平均值，即：Xs0=ΣXtpi/5，Ys0=ΣYtpi/5。计算旋转矩阵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(κ)c2=-Sin(φ)*Sin(κ)+Cos(φ)*Sin(ω)*Cos(κ)c3=Cos(φ)*Cos(ω)逐点计算像点坐标的近似值：利用未知数的近似值按共线方程式计算掌握点像点坐标的近似值〔x〕，〔y〕组成误差方程式：按下式逐点计算误差方程式的系数和常数项：系数：zba=(a3*(Xa-Xs)+b3*(Ya-Ys)+c3*(Za-zs))a11=(a1*f+a3*x)/zbaa12=(b1*f+b3*x)/zbaa13=(c1*f+c3*x)/zbaa21=(a2*f+a3*y)/zbaa22=(b2*f+b3*y)/zbaa23=(c2*f+c3*y)/zbaa14=y*Sin(ω)-(x*(x*Cos(κ)-y*Sin(κ))/f+f*Cos(κ))*Cos(ω)a15=-f*Sin(κ)-x*(x*Sin(κ)+y*Cos(κ))/fa16=ya24=-x*Sin(ω)-(y*(x*Cos(κ)-y*Sin(κ))/f-f*Sin(κ))*Cos(ω)a25=-f*Cos(κ)-y*(x*Sin(κ)+y*Cos(κ))/fa26=-x常数项：lx=x-〔x〕ly=y-〔y〕用矩阵形式表示为：V=AX-l式中：V=[VxVy]A=a11a12a13a14a15a16a21a22a23a24a25a26X=[dXsdYsdZsdφdωdκ]L=[lxly]’全部掌握点，列出误差方程式，构成总误差方程式：V=AX-l式中：V=[V1V2V3V4V5]A=[A1A2A3A4A5]L=[l1l2l3l4l5]组成法方程式：计算法方程的系数矩阵Na=A’A与常数项U=A’l。解求外方位元素：依据法方程，按下式解求外方位元素改正数：X=Nani*U并与相应的近似值求和，得到外方位元素的近似值。检查计算是否收敛：重复4-8的计算，直到满足要求为止。迭代限差为：dXs<1m,dYs<1m,dZs<1m,dφ<0.00001,dω<0.00001,dκ<0.00001最终得出的六个外方位元素的解为：Xs=Xs0+dXs1+dXs2+……Ys=Ys0+dYs1+dYs2+……Zs=Zs0+dZs1+dZs2+……φ=φ0+dφ1+dφ2+……ω=ω0+dω1+dω2+……κ=κ0+dκ1+dκ2+……五、试验程序：”读取数据PrivateSubread_ClickDimi%CommonDialog1.ShowOpenOpenCommonDialog1.FileNameForInputAs#1Fori=1To6LineInput#1,sIfi<>1Thens=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)EndIfNextiEndSub”计算PrivateSubjisuan_ClickDimi%,j%Dims1AsSingle,s2AsSingle,wAsSingle,EZAsSingle,HAsSingle,Xs0AsSingle,Ys0AsSingle,Zs0AsSingle,faiAsSingle,omgAsSingle,kabAsSingleDima1AsSingle,a2AsSingle,a3AsSingle,b1AsSingle,b2AsSingle,b3AsSingle,c1AsSingle,c2AsSingle,c3AsSingleDimA(10,6)AsSingle,L(10,1)AsSingle,P(5,2,6)AsSingleDimat(6,10)AsSingle,BB(6,6)AsSingle,CC(6,6)AsSingle,DD(6,10)AsSingle,EE(6,1)AsSingleDimx0(5)AsSingle,y0(5)AsSingle,lx(5)AsSingle,ly(5)AsSingle,zba(5)AsSinglef=0.210681”像片比例尺1/ms1=Sqr((M(1)-M(2))^2+(N(1)-N(2))^2)s2=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=1DoWhilelj<6”R矩阵系数a1=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)*Sin(omg)*Cos(kab)c3=Cos(fai)*Cos(omg)”像点坐标的近似值Fori=1To5x0(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*(Xtp(i)-Xs0)+b3*(Ytp(i)-Ys0)+c3*(Ztp(i)-Zs0))lx(i)=M(i)-x0(i)ly(i)=N(i)-y0(i)Nexti”误差方程的系数Fori=1To5zba(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)=(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))/fP(i,1,6)=N(i)P(i,2,1)=(a2*f+a3*N(i))/zba(i)P(i,2,2)=(b2*f+b3*N(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))/fP(i,2,6)=-M(i)Nexti”误差方程常数项Fori=1To6A(1,i)=P(1,1,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)NextiFori=1To2L(2*i-1,1)=lx(i)L(2*i,1)=ly(i)NextiFori=1To6Forj=1To10at(i,j)=A(j,i)NextjNexti”求解法方程的解JZCat,A,6,10,6,BBRectBB,6,CCJZCCC,at,6,6,10,DDJZCDD,L,6,10,1,EEIfEE(1,1)<1AndEE(2,1)<1AndEE(3,1)<1AndEE(4,1)<0.00001AndEE(5,1)<0.00001AndEE(6,1)<0.00001ThenExitDoElseXs0=Xs0+EE(1,1)Ys0=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+1EndIfLoopZs0=w*f*1000+EZText1.Text&vbCrLf&“Xs=“&Format(Xs0,“0.000000“)&vbCrLf&“Ys=“&Format(Ys0,“0.000000“)&vbCrLf&“Zs=“&Format(Zs0,“0.000000“)&vbCrLf&“fai=“&Format(fai,“0.000000“)&vbCrLf&“omg=“&Format(omg,“0.000000“)&vbCrLf&“kab=“&Format(kab,“0.000000“)Text1.Text=Text1&vbCrLf&“f=210.681000“&vbCrLf&“x0=0.000000“&vbCrLf&“y0=0.000000“&vbCrLf&“H=“&Format(H,“0.000000“)EndSub六、试验结果1504外方位元素：Xs=501257.416158Ys=543170.797562Zs=911.522700φ=0.040643ω=-0.014540κ=-0.013461x0=0.000000y0=0.000000H=7.966000

1505外方位元素：Xs=500934.463445Ys=543180.119596Zs=910.372341φ=0.028448ω=-0.012638κ=-0.050018x0=0.000000y0=0.000000H=7.966000试验三空间前方交会编程试验一、目的要求：会的过程，完成所给立体像对〔1504、1505〕上假设干对同名点对应的地面物方点的坐标计算。二、仪器用具：计算机，编程软件〔VB〕，通用图像处理软件photoshop。三、试验数据试验数据为一立体像对〔1504、1505〕两像片影像数据，以及它们的内、外方位元素，具体见表5-1。15041505h0-121.642023-121.766062h10.0279760.027976内方位元素h2k00.000607115.3218280.000523115.125954k10.0006000.000516k2-0.027981-0.027982Xs500934.463445501257.416158Ys543180.119596543170.797562外方位元素Zsφ910.3723410.028448911.5227000.040643ω-0.012638-0.014540κ-0.050018-0.013461表5-1四、试验内容：利用所给立体像对两张像片的内、外方位元素，编写空间前方交会程序，依据所给像对〔1504、1505〕中假设干同名像点在左右像片上的坐标，解求其对应的地面点的物方坐标，实现空间前方交会的过程。数学模型：基于点投影系数的空间前方交会式：B X XB X XBs2 s1BB Z B Z ZX x 2 2 2 2 2a (X ,Y ,Z )2 2 2 2 2

2 R2yRZ2

fX x 1 1 1 1 1a(X ,Y ,Z )1 1 1 1 1

1 R1yRZ1

fZXZ1N BX 2ZXZ1

BX Z2Z21 2 2Z1Z1N BZ1

BXZXZ21 2 2Z1ZXZ211UNX11

U”NX22”22122VNY VNY12222WN1Z1 W”NZ22XX NXS1 1 1S1 1 1 S1 2 1S1 1 1 S1 2 Y 2

(Y NY)(Y NY)ZZS1N1Z1dY(Y NY)(Y NY)S1 1 1 S2 2 2同名点量测：利用photoshop并记录下同名点在左右影像上的像素坐标〔i，j〕，〔i，j〕，具体形式如表5-2。l l r r点号点号15051504i〔横〕 j〔综〕i〔横〕j〔综〕17682486249414857276494895490648903759549474850494147562498148144974575095034475850266747450684723505877387515446315143求解过程：

表5-2利用上面4.2中量测的同名像点的像素坐标以及对应像片的内定向参数，计算出同名点在左外方位元素、以及4.1中的前方交会计算式，计算该像点对应地面点的坐标〔X,Y,Z〕，并计算地面点在Y方向上的剩余上下视差ΔY。五、试验程序：”读取同名点数据PrivateSubread_ClickDimi%CommonDialog1.ShowOpenOpenCommonDialog1.FileNameForInputAs#1Fori=1To10LineInput#1,SIfi<>1ThenS=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)EndIfNextiText1.Text=“1504“&“ “&“1505“&vbCrLfFori=1To9Text1.Text=Text1&P1(i,1)&““&Q1(i,1)&“ “&P2(i,1)&““&Q2(i,1)&vbCrLfNextiEndSub”计算空间坐标PrivateSubjisuan_ClickDimBxAsSingle,ByAsSingle,BzAsSingleDimR1(3,3)AsSingle,R2(3,3)AsSingleDimx1(9)AsSingle,y1(9)AsSingle,x2(9)AsSingle,y2(9)AsSingleDima(9,3,1)AsSingle,b(9,3,1)AsSingle,c(9,3,1)AsSingle,d(9,3,1)AsSingleDimM(9)AsSingle,N(9)AsSingle,R(9)AsSingle,S(9)AsSingle,T(9)AsSingleDimk1(3,1)AsSingle,k2(3,1)AsSingle,k3(3,1)AsSingle,k4(3,1)AsSingle,k5(3,1)AsSingle,k6(3,1)AsSingle,k7(3,1)AsSingle,k8(3,1)AsSingle,k9(3,1)AsSingleDiml1(3,1)AsSingle,l2(3,1)AsSingle,l3(3,1)AsSingle,l4(3,1)AsSingle,l5(3,1)AsSingle,l6(3,1)AsSingle,l7(3,1)AsSingle,l8(3,1)AsSingle,l9(3,1)AsSingleDimg1(3,1)AsSingle,g2(3,1)AsSingle,g3(3,1)AsSingle,g4(3,1)AsSingle,g5(3,1)AsSingle,g6(3,1)AsSingle,g7(3,1)AsSingle,g8(3,1)AsSingle,g9(3,1)AsSingleDimh1(3,1)AsSingle,h2(3,1)AsSingle,h3(3,1)AsSingle,h4(3,1)AsSingle,h5(3,1)AsSingle,h6(3,1)AsSingle,h7(3,1)AsSingle,h8(3,1)AsSingle,h9(3,1)AsSingle”内外方位元素h01=-121.766062h11=0.027976h21=0.000523k01=115.125954k11=0.000516k21=-0.027982h02=-121.642023h12=0.027976h22=0.000607k02=115.321828k12=