摄影测量-空间前交、后交

空间后交-前交程序设计

(实验报告)

姓名： 班级： 学号： 时间：

空间后交-前交程序设计

一、实验目的

用 C 、VB或MATLAB语言编写空间后方交会-空间前方交会程序 ⑴提交实习报告：程序框图、程序源代码、计算结果、体会 ⑵计算结果：像点坐标、地面坐标、单位权中误差、外方位元素及其精度

二、实验数据

f=150.000mm，x0=0，y0=0

三、实验思路

1.利用空间后方交会求左右像片的外方位元素

(1).获取m(于像片中选取两点，于地面摄影测量坐标系中选取同点，分别计算距离，距离比值即为m)，x,y,f,X,Y,Z (2).确定未知数初始值Xs,Ys,Zs,q,w,k (3).计算旋转矩阵R

(4).逐点计算像点坐标的近似值（x），（y）

(5).组成误差方程式 (6).组成法方程式 (7).解求外方位元素

(8).检查是否收敛，即将求得的外方位元素的改正数与规定限差比较，小于限差即终止；否则用新的近似值重复步骤（3）-（7） 2.利用求出的外方位元素进行空间前交，求出待定点地面坐标 (1).用各自像片的角元素计算出左、右像片的方向余弦值，组成旋转矩阵R1,R2

(2).根据左、右像片的外方位元素 ，计算摄影基线分量Bx，By，Bz (3).计算像点的像空间辅助坐标（X1,Y1,Z1）和（X2,Y2,Z2） (4).计算点投影系数N1和N2 (5).计算未知点的地面摄影测量坐标

四、实验过程

⑴程序框图

函数AandL

%求间接平差时需要的系数

%%%已知

%a=像点坐标x,b=像点坐标y，f内方位元素主距 %φ=q,ψ=w,κ=k %像空间坐标系X,Y,Z

%地面摄影测量坐标系Xs,Ys,Zs

function [A1,L1,A2,L2]=AandL(a,b,f,q,w,k,X,Y,Z,Xs,Ys,Zs) %%%%%%%%%%%选择矩阵元素

a1=cos(q)*cos(k)-sin(q)*sin(w)*sin(k); a2=-cos(q)*sin(k)-sin(q)*sin(w)*cos(k); a3=-sin(q)*cos(w); b1=cos(w)*sin(k); b2=cos(w)*cos(k); b3=-sin(w);

c1=sin(q)*cos(k)+cos(q)*sin(w)*sin(k); c2=-sin(q)*sin(k)+cos(q)*sin(w)*cos(k); c3=cos(q)*cos(w);

%%%%%%%共线方程的分子分母

X_=a1*(X-Xs)+b1*(Y-Ys)+c1*(Z-Zs); Y_=a2*(X-Xs)+b2*(Y-Ys)+c2*(Z-Zs); Z_=a3*(X-Xs)+b3*(Y-Ys)+c3*(Z-Zs); %%%%%%%近似值 x=-f*X_/Z_; y=-f*Y_/Z_;

%%%%%%%A组成L组成 a11=1/Z_*(a1*f+a3*x); a12=1/Z_*(b1*f+b3*x); a13=1/Z_*(c1*f+c3*x); a21=1/Z_*(a2*f+a3*y); a22=1/Z_*(b2*f+b3*y); a23=1/Z_*(c2*f+c3*y);

a14=y*sin(w)-(x/f*(x*cos(k)-y*sin(k))+f*cos(k))*cos(w);

a15=-f*sin(k)-x/f*(x*sin(k)+y*cos(k)); a16=y;

a24=-x*sin(w)-(y/f*(x*cos(k)-y*sin(k))-f*sin(k))*cos(w); a25=-f*cos(k)-y/f*(x*sin(k)+y*cos(k)); a26=-x; lx=a-x; ly=b-y;

%%%%%%%%%组成一个矩阵，并返回 A1=[a11,a12,a13,a14,a15,a16]; A2=[a21,a22,a23,a24,a25,a26]; L1=lx; L2=ly;

函数deg2dms

%%%%%%%%角度转度分秒 function y=deg2dms(x) a=floor(x);

b=floor((x-a)*60); c=(x-a-b/60)*3600; y=a+(b/100)+(c/10000);

函数dms2deg

%%%%%度分秒转度

function y=dms2deg(x) a=floor(x);

b=floor((x-a)*100); c=(x-a-b/100)*10000; y=a+b/60+c/3600;

函数ok

%%%%%%%%%%%

%%%目的是为了保证各取的值的有效值 %%xy为n*1,a为1*n

function result=ok(xy,a) format short g i=size(xy,1); for n=1:i

o=xy(n)-floor(xy(n,1));

o=round(o*(10^a(n)))/(10^a(n)); xy(n,1)=floor(xy(n,1))+o; end

format long g result=xy;

函数rad2dmsxy

%%%%求度分秒表现形式的三个外方位元素，三个角度 function xydms=rad2dmsxy(xy) [a,b,c,d,e,f]=testvar(xy); d=deg2dms(rad2deg(d)); e=deg2dms(rad2deg(e)); f=deg2dms(rad2deg(f)); xydms=[a,b,c,d,e,f]';

函数spacehoujiao

%%%%%%%空间后交 %%% f

%%输入p（2*n,1）

%%像点坐标x,y,X,Y,Z,均为（n,1）

function [xy,m,R]=spacehoujiao(p,x,y,f,X,Y,Z) format long;

%%%%%权的矢量化,这是等精度时的，如果非，将函数参数改为P P=diag(p);

%%求n

j=size(X,2); %%初始化

Xs=0;Ys=0;Zs=0; for n=1:j Xs=Xs+X(n); Ys=Ys+Y(n); Zs=Zs+Z(n); end

Sx=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2);%%%%两像点之间距离

Sd=sqrt((X(2)-X(1))^2+(Y(2)-Y(1))^2);%%%%两地面控制点之间距离 m=Sd/Sx; %%%%图像比例系数 Xs=Xs/j; Ys=Ys/j; Zs=m*f+Zs/j; m0=0;

q=0;w=0;k=0;i=0; a=rand(2*j,6); l=rand(2*j,1); %%%%

for n=1:j

[a(2*n-1,:),l(2*n-1,1),a(2*n,:),l(2*n,1)]=AandL(x(n),y(n),f,q,w,k,X(n),Y(n),Z(n),Xs,Ys,Zs); end

det=inv(a'*P*a)*transpose(a)*P*l; %%%%%%%%%循环体 while 1

%%%%%%%%%%%%%%%%

[dXs,dYs,dZs,dq,dw,dk]=testvar(det); detXs=abs(dXs); detYs=abs(dYs); detZs=abs(dZs); detq=abs(dq); detw=abs(dw); detk=abs(dk); %%%%%%%%% if

((detXs<0.01)&&(detYs<0.01)&&(detZs<0.01)&&(detq<pi/648000)&&(detw<pi/648000)&&(detk<pi/648000)) break; else

V=(a*det-l);Q=inv(a'*P*a);

m0=m0+sqrt((V'*P*V)/(2*j-6));%%m0需要每次的改正数算出来相加 %%%

Xs=Xs+dXs; Ys=Ys+dYs; Zs=Zs+dZs; q=q+dq; w=w+dw; k=k+dk; %%%

for n=1:j

[a(2*n-1,:),l(2*n-1,1),a(2*n,:),l(2*n,1)]=AandL(x(n),y(n),f,q,w,k,X(n),Y(n),Z(n),Xs,Ys,Zs); end

det=inv(a'*P*a)*transpose(a)*P*l; i=i+1; %%%% end %%% end

[dXs,dYs,dZs,dq,dw,dk]=testvar(det); detXs=abs(dXs); detYs=abs(dYs); detZs=abs(dZs); detq=abs(dq); detw=abs(dw); detk=abs(dk); V=(a*det-l); Q=inv(a'*P*a);

m0=m0+sqrt((V'*P*V)/(2*n-6)); %%%

Xs=Xs+dXs; Ys=Ys+dYs; Zs=Zs+dZs; q=q+dq; w=w+dw; k=k+dk;

%%%%%%%%%%%%%可以输出迭代次数的i

%%%%%%%%%%%%Xs,Ys,Zs,q,w,k,i,dXs,dYs,dZs,dq,dw,dk,detXs,detYs,detZs %%%%%%%%%%%精度 mo=m0*sqrt(Q);

m=[mo(1,1),mo(2,2),mo(3,3),mo(4,4),mo(5,5),mo(6,6)]'; [mXs,mYs,mZs,mq,mw,mk]=testvar(m); %%%%%%%%%输出

xy=[Xs,Ys,Zs,q,w,k]';%%输出（6,1）的外方位元素

m=[m0,mXs,mYs,mZs,mq,mw,mk]';%%单位误差，各元素中误差 R=xyR(xy);%%旋转矩阵

函数spaceqianjiao

%%空间前交 %输入f

%输入x1,y1,x2,y2,R1,R2,xy1,xy2 （n,1） %输出X,Y,Z （n,1）

function [X,Y,Z]=spaceqianjiao(x1,y1,x2,y2,f,R1,R2,xy1,xy2) i=size(x1,2);

[Xs1,Ys1,Zs1,q1,w1,k1]=testvar(xy1); [Xs2,Ys2,Zs2,q2,w2,k2]=testvar(xy2); for n=1:i

[X1(n),Y1(n),Z1(n)]=testvar(R1*[x1(n),y1(n),-f]'); [X2(n),Y2(n),Z2(n)]=testvar(R2*[x2(n),y2(n),-f]'); Bx=Xs2-Xs1; By=Ys2-Ys1; Bz=Zs2-Zs1;

N1=(Bx*Z2(n)-Bz*X2(n))/(X1(n)*Z2(n)-X2(n)*Z1(n)); N2=(Bx*Z1(n)-Bz*X1(n))/(X1(n)*Z2(n)-X2(n)*Z1(n)); X(n)=Xs1+N1*X1(n); Z(n)=Zs1+N1*Z1(n);

Y(n)=0.5*((Ys1+N1*Y1(n))+(Ys2+N2*Y2(n))); en …… 此处隐藏：4907字，全部文档内容请下载后查看。喜欢就下载吧 ……