摄影测量学空间后方交会实验报告

摄影测量学实验报告实验一、单像空间后方交会学院：建测学院班级：测绘082姓名：肖澎学号：15一．实验目的1.深入了解单像空间后方交会的计算过程；2.加强空间后方交会基本公式和误差方程式，法线方程式的记忆；3.通过上机调试程序加强动手能力的培养。二．实验原理以单幅影像为基础，从该影像所覆盖地面范围内若干控制点和相应点的像坐标量测值出发，根据共线条件方程，求解该影像在航空摄影时刻的相片外方位元素。三．实验内容1.程序图框图2.实验数据（1）已知航摄仪内方位元素f＝153.24mm，Xo＝Yo＝0。限差0.1秒（2）已知4对点的影像坐标和地面坐标：影像坐标地面坐标x（mm）y（mm）X(m)Y（m）Z(m)1-86.15-68.9936589.4125273.322195.172-53.4082.2137631.0831324.51728.693-14.78-76.6339100.9724934.982386.50410.4664.4340426.5430319.81757.313.实验程序usingSystem;usingSystem.Collections.Generic;usingSystem.Linq;usingSystem.Text;namespaceConsoleApplication3{classProgram{staticvoidMain(){//输入比例尺，主距，参与平参点的个数Console.WriteLine(请输入比例尺分母m：\r);stringm1=Console.ReadLine();doublem=(double)Convert.ToSingle(m1);Console.WriteLine(请输入主距f：\r);stringf1=Console.ReadLine();doublef=(double)Convert.ToSingle(f1);Console.WriteLine(请输入参与平差控制点的个数n：\r);stringn1=Console.ReadLine();intn=(int)Convert.ToSingle(n1);//像点坐标的输入代码double[]arr1=newdouble[2*n];//1.像点x坐标的输入for(inti=0;in;i++){Console.WriteLine(请输入已进行系统误差改正的像点坐标的x{0}值：\r,i+1);stringu=Console.ReadLine();for(intj=0;jn;j+=2){arr1[j]=(double)Convert.ToSingle(u);}}//2.像点y坐标的输入for(inti=0;in;i++){Console.WriteLine(请输入已进行系统误差改正的像点坐标的y{0}值：\r,i+1);stringv=Console.ReadLine();for(intj=1;jn;j+=2){arr1[j]=(double)Convert.ToSingle(v);}}//控制点的坐标输入代码double[,]arr2=newdouble[n,3];//1.控制点X坐标的输入for(intj=0;jn;j++){Console.WriteLine(请输入控制点在地面摄影测量坐标系的坐标的X{0}值：\r,j+1);stringu=Console.ReadLine();arr2[j,0]=(double)Convert.ToSingle(u);}//2.控制点Y坐标的输入for(intk=0;kn;k++){Console.WriteLine(请输入控制点在地面摄影测量坐标系的坐标的Y{0}值：\r,k+1);stringv=Console.ReadLine();arr2[k,1]=(double)Convert.ToSingle(v);}//3.控制点Z坐标的输入for(intp=0;pn;p++){Console.WriteLine(请输入控制点在地面摄影测量坐标系的坐标的Z{0}值：\r,p+1);stringw=Console.ReadLine();arr2[p,2]=(double)Convert.ToSingle(w);}//确定外方位元素的初始值//1.确定Xs的初始值：doubleXs0=0;doublesumx=0;for(intj=0;jn;j++){doubleh=arr2[j,0];sumx+=h;}Xs0=sumx/n;//2.确定Ys的初始值：doubleYs0=0;doublesumy=0;for(intj=0;jn;j++){doubleh=arr2[j,1];sumy+=h;}Ys0=sumy/n;//3.确定Zs的初始值：doubleZs0=0;doublesumz=0;for(intj=0;j=n-1;j++){doubleh=arr2[j,2];sumz+=h;}Zs0=sumz/n;doubleΦ0=0;doubleΨ0=0;doubleK0=0;Console.WriteLine(Xs0,Ys0,Zs0,Φ0,Ψ0,K0的值分别是：{0},{1},{2},{3},{4},{5},Xs0,Ys0,Zs0,0,0,0);//用三个角元素的初始值按（3-4-5）计算各方向余弦值，组成旋转矩阵,此时的旋转矩阵为单位矩阵I：double[,]arr3=newdouble[3,3];for(inti=0;i3;i++){arr3[i,i]=1;}doublea1=arr3[0,0];doublea2=arr3[0,1];doublea3=arr3[0,2];doubleb1=arr3[1,0];doubleb2=arr3[1,1];doubleb3=arr3[1,2];doublec1=arr3[2,0];doublec2=arr3[2,1];doublec3=arr3[2,2];/*利用线元素的初始值和控制点的地面坐标，代入共线方程（3-5-2），*逐点计算像点坐标的近似值*///1.定义存放像点近似值的数组double[]arr4=newdouble[2*n];//----------近似值矩阵//2.逐点像点坐标计算近似值//a.计算像点的x坐标近似值（x）for(inti=0;i2*n;i+=2){for(intj=0;jn;j++){arr4[i]=-f*(a1*(arr2[j,0]-Xs0)+b1*(arr2[j,1]-Ys0)+c1*(arr2[j,2]-Zs0))/(a3*(arr2[j,0]-Xs0)+b3*(arr2[j,1]-Ys0)+c3*(arr2[j,2]-Zs0));}}//b.计算像点的y坐标近似值（y）for(inti=1;i2*n;i+=2){for(intj=0;jn;j++){arr4[i]=-f*(a2*(arr2[j,0]-Xs0)+b2*(arr2[j,1]-Ys0)+c2*(arr2[j,2]-Zs0))/(a3*(arr2[j,0]-Xs0)+b3*(arr2[j,1]-Ys0)+c3*(arr2[j,2]-Zs0));}}//逐点计算误差方程式的系数和常数项，组成误差方程：double[,]arr5=newdouble[2*n,6];//------------系数矩阵（A）//1.计算dXs的系数for(inti=0;i2*n;i+=2){arr5[i,0]=-1/m;//-f/H==-1/m}//2.计算dYs的系数for(inti=1;i2*n;i+=2){arr5[i,1]=-1/m;//-f/H==-1/m}//3.a.计算误差方程式Vx中dZs的系数for(inti=0;i2*n;i+=2){arr5[i,2]=-arr1[i]/m*f;}//3.b.计算误差方程式Vy中dZs的系数for(inti=1;i2*n;i+=2){arr5[i,2]=-arr1[i]/m*f;}//4.a.计算误差方程式Vx中dΦ的系数for(inti=0;i2*n;i+=2){arr5[i,3]=-f*(1+arr1[i]*arr1[i]/f*f);}//4.a.计算误差方程式Vy中dΦ的系数for(inti=1;i2*n;i+=2){arr5[i,3]=-arr1[i-1]*arr1[i]/f;}//5.a.计算误差方程式Vx中dΨ的系数for(inti=0;i2*n;i+=2){arr5[i,4]=-arr1[i]*arr1[i+1]/f;}//5.b.计算误差方程式Vy中dΨ的系数for(inti=1;i2*n;i+=2){arr5[i,4]=-f*(1+arr1[i]*arr1[i]/f*f);}//6.a.计算误差方程式Vx中dk的系数for(inti=0;i2*n;i+=2){arr5[i,5]=arr1[i+1];}//6.b.计算误差方程式Vy中dk的系数for(inti=1;i2*n;i+=2){arr5[i,5]=-arr1[i-1];}//定义外方位元素组成的数组double[]arr6=newdouble[6];//--------------------外方位元素改正数矩阵（X）//定义常数项元素组成的数组double[]arr7=newdouble[2*n];//-----------------常数矩阵（L）//计算lx的值for(inti=0;i2*n;i+=2){arr7[i]=arr1[i]-arr4[i];//将近似值矩阵的元素代入}//计算ly的值for(inti=1;i=2*(n-1);i+=2){arr7[i]=arr1[i]-arr4[i];//将近似值矩阵的元素代入}/*对于所有像点的坐标观测值，一般认为是等精度量测，所以权阵P为单位阵.所以X=(ATA)-1ATL*///1.计算ATdouble[,]arr5T=newdouble[6,2*n];for(inti=0;i6;i++){for(intj=0;j2*n;j++){arr5T[i,j]=arr5[j,i];}}//A的转置与A的乘积,存放在arr5AA中double[,]arr5AA=newdouble[6,6];for(inti=0;i6;i++){for(intj=0;j6;j++){arr5AA[i,j]=0;for(intl=0;l2*n;l++){arr5AA[i,j]+=arr5T[i,l]*arr5[l,j];}}}nijuzhen(arr5AA);//arr5AA经过求逆后变成原矩阵的逆矩阵//arr5AA*arr5T存在arr5AARATdouble[,]arr5AARAT=newdouble[6,2*n];for(inti=0;i6;i++){for(intj=0;j2*n;j++){arr5AARAT[i,j]=0;for(intp=0;p6;p++){arr5AARAT[i,j]+=arr5AA[i,p]*arr5T[p,j];}}}//计算arr5AARATxL，存在arrX中double[]arrX=newdouble[6];for(inti=0;i6;i++){for(intj=0;j1;j++){arrX[i]=0;for(intvv=0;vv6;vv++){arrX[i]+=arr5AARAT[i,vv]*arr7[vv];}}}//计算外方位元素值doubleXs,Ys,Zs,Φ,Ψ,K;Xs=Xs0+arrX[0];Ys=Ys0+arrX[1];Zs=Zs0+arrX[2];Φ=Φ0+arrX[3];Ψ=Ψ0+arrX[4];K=K0+arrX[5];for(inti=0;i=2;i++){Xs+=arrX[0];Ys+=arrX[1];Zs+=arrX[2];Φ+=arrX[3];Ψ+=arrX[4];K+=arrX[5];}Console.WriteLine(Xs,Ys,Zs,Φ,Ψ,K的值分别是：