《摄影测量学》综合实习一：解析空中三角测量程序设计
附上实习报告供参考
第一个像对的坐标
648790 735260 -230980 788550
113860 595800 -771380 679190
241050 586260 -644640 661950
578030 223960 -327700 277590
256140 187820 -648360 260160
606230 -104550 -315660 -51440
278340 -565020 -660020 -487200
179220 -757800 -765430 -670400
478510 -637940 -466930 -569840
975930 -700850 16690 -658940
917380 -16160 -4600 18540
803150 758570 -76440 802960
94870 709670 -785850 795830
35030 -34550 -877140 50930
19790 -843710 -925660 -745370
第二个像对的坐标
426790 -634450 -568240 -631500
323920 637060 -601170 642970
16690 -658940 -980300 -630600
-4600 18540 -963100 51030
-76410 802880 -989450 836700
921460 679700 -8790 641530
981740 -884160 -29060 -915900
965780 38900 1460 5490
第三个像对的坐标
931930 729560 -26020 680090
907100 -785750 20330 -836630
753660 -838360 -134160 -896670
397770 -725220 -500070 -791380
49010 50160 -875490 -17320
918060 48330 -10150 11950
890340 544420 -58540 499270
624000 444520 -317940 392160
-8690 641540 -945780 560530
-29020 -915940 -930070 -1004070
1460 5460 -921980 -63410
控制点的坐标
24204.689 46604.652 46.251
24343.683 45111.263 48.198
24965.047 44309.253 49.340
22167.904 46432.515 46.457
f＝153.033mm
bx＝200mm
如何插入一段漂亮的代码片
用matlab进行程序设计的代码
%进行相对定向代码
function [u,v,fai,omig,kb,X1,Y1,Z1,N1,N2,bz,bx,by,X2,Y2,Z2]=XDDX(x1,y1,x2,y2,fai0,omig0,kb0)
clc
bx = 200/1000;
f = 153.033/1000;
%g构造最小二乘方程
A = zeros(length(x1),5);
XXX = zeros(5,1);
u=0*pi/180;
v=0*pi/180;
fai=fai0;
omig =omig0;
kb=kb0;
a1=cos(fai)*cos(kb)-sin(fai)*sin(omig)*sin(kb);
a2=-cos(fai)*sin(kb)-sin(fai)*sin(omig)*cos(kb);
a3=-sin(fai)*cos(omig);
b1=cos(omig)*sin(kb);
b2=cos(omig)*cos(kb);
b3=-sin(omig);
c1=sin(fai)*cos(kb)+cos(fai)*sin(omig)*sin(kb);
c2=-sin(fai)*sin(kb)+cos(fai)*sin(omig)*cos(kb);
c3=cos(fai)*cos(omig);
R1=[a1 a2 a3;b1 b2 b3;c1 c2 c3];
X1 = zeros(length(x1),1);
Y1 = zeros(length(x1),1);
Z1 = zeros(length(x1),1);
for i=1:length(x1)
temp = zeros(3,1);
temp= R1*[x1(i);y1(i);-f];
X1(i) = temp(1);
Y1(i) = temp(2);
Z1(i)= temp(3);
end
by = bx*u;
bz = bx*v;
R2 = zeros(3,3);
X2 = zeros(length(x1),1);
Y2 = zeros(length(x1),1);
Z2 = zeros(length(x1),1);
R2 = [1 0 0;0 1 0;0 0 1];
for i=1:length(x2)
temp = zeros(3,1);
temp= R2*[x2(i);y2(i);-f];
X2(i) = temp(1);
Y2(i) = temp(2);
Z2(i)= temp(3);
end
N1 =zeros(length(x1),1);
N2 = zeros(length(x2),1);
Q =zeros(length(x1),1);
for i=1:length(x1)
N1(i) = (bx*Z2(i)-bz*X2(i))/(X1(i)*Z2(i)-X2(i)*Z1(i));
N2(i) = (bx*Z1(i)-bz*X1(i))/(X1(i)*Z2(i)-X2(i)*Z1(i));
Q(i) = N1(i)*Y1(i)-N2(i)*Y2(i)-by;
end
%组成误差方程
for i=1:length(x1)
A(i,1) = bx;
A(i,2) = -(Y2(i)/Z2(i))*bx;
A(i,3) = -(X2(i)*Y2(i)/Z2(i))*N2(i);
A(i,4) = -(Z2(i)+Y2(i)*Y2(i)/Z2(i))*N2(i);
A(i,5) = X2(i)*N2(i);
end
format longg
XXX = (inv(A'*A))*(A'*Q);
u=u+XXX(1);
v=v+XXX(2);
fai=fai+XXX(3);
omig =omig+XXX(4);
kb=kb+XXX(5);
n=1;
while(abs(XXX(1))>=0.00003||abs(XXX(2))>=0.00003||abs(XXX(3))>=0.00003||abs(XXX(4))>=0.00003||abs(XXX(5))>=0.00003)
n=n+1;
by = u*bx;
bz = v*bx;
a1=cos(fai)*cos(kb)-sin(fai)*sin(omig)*sin(kb);
a2=-cos(fai)*sin(kb)-sin(fai)*sin(omig)*cos(kb);
a3=-sin(fai)*cos(omig);
b1=cos(omig)*sin(kb);
b2=cos(omig)*cos(kb);
b3=-sin(omig);
c1=sin(fai)*cos(kb)+cos(fai)*sin(omig)*sin(kb);
c2=-sin(fai)*sin(kb)+cos(fai)*sin(omig)*cos(kb);
c3=cos(fai)*cos(omig);
R2=[a1 a2 a3;b1 b2 b3;c1 c2 c3];
for i=1:length(x1)
temp = zeros(3,1);
temp= R2*[x2(i);y2(i);-f];
X2(i) = temp(1);
Y2(i) = temp(2);
Z2(i)= temp(3);
end
for i=1:length(x1)
N1(i) = (bx*Z2(i)-bz*X2(i))/(X1(i)*Z2(i)-X2(i)*Z1(i));
N2(i) = (bx*Z1(i)-bz*X1(i))/(X1(i)*Z2(i)-X2(i)*Z1(i));
Q(i) = N1(i)*Y1(i)-N2(i)*Y2(i)-by;
end
for i=1:length(x1)
A(i,1) = bx;
A(i,2) = -(Y2(i)/Z2(i))*bx;
A(i,3) = -(X2(i)*Y2(i)/Z2(i))*N2(i);
A(i,4) = -(Z2(i)+Y2(i)*Y2(i)/Z2(i))*N2(i);
A(i,5) = X2(i)*N2(i);
end
XXX = (inv(A'*A))*(A'*Q);
u=u+XXX(1);
v=v+XXX(2);
fai=fai+XXX(3);
omig =omig+XXX(4);
kb=kb+XXX(5);
end
%可以单独设置一个脚本文件 进行归化系数的求解
function k =GHXS(N11,Z11,bz1,N21,Z21)
k = ((N11*Z11)-bz1)/(N21*Z21);
%%%%主函数
[a,b,c,d]=textread('像对1.txt','%n%n%n%n');
tempzz =textread('像对1.txt');
[u1,v1,fai1,omig1,kb1,X11,Y11,Z11,N11,N12,bz1,bx1,by1,X12,Y12,Z12]=XDDX(a/1000000,b/1000000,c/1000000,d/1000000,0,0,0)
[a2,b2,c2,d2]=textread('像对2.txt','%n%n%n%n');
[u2,v2,fai2,omig2,kb2,X21,Y21,Z21,N21,N22,bz2,bx2,by2,X22,Y22,Z22]=XDDX(a2/1000000,b2/1000000,c2/1000000,d2/1000000,fai1,omig1,kb1);
[a3,b3,c3,d3]=textread('像对3.txt','%n%n%n%n');
[u3,v3,fai3,omig3,kb3,X31,Y31,Z31,N31,N32,bz3,bx3,by3,X32,Y32,Z32]=XDDX(a3/1000000,b3/1000000,c3/1000000,d3/1000000,fai2,omig2,kb2);
%求解比例尺规划系数
[temp1] = GHXS(N11(10),Z11(10),bz1,N21(3),Z22(3));
[temp2] = GHXS(N11(11),Z11(11),bz1,N21(4),Z22(4));
[temp3] = GHXS(N11(12),Z11(12),bz1,N21(5),Z22(5));
k2=(temp1+temp2+temp3)/3;
[temp1] = GHXS(N21(3),Z22(3),bz2,N31(9),Z32(9));
[temp2] = GHXS(N21(4),Z22(4),bz2,N31(10),Z32(10));
[temp3] = GHXS(N21(5),Z22(5),bz2,N31(11),Z32(11));
k3=(temp1+temp2+temp3)/3;
%各个基分量乘以归化系数
bx2 = bx2*k2;
by2 =by2*k2;
bz2 = bz2*k2;
bx3 = bx3*k3;
by3 =by3*k3;
bz3 = bz3*k3;
%%%%%%%%%%%%%%%%%%
%计算比例尺m
tempxx=[tempzz(1,:);tempzz(6,:);tempzz(7,:)];
tempx = tempxx(:,1)/1000000;
tempy = tempxx(:,2)/1000000;
s=0;
S=0;
for i=1:2
j=i+1;
sij=sqrt((tempx(i)-tempx(j))^2+(tempy(i)-tempy(j))^2);
Sij=sqrt((kzx(i)-kzx(j))^2+(kzy(i)-kzy(j))^2);
s=s+sij;
S=S+Sij;
end
m =S/s
%%%%%%%%%%%%%%%%%%%%%
%模型点摄影测量坐标系计算统一到第一个模型的左像片中
%第一个模型左像片摄影中心摄影测量坐标
XPS1 = 0;
YPS1 = 0;
ZPS1 =m*153.033/1000;
%%%%
%第一个模型右像片摄影中心摄影测量坐标
XPS2 =XPS1+ m*bx1;
YPS2 =YPS1+ m*by1;
ZPS2 =m*153.033/1000+bz1;
%%%%
%第一个模型像点摄影测量坐标
Xp1 = m*N11.*X11;
Yp1= m* 0.5*(N11.*Y11+N12.*Y12+by1);
Zp1=m*N11.*Z11+m*153.033/1000;
%%%
%第二个模型右相片摄影中心摄影测量坐标
XPS3 = XPS2 +k2*m*bx2;
YPS3 = YPS2 +k2*m*by2;
ZPS3 = ZPS2 +k2*m*bz2;
%%%
%第二个模型像点的摄影测量坐标
Xp2 = XPS2+k2*m*N21.*X21;
Yp2 = 0.5*(YPS2+k2*m*N21.*Y21+YPS3+k2*m*N22.*Y22)
Zp2 = ZPS2 +k2*m*N21.*Z21;
%%%
%第三个模型右像片摄影中心摄影测量坐标
XPS4 = XPS3+k3*m*bx3;
YPS4 = YPS3+k3*m*by3;
ZPS4 = ZPS3+k3*m*bz3;
%%%
%第三个模型像点摄影测量坐标
Xp3 = XPS3+k3*m*N31.*X31;
Yp3 = 0.5*(YPS3+k3*m*N31.*Y31+YPS4+k3*m*N32.*Y32)
Zp3 = ZPS3 +k2*m*N31.*Z31;
%%%%%%%%%%%
%%%%%%%%%%%%%%%%摄影测量坐标系转化完毕
%利用重心化进行绝对定向
%根据已知控制点确定大地坐标系与地面摄影测量坐标系之间的转换系数a,b
change_Xt = kzx(1)-kzx(2);
change_Yt =kzy(1)-kzy(2);
change_Xp = Xp1(1)-Xp1(6);
change_Yp = Yp1(1)-Yp1(6);
a = (change_Xp*change_Yt+change_Yp*change_Xt)/(change_Xt*change_Xt+change_Yt*change_Yt);
b = (change_Xp*change_Xt-change_Yp*change_Yt)/(change_Xt*change_Xt+change_Yt*change_Yt);
lameda = sqrt(a^2+b^2);
%%%%%%%%%%%%%%%%%%%
%地面摄影测量坐标系重心化
Xt=[24204.689,24343.683,24965.047,22167.904];
Yt=[46604.652,45111.263,44309.253,46432.515];
Zt=[46.251,48.198,49.340,46.457];
%大地坐标系转化为地面摄影测量坐标系
for i=1:4
Xtp(i)=b*(Xt(i)-Xt(1))+a*(Yt(i)-Yt(1));
Ytp(i)=a*(Xt(i)-Xt(1))-b*(Yt(i)-Yt(1));
Ztp(i)=lamda*(Zt(i)-Zt(1));
end
%%%%%%
%重心化地面摄影测量坐标系
Xtpg = mean(Xtp);
Ytpg = mean(Ytp);
Ztpg = mean(Ztp);
Xtpba =Xtp - Xtpg;
Ytpba =Ytp - Ytpg;
Ztpba =Ztp - Ztpg;
%%%%
%摄影测量坐标系重心化
Xpg =[Xp1(1) Xp1(6) Xp1(7) Xp2(6)];
Ypg = [Yp1(1) Yp1(6) Yp1(7) Yp2(6)];
Zpg = [Zp1(1) Zp1(6) Zp1(7) Zp2(6)];
Xpg =mean(Xpg);
Ypg =mean(Ypg);
Zpg =mean(Zpg);
Xpba(1)=Xp1(1)-Xpg;
Xpba(2)=Xp1(6)-Xpg;
Xpba(3)=Xp1(7)-Xpg;
Xpba(4)=Xp2(6)-Xpg;
Ypba(1)=Yp1(1)-Ypg;
Ypba(2)=Yp1(6)-Ypg;
Ypba(3)=Yp1(7)-Ypg;
Ypba(4)=Yp2(6)-Ypg;
Zpba(1)=Zp1(1)-Zpg;
Zpba(2)=Zp1(6)-Zpg;
Zpba(3)=Zp1(7)-Zpg;
Zpba(4)=Zp2(6)-Zpg;
%%%%%%%
%用法方程求解绝对定向元素
A=[1,0,0,Xpba(1),-Zpba(1),0,-Ypba(1);
0,1,0,Ypba(1),0,-Zpba(1),Xpba(1);
0,0,1,Zpba(1),Xpba(1),Ypba(1),0;
1,0,0,Xpba(2),-Zpba(2),0,-Ypba(2);
0,1,0,Ypba(2),0,-Zpba(2),Xpba(2);
0,0,1,Zpba(2),Xpba(2),Ypba(2),0;
1,0,0,Xpba(3),-Zpba(3),0,-Ypba(3);
0,1,0,Ypba(3),0,-Zpba(3),Xpba(3);
0,0,1,Zpba(3),Xpba(3),Ypba(3),0;
1,0,0,Xpba(4),-Zpba(4),0,-Ypba(4);
0,1,0,Ypba(4),0,-Zpba(4),Xpba(4);
0,0,1,Zpba(4),Xpba(4),Ypba(4),0;
];
for i=1:4
lx(i)=Xtpba(i)-Xpba(i);
ly(i)=Ytpba(i)-Ypba(i);
lz(i)=Ztpba(i)-Zpba(i);
end
L=[lx(1);ly(1);lz(1);lx(2);ly(2);lz(2);lx(3);ly(3);lz(3);lx(4);ly(4);lz(4)];
deta7=inv(A'*A)*A'*L;
%%%绝对定向结束，求得七个绝对定向元素
deta7x=deta7(1);
deta7y=deta7(2);
deta7z=deta7(3);
deta7lameda=deta7(4)+1;
deta7p=deta7(5);
deta7w=deta7(6);
deta7k=deta7(7);
%%%%%%%
%%%%%下一步即把像点由摄影测量坐标系转化为大地坐标系过度为地面摄影测量坐标系
a1=cos(deta7p)*cos(deta7k)-sin(deta7p)*sin(deta7w)*sin(deta7k);
a2=-cos(deta7p)*sin(deta7k)-sin(deta7p)*sin(deta7w)*cos(deta7k);
a3=-sin(deta7p)*cos(deta7w);
b1=cos(deta7w)*sin(deta7k);
b2=cos(deta7w)*cos(deta7k);
b3=-sin(deta7w);
c1= sin(deta7p)*cos(deta7k)+cos(deta7p)*sin(deta7w)*sin(deta7k);
c2=-sin(deta7p)*sin(deta7k)+cos(deta7p)*sin(deta7w)*cos(deta7k);
c3=cos(deta7p)*cos(deta7w);
R=[a1,a2,a3;b1,b2,b3;c1,c2,c3];
%经过概略定向后得控制点的地面摄影测量坐标系
XbaYbaZba1=deta7lameda*R*[Xpba(1);Ypba(1);Zpba(1)]+[deta7x;deta7y;deta7z];
XbaYbaZba2=deta7lameda*R*[Xpba(2);Ypba(2);Zpba(2)]+[deta7x;deta7y;deta7z];
XbaYbaZba3=deta7lameda*R*[Xpba(3);Ypba(3);Zpba(3)]+[deta7x;deta7y;deta7z];
XbaYbaZba4=deta7lameda*R*[Xpba(4);Ypba(4);Zpba(4)]+[deta7x;deta7y;deta7z];
%%%%进行非线性化改正
B=[1,XbaYbaZba1(1),XbaYbaZba1(2),XbaYbaZba1(1)^2,XbaYbaZba1(1)*XbaYbaZba1(2);
1,XbaYbaZba2(1),XbaYbaZba2(2),XbaYbaZba2(1)^2,XbaYbaZba2(1)*XbaYbaZba2(2);
1,XbaYbaZba3(1),XbaYbaZba3(2),XbaYbaZba3(1)^2,XbaYbaZba3(1)*XbaYbaZba3(2);
1,XbaYbaZba4(1),XbaYbaZba4(2),XbaYbaZba4(1)^2,XbaYbaZba4(1)*XbaYbaZba4(2)];
lxx=[Xtpba(1)-XbaYbaZba1(1);Xtpba(2)-XbaYbaZba2(1);Xtpba(3)-XbaYbaZba3(1);Xtpba(4)-XbaYbaZba4(1)];
lyy=[Ytpba(1)-XbaYbaZba1(2);Ytpba(2)-XbaYbaZba2(2);Ytpba(3)-XbaYbaZba3(2);Ytpba(4)-XbaYbaZba4(2)];
lzz=[Ztpba(1)-XbaYbaZba1(3);Ztpba(2)-XbaYbaZba2(3);Ztpba(3)-XbaYbaZba3(3);Ztpba(4)-XbaYbaZba4(3)];
AAA=inv(B'*B)*B'*lxx;
BBB=inv(B'*B)*B'*lyy;
CCC=inv(B'*B)*B'*lzz;
a00=AAA(1);a11=AAA(2);a22=AAA(3);a33=AAA(4);a44=AAA(5);
b00=BBB(1);b11=BBB(2);b22=BBB(3);b33=BBB(4);b44=BBB(5);
c00=CCC(1);c11=CCC(2);c22=CCC(3);c33=CCC(4);c44=CCC(5);
Xba=[XbaYbaZba1(1),XbaYbaZba2(1),XbaYbaZba3(1),XbaYbaZba4(1)];
Yba=[XbaYbaZba1(2),XbaYbaZba2(2),XbaYbaZba3(2),XbaYbaZba4(2)];
Zba=[XbaYbaZba1(3),XbaYbaZba2(3),XbaYbaZba3(3),XbaYbaZba4(3)];
for i=1:4
Xtpnew(i)=(Xtpg*10000+Xba(i)*10000+a00+a11*Xba(i)+a22*Yba(i)+a33*Xba(i)*Xba(i)+a44*Xba(i)*Yba(i))/10000;
Ytpnew(i)=(Ytpg*10000+Yba(i)*10000+b00+b11*Xba(i)+b22*Yba(i)+b33*Xba(i)*Xba(i)+b44*Xba(i)*Yba(i))/10000;
Ztpnew(i)=(Ztpg*10000+Zba(i)*10000+c00+c11*Xba(i)+c22*Yba(i)+c33*Xba(i)*Xba(i)+c44*Xba(i)*Yba(i))/10000;
end
%求大地坐标，由摄影测量坐标系到地面摄影测量坐标系到大地坐标系
dian1=[];%相当于Xba
dian2=[];
dian3=[];
%转化为以重心点为原点的摄影测量坐标系
for i=1:15
Xpba1(i)=Xp1(i)-Xpg;
Ypba1(i)=Yp1(i)-Ypg;
Zpba1(i)=Zp1(i)-Zpg;
dian=deta7lameda*R*[Xpba1(i);Ypba1(i);Zpba1(i)]+[deta7x;deta7y;deta7z];
dian1=[dian1,dian];
end
for i=1:8
Xpba2(i)=Xp2(i)-Xpg;
Ypba2(i)=Yp2(i)-Ypg;
Zpba2(i)=Zp2(i)-Zpg;
dian=deta7lameda*R*[Xpba2(i);Ypba2(i);Zpba2(i)]+[deta7x;deta7y;deta7z];
dian2=[dian2,dian];
end
for i=1:11
Xpba3(i)=Xp3(i)-Xpg;
Ypba3(i)=Yp3(i)-Ypg;
Zpba3(i)=Zp3(i)-Zpg;
dian=deta7lameda*R*[Xpba3(i);Ypba3(i);Zpba3(i)]+[deta7x;deta7y;deta7z];
dian3=[dian3,dian];
end
dadizuobiao1=[];
dadizuobiao2=[];
dadizuobiao3=[];
%%%进行两步，第一步由摄影测量坐标系到地面摄影测量坐标系 第二步由地面摄影测量坐标系到大地坐标系
for i=1:15
Xtpdian1(i)=(Xtpg*10000+dian1(1,i)*10000+a00+a11*dian1(1,i)+a22*dian1(2,i)+a33*dian1(1,i)*dian1(1,i)+a44*dian1(1,i)*dian1(2,i))/10000;
Ytpdian1(i)=(Ytpg*10000+dian1(2,i)*10000+b00+b11*dian1(1,i)+b22*dian1(2,i)+b33*dian1(1,i)*dian1(1,i)+b44*dian1(1,i)*dian1(2,i))/10000;
Ztpdian1(i)=(Ztpg*10000+dian1(3,i)*10000+c00+c11*dian1(1,i)+c22*dian1(2,i)+c33*dian1(1,i)*dian1(1,i)+c44*dian1(1,i)*dian1(2,i))/10000;
dadizuobiao=1/lameda/lameda*[b,a,0;a,-b,0;0,0,lameda]*[Xtpdian1(i);Ytpdian1(i);Ztpdian1(i)]+[Xt(1);Yt(1);Zt(1)]
dadizuobiao1=[dadizuobiao1,dadizuobiao];
end
for i=1:8
Xtpdian2(i)=(Xtpg*10000+dian2(1,i)*10000+a00+a11*dian2(1,i)+a22*dian2(2,i)+a33*dian2(1,i)*dian2(1,i)+a44*dian2(1,i)*dian2(2,i))/10000;
Ytpdian2(i)=(Ytpg*10000+dian2(2,i)*10000+b00+b11*dian2(1,i)+b22*dian2(2,i)+b33*dian2(1,i)*dian2(1,i)+b44*dian2(1,i)*dian2(2,i))/10000;
Ztpdian2(i)=(Ztpg*10000+dian2(3,i)*10000+c00+c11*dian2(1,i)+c22*dian2(2,i)+c33*dian2(1,i)*dian2(1,i)+c44*dian2(1,i)*dian2(2,i))/10000;
dadizuobiao=1/lameda/lameda*[b,a,0;a,-b,0;0,0,lameda]*[Xtpdian2(i);Ytpdian2(i);Ztpdian2(i)]+[Xt(1);Yt(1);Zt(1)]
dadizuobiao2=[dadizuobiao2,dadizuobiao];
end
for i=1:11
Xtpdian3(i)=(Xtpg*10000+dian3(1,i)*10000+a00+a11*dian3(1,i)+a22*dian3(2,i)+a33*dian3(1,i)*dian3(1,i)+a44*dian3(1,i)*dian3(2,i))/10000;
Ytpdian3(i)=(Ytpg*10000+dian3(2,i)*10000+b00+b11*dian3(1,i)+b22*dian3(2,i)+b33*dian3(1,i)*dian3(1,i)+b44*dian3(1,i)*dian3(2,i))/10000;
Ztpdian3(i)=(Ztpg*10000+dian3(3,i)*10000+c00+c11*dian3(1,i)+c22*dian3(2,i)+c33*dian3(1,i)*dian3(1,i)+c44*dian3(1,i)*dian3(2,i))/10000;
dadizuobiao=1/lameda/lameda*[b,a,0;a,-b,0;0,0,lameda]*[Xtpdian3(i);Ytpdian3(i);Ztpdian3(i)]+[Xt(1);Yt(1);Zt(1)]
dadizuobiao3=[dadizuobiao3,dadizuobiao];
end
9104
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
