摄影测量实习空间后方交会

1、 摄影测量实习空间后方交会. 后方交会解算程序。（1） 原程序：D=input(请输入地面点坐标:)Y=input(请输入影像坐标点:)/1000f=input(请输入主距:)/1000C=input(请输入像主点坐标初始值:)J=input(请输入外方位角初始值:)i=size(Y,1);Y=Y;Y(1,:);D=D;D(1,:); Xs=0;Ys=0;for j=1:i m=sqrt(D(j,1)-D(j+1,1).2+(D(j,2)-D(j+1,2).2)/(Y(j,1)-Y(j+1,1).2+(Y(j,2)-Y(，2).2)/(i-1); Xs=Xs+D(j,1)/i; Ys=Ys+D

2、(j,2)/i;end Zs=m*f;for j=1:1000R(1,1)=cos(J(1)*cos(J(3)-sin(J(1)*sin(J(2)*sin(J(3);R(1,2)=-cos(J(1)*sin(J(3)-sin(J(1)*sin(J(2)*cos(J(3);R(1,3)=-sin(J(1)*cos(J(2);R(2,1)=cos(J(2)*sin(J(3);R(2,2)=cos(J(2)*cos(J(3);R(2,3)=-sin(J(2);R(3,1)=sin(J(1)*cos(J(3)+cos(J(1)*sin(J(2)*sin(J(3);R(3,2)=-sin(J(1)*si

3、n(J(3)+cos(J(1)*sin(J(2)*cos(J(3);R(3,3)=cos(J(1)*cos(J(2);for n=1:iXp(n)=R(1,1)*(D(n,1)-Xs)+R(2,1)*(D(n,2)-Ys)+R(3,1)*(D(n,3)-Zs);Yp(n)=R(1,2)*(D(n,1)-Xs)+R(2,2)*(D(n,2)-Ys)+R(3,2)*(D(n,3)-Zs);Zp(n)=R(1,3)*(D(n,1)-Xs)+R(2,3)*(D(n,2)-Ys)+R(3,3)*(D(n,3)-Zs);xc(n)=-f*Xp(n)/Zp(n);yc(n)=-f*Yp(n)/Zp(n);A

4、(2*n-1,1)=(R(1,1)*f+R(1,3)*xc(n)/Zp(n);A(2*n-1,2)=(R(2,1)*f+R(2,3)*xc(n)/Zp(n);A(2*n-1,3)=(R(3,1)*f+R(3,3)*xc(n)/Zp(n);A(2*n,1)=(R(1,2)*f+R(1,3)*yc(n)/Zp(n);A(2*n,2)=(R(2,2)*f+R(2,3)*yc(n)/Zp(n);A(2*n,3)=(R(3,2)*f+R(3,3)*yc(n)/Zp(n);A(2*n-1,4)=yc(n)*sin(J(2)-(xc(n)/f*(xc(n)*cos(J(3)-yc(n)*sin(J(3)+f

5、*cos(J(3)*cos(J(2);A(2*n-1,5)=-f*sin(J(3)-xc(n)*(xc(n)*sin(J(3)+yc(n)*cos(J(3)/f;A(2*n-1,6)=yc(n);A(2*n,4)=-xc(n)*sin(J(2)- (yc(n)/f*(xc(n)*cos(J(3)-yc(n)*sin(J(3)-f*sin(J(3)*cos(J(2);A(2*n,5)= -f*cos(J(3)-yc(n)*(xc(n)*sin(J(3)+yc(n)*cos(J(3)/f;A(2*n,6)=-xc(n);L(2*n-1,1)=Y(n,1)-xc(n);L(2*n,1)=Y(n,2)

6、-yc(n);end X=inv(A*A)*A*L; Xs=Xs+X(1); Ys=Ys+X(2); Zs=Zs+X(3); J(1)=J(1)+X(4); J(2)=J(2)+X(5); J(3)=J(3)+X(6); M=max(abs(X); if M=0.001 Xs Ys Zs R break; endend二运行实例（课本p-44 习题24）27.已知四队点的影像坐标和地点坐标如下：影像坐标地面坐标 x(mm)y(mm0X(mm)Y(mm)Z(mm)186.1568.9936589.4125273.322195.17253.4082.2137631.0831324.51728.69

7、314.7876.6339100.9724934.982386.50410.4664.4340426.5430319.81757.31 试计算近似垂直摄影情况下空间后方交会的解。假设内方位元素已知：f=153.24mm,x0=y0=0.（3）运行结果：请输入地面点坐标:36589.41 25273.32 2195.17;37631.08 31324.51 728.69;39100.97 24934.98 2386.50;40426.54 30319.91 757.31D = 1.0e+004 * 3.6589 2.5273 0.2195 3.7631 3.1325 0.0729 3.9101

8、2.4935 0.2387 4.0427 3.0320 0.0757请输入影像坐标点:-86.15 -68.99;-53.40 82.21;-14.78 -76.63;10.46 64.43Y = -0.0862 -0.0690 -0.0534 0.0822 -0.0148 -0.0766 0.0105 0.0644请输入主距:153.24f = 0.1532请输入像主点坐标初始值:0 0 C = 0 0 0请输入外方位角初始值:0 0 0J = 0 0 0Xs = 3.9795e+004Ys = 2.7477e+004Zs = 7.5728e+003R = 0.9977 0.0675 0.0

9、039 -0.0675 0.9977 -0.0021 -0.0041 0.0018 1.0000二Moravec点特征提取算法程序。（1）原程序G=input(请输入搜索区灰度矩阵：)w=input(请输入窗口大小（象元表示;）：)s=size(G);if max(s)y)（2）运行实例（自由编辑一组灰度值）45 56 87 98 98 26 34 54 55;23 54 67 87 42 32 33 65 44;43 87 98 44 22 44 33 55 54 ;26 34 76 43 75 83 84 98 23;25 65 87 87 33 32 43 54 63（3）运行结果：请输

10、入搜索区灰度矩阵：45 56 87 98 98 26 34 54 55;23 54 67 87 42 32 33 65 44;43 87 98 44 22 44 33 55 54 ;26 34 76 43 75 83 84 98 23;25 65 87 87 33 32 43 54 63G = 45 56 87 98 98 26 34 54 55 23 54 67 87 42 32 33 65 44 43 87 98 44 22 44 33 55 54 26 34 76 43 75 83 84 98 23 25 65 87 87 33 32 43 54 63请输入窗口大小（象元表示;）：4w

11、= 4L = Columns 1 through 5 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN 1966 1405 2481 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN Columns 6 through 9 NaN NaN NaN NaN NaN NaN NaN NaN 1573 1090 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN请输入经验阈值:1500y = 1500r = 3 3 3c = 3 5 6二Marover算子的点特征提取原程序:clear;G=inp

12、ut(请输入搜索区灰度矩阵：)w=input(请输入窗口大小（象元表示;（如5*5，则输入5）：)s=size(G);if max(s)y);for i=1:size(r) Z(i,:)=r(i),c(i);enddisp(特征点坐标表示（如点（2，3）则表示为2 3）如下：);Z运行实例：已知一灰度距阵（计算中的边界点可视为2）如下：2 2 2 2 22 2 2 2 22 10 2 6 22 2 2 2 2 2 2 2 2 2试用morver算子提取其中的特征点。运行结果： 请输入搜索区灰度矩阵：2 2 2 2 2;2 2 2 2 2;2 10 2 6 2;2 2 2 2 2;2 2 2 2

13、 2G = 2 2 2 2 2 2 2 2 2 2 2 10 2 6 2 2 2 2 2 2 2 2 2 2 2请输入窗口大小（象元表示;（如5*5，则输入5）：3w = 3L = NaN NaN NaN NaN NaN NaN 0 0 0 NaN NaN 128 0 32 NaN NaN 0 0 0 NaN NaN NaN NaN NaN NaN请输入阈值:32y = 32特征点坐标表示（如点（2，3）则表示为2 3）如下：Z = 3 2三Wong-Trinder圆点定位算子原程序:W=input(请输入灰度矩阵:)s=size(W);T=(min(min(W)+mean(mean(W)/2

14、; G=WT;M=sum(sum(G.*W);M=sum(sum(G.*W); I=repmat(1:s(2),s(1),1); J=repmat(1:s(1),1,s(2); x=round(sum(sum(I.*G.*W)/M); y=round(sum(sum(J.*G.*W)/M); disp(目标重心坐标(x,y);disp(x y);运行实例：运用MATLAB软件自动生成一个9*8的距阵，作为目标矩阵。运行结果：;请输入灰度矩阵:fix(rand(9,8)*20)W = 3 13 4 0 15 13 11 1 2 14 8 6 7 8 3 2 13 10 13 8 16 16 3

15、2 1 6 7 5 15 15 8 13 10 2 14 3 7 3 14 9 10 11 7 16 4 17 16 3 17 5 13 8 15 19 15 9 9 0 14 17 18 10 6 2 7 15 8 7 6 17 10 1目标重心坐标(x,y)4 5四移动曲面曲面原程序：clear;D=xlsread(hu.xlsx) %DEM高程数据C=input(请输入待求高程点坐标x y:) %待求高程坐标w=input(请输入拟合点下标(例如输入1 2 3 4 5 6 7):) %按半径大小选取所需点 D1=D(:,1)-C(1) D(:,2)-C(2) D(:,3); %坐标重心化

16、W1=D1(:,1).2 +D1(:,2).2; %求距离平方大小并排序W2,index=sort(W1); for i=1:size(index) W3(i,1:3)=D1(index(i),1:3); %经排序后的坐标end k=0;for i=1:size(D,1) if size(find(w=i),2)=0 %未代入的点即检验点 k=k+1; Q(k,1:3)=W3(i,1:3); endendW4=W3(w,:);n=size(w,2);M=W4(:,1).2 W4(:,1).*W4(:,2) W4(:,2).2 W4(:,1) W4(:,2) ones(n,1) %求MZ=W4(

17、:,3) %求ZP=diag(1./W4(:,1).2+W4(:,2).2) %求权阵PX=(M*P*M)M*P*Z;disp(曲面方程式各系数:);disp( A B C D E F);disp(X);X(6)Zp=X(1).*Q(1:k,1).2+X(2).*Q(1:k,1).*Q(1:k,2)+X(3).*Q(1:k,2).2+X(4).*Q(1:k,1)+X(5).*Q(1:k,2)+X(6); %求实际精度disp(检测点理论高程与实际高程差值:);Zp(:,1)-Q(:,3) %理论与实际差值O=sum(Zp(:,1)-Q(:,3).2 %求精度plot(D(:,1),D(:,2)

18、,*); %画出实际图形grid on;hold on;plot(C(1),C(2),h);ima=sprintf(x-%d)2+(y-%d)2-%d,C(1),C(2),W2(n);ezplot(ima,C(1)-sqrt(W2(n),C(1)+sqrt(W2(n),C(2)-sqrt(W2(n),C(2)+sqrt(W2(n);axis equal;axis square;（2）运行实例已知点的坐标为：点号XYZ110211015210911318310511519410310317510810521610510815711510420811810815911611317101131182

19、2求待定点(110 110)的坐标（3） 运行结果：文件（sun.xls）的内容数据为：D = 102 110 15 109 113 18 105 115 19 103 103 17 108 105 21 105 108 15 115 104 20 118 108 15 116 113 17 113 118 22请输入待求高程点坐标x y:110 110C = 110 110请输入拟合点下标(例如输入1 2 3 4 5 6 7):1:6 7 注释:可以按半径大小序号动态选点插值.w = 1 2 3 4 5 6 7M = 1 -3 9 -1 3 1 4 10 25 -2 -5 1 25 10 4

20、 -5 -2 1 36 18 9 6 3 1 25 -25 25 -5 5 1 25 -30 36 5 -6 1 64 0 0 -8 0 1Z = 18 21 15 17 19 20 15P = 0.1000 0 0 0 0 0 0 0 0.0345 0 0 0 0 0 0 0 0.0345 0 0 0 0 0 0 0 0.0222 0 0 0 0 0 0 0 0.0200 0 0 0 0 0 0 0 0.0164 0 0 0 0 0 0 0 0.0156曲面方程式各系数: A B C D E F -0.0368 0.0395 0.1984 -0.0288 0.1046 15.8433ans = 15.8433检测点理论高程与实际高程差值:ans = -1.7874 7.9072 8.1669O