摄影测量实习空间后方交会
—. 后方交会解算程序。 （1） 原程序：
D=input('请输入地面点坐标:') Y=input('请输入影像坐标点:')/1000 f=input('请输入主距:')/1000
C=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(j,2)/i; end Zs=m*f; for j=1:1000
R(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))*sin(J(3))+cos(J(1))*sin(J(2))*cos(J(3)); R(3,3)=cos(J(1))*cos(J(2)); for n=1:i
Xp(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(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*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)-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; end end
二．运行实例（课本p-44 习题24） 27.已知四队点的影像坐标和地点坐标如下： x(mm) 1 —86.15 y(mm0 —68.99 影像坐标 X(mm) 365.41 地面坐标 Y(mm) 25273.32 Z(mm) 2195.17 2 —53.40 82.21 37631.08 31324.51 728.69 3 —14.78 —76.63 39100.97 24934.98 2386.50 4 10.46 .43 40426. 30319.81 757.31 试计算近似垂直摄影情况下空间后方交会的解。假设内方位元素已知：
f=153.24mm,x0=y0=0.
（3）运行结果：请输入地面点坐标:[365.41 25273.32 2195.17;37631.08 31324.51
728.69;39100.97 24934.98 2386.50;40426. 30319.91 757.31]
D =
1.0e+004 *
3.65 2.5273 0.2195 3.7631 3.1325 0.0729 3.9101 2.4935 0.2387 4.0427 3.0320 0.0757
请输入影像坐标点:[-86.15 -68.99;-53.40 82.21;-14.78 -76.63;10.46 .43] Y =
-0.0862 -0.0690 -0.0534 0.0822 -0.0148 -0.0766 0.0105 0.04
请输入主距:[153.24] f =
0.1532
请输入像主点坐标初始值:[0 0 ] C =
0 0 0
请输入外方位角初始值:[0 0 0] J =
0 0 0 Xs =
3.9795e+004 Ys =
2.7477e+004 Zs =
7.5728e+003 R =
0.9977 0.0675 -0.0675 0.9977 -0.0041 0.0018
0.0039 -0.0021 1.0000
二．Moravec点特征提取算法程序。
（1）原程序
G=input('请输入搜索区灰度矩阵：') w=input('请输入窗口大小（象元表示;）：') s=size(G); if max(s) v1=sum(diff(G(r,c-k:c+k)).^2); v2=sum(diff(diag(G(r-k:r+k,c-k:c+k))).^2); v3=sum(diff(G(r-k:r+k,c)).^2); v4=sum(diff(diag(flipud(G(r-k:r+k,c-k:c+k)))).^2); iv(r-k,c-k)=min([v1 v2 v3 v4]); end end L=nan*ones(s); L(k+1:s(1)-k,k+1:s(2)-k)=iv y=input('请输入经验阈值:') [r,c]=find(L>y) （2）运行实例（自由编辑一组灰度值） [45 56 87 98 98 26 34 55;23 67 87 42 32 33 65 44;43 87 98 44 22 44 33 55 ;26 34 76 43 75 83 84 98 23;25 65 87 87 33 32 43 63] （3）运行结果：请输入搜索区灰度矩阵：[45 56 87 98 98 26 34 55;23 67 87 42 32 33 65 44;43 87 98 44 22 44 33 55 ;26 34 76 43 75 83 84 98 23;25 65 87 87 33 32 43 63] G = 45 56 87 98 98 26 34 55 23 67 87 42 32 33 65 44 43 87 98 44 22 44 33 55 26 34 76 43 75 25 65 87 87 33 请输入窗口大小（象元表示;）：4 w = 4 L = Columns 1 through 5 NaN NaN NaN NaN NaN NaN NaN NaN 1966 NaN NaN NaN NaN NaN NaN 83 84 32 43 NaN NaN 1405 NaN NaN 98 NaN NaN 2481 NaN NaN 23 63 Columns 6 through 9 NaN NaN NaN NaN NaN NaN NaN NaN 1573 1090 NaN NaN NaN NaN NaN NaN NaN NaN 请输入经验阈值:1500 y = 1500 r = 3 3 3 c = 3 5 6 NaN NaN 二．Marover算子的点特征提取 原程序: clear; G=input('请输入搜索区灰度矩阵：') w=input('请输入窗口大小（象元表示;（如5*5，则输入5）：') s=size(G); if max(s) v1=sum(diff(G(r,c-k:c+k)).^2); v2=sum(diff(diag(G(r-k:r+k,c-k:c+k))).^2); v3=sum(diff(G(r-k:r+k,c)).^2); v4=sum(diff(diag(flipud(G(r-k:r+k,c-k:c+k)))).^2); iv(r-k,c-k)=min([v1 v2 v3 v4]); end end L=nan*ones(s); L(k+1:s(1)-k,k+1:s(2)-k)=iv y=input('请输入阈值:') [r,c]=find(L>y); for i=1:size(r) Z(i,:)=[r(i),c(i)]; end disp('特征点坐标表示（如点（2，3）则表示为2 3）如下：'); Z 运行实例： 已知一灰度距阵（计算中的边界点可视为2）如下： 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 试用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 2] G = 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）：3 w = 3 L = 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 请输入阈值:32 y = 32 特征点坐标表示（如点（2，3）则表示为2 3）如下： Z = 3 2 三Wong-Trinder圆点定位算子 原程序: W=input('请输入灰度矩阵:') s=size(W); T=(min(min(W))+mean(mean(W)))/2; G=W>T; 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 2 1 6 7 5 15 10 2 14 3 7 10 11 7 16 4 17 5 13 8 15 9 0 14 17 18 7 15 8 7 6 目标重心坐标(x,y) 4 5 四移动曲面曲面 原程序：clear; D=xlsread('hu.xlsx') C=input('请输入待求高程点坐标[x y]:') w=input('请输入拟合点下标(例如输入[1 2 3 4 5 6 7]):') 点 D1=[D(:,1)-C(1) D(:,2)-C(2) D(:,3)]; 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); end end W4=W3(w,:); 15 8 13 3 14 9 16 3 19 15 9 10 6 2 17 10 1 %DEM高程数据 %待求高程坐标 %按半径大小选取所需 %坐标重心化 %求距离平方大小并排 %经排序后的坐标 %未代入的点即检验点 17 n=size(w,2); M=[W4(:,1).^2 W4(:,1).*W4(:,2) W4(:,2).^2 W4(:,1) W4(:,2) ones(n,1)] %求M Z=[W4(:,3)] %求Z P=diag(1./[W4(:,1).^2+W4(:,2).^2]) %求权阵P X=(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),'*'); %画出实际图形 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）运行实例 已知点的坐标为： 点号 1 2 3 4 5 6 7 8 9 10 求待定点(110 110)的坐标 （3） 运行结果： X 102 109 105 103 108 105 115 118 116 113 Y 110 113 115 103 105 108 104 108 113 118 Z 15 18 19 17 21 15 20 15 17 22 文件（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 110] C = 110 110 请输入拟合点下标(例如输入[1 2 3 4 5 6 7]):[1:6 7] 注释:可以按半径大小序号动态选点插值. w = 1 2 3 4 5 6 7 M = 1 -3 9 -1 3 1 4 10 25 -2 -5 1 25 10 4 -5 -2 1 36 18 9 6 3 1 25 -25 25 -5 5 1 25 -30 36 5 -6 1 0 0 -8 0 1 Z = 18 21 15 17 19 20 15 P = 0.1000 0 0 0 0.0345 0 0 0 0.0345 0 0 0 0 0 0 0 0 0 0 0 0 曲面方程式各系数: A B C -0.0368 0.0395 0.1984 ans = 15.8433 检测点理论高程与实际高程差值: ans = -1.7874 7.9072 0 0 0 0 0 0 0.0222 0 0 0.0200 0 0 0 0 D E -0.0288 0.1046 0 0 0 0 0 0 0 0 0 0 0.01 0 0 0.0156 F 15.8433 8.1669 O = 44.1392 图像为： 再次运行检验上次的结果： 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 110] C = 110 110 请输入拟合点下标(例如输入[1 2 3 4 5 6 7]):[1:4 7 8] w = 1 2 3 4 7 8 M = 1 -3 9 -1 3 1 4 10 25 -2 -5 1 25 10 4 -5 -2 1 36 18 9 6 3 1 0 0 -8 0 1 -16 4 8 -2 1 Z = 18 21 15 17 15 15 P = 0.1000 0 0 0 0 0 0 0.0345 0 0 0 0 0 0 0.0345 0 0 0 0 0 0 0.0222 0 0 0 0 0 0 0.0156 0 0 0 0 0 0 0.0147 曲面方程式各系数: A B C D E F 0.0068 -0.0313 0.3435 -0.0829 0.2740 13.9030 ans = 13.9030 检测点理论高程与实际高程差值: ans = 6.2263 5.3176 15.1397 11.1950 O = 105.3957 结果分析： 待求点 (110,110)的高程为: 13.9030 实际精度为: 105.3957 用圆圈圈住的点为检测点. 可见拟合结果和拟合点的个数及拟合点的与待求点的位置有很大关系. 因篇幅问题不能全部显示，请点此查看更多更全内容
Copyright © 2019- sarr.cn 版权所有 赣ICP备2024042794号-1
违法及侵权请联系：TEL:199 1889 7713 E-MAIL:2724546146@qq.com
本站由北京市万商天勤律师事务所王兴未律师提供法律服务