摄影测量学后方交会matlab实习报告

摄影测量原理

单张影像后方交会实习

目录

一实习目的 (3)

二实习原理 (3)

1. 间接平差 (3)

2. 共线方程 (3)

3. 单向空间后方交会 (4)

三计算流程 (4)

1. 求解步骤 (4)

2.计算机框图 (4)

四程序实现 (5)

五结果分析 (6)

1.外方位元素 (6)

2.误差 (6)

3.旋转矩阵R (7)

六实习体会 (7)

1. 平台的选择 (7)

2.问题的解决 (7)

3.心得体会 (8)

七代码展示 (8)

一实习目的

为了增强同学们对后方交会公式的理解，培养同学们对迭代循环编程的熟悉感，本次摄影测量课间实习内容定为用C语言或其他程序编写单片空间后方交会程序，最终输出像点坐标、地面坐标、单位权中误差、外方位元素及其精度。

已知四对点的影像坐标和地面坐标如下。内方位元素fk=153.24mm，x0=y0=0。

本次实习，我使用了matlab2014进行后方交会程序实现。结果与参考答案一致，精度良好。

二实习原理

题干中有四个控制点在地面摄影测量坐标系中的坐标和对应的像点坐标，由此可列出8个误差方程，存在2个多余观测（n=2）。故可利用间接平差的最小二乘法则求解。

由于共线方程是非线性函数模型，为了方便计算，需要将其“线性化”。但如果仅取泰勒级数展开式的一次项，未知数的近似值改正是不精确的。因此必须采用迭代趋近法计算，直到外方位元素的改正值小于限差。

1.间接平差

间接平差为平差计算最常用的方法。在确定多个未知量的最或然值时，选择它们之间不存在任何条件关系的独立量作为未知量组成用未知量表达测量的函数关系、列出误差方程式，按最小二乘法原理求得未知量的最或然值的平差方法。

在一个间接平差问题中，当所选的独立参数X个数与必要观测值t个数相等时，可将每个观测值表达成这t个参数的函数，组成观测方程。

函数模型为：L = BX + d。

2.共线方程

共线方程是中心投影构像的数学基础，也是各种摄影测量处理方法的重要理论基础。

式中：

x，y 为像点的像平面坐标；

x0，y0，f 为影像的内方位元素；

XS，YS，ZS 为摄站点的物方空间坐标；

XA，YA，ZA 为物方点的物方空间坐标；

ai，bi，ci （i = 1,2,3）为影像的3 个外方位角元素组成的9 个方向余弦。

3.单向空间后方交会

利用至少三个已知地面控制点的坐标A、B、Z与其影像上对应的三个像点的影像坐标a、b、c，根据共线方程，反求该像点的外方位元素。这种解算方法以单张相片为基础，亦称单像空间后方交会。

三计算流程

1.求解步骤

（1）获取已知数据。

（2）量测控制点的像点坐标并进行必要的影像坐标系统误差改正，得到像点坐标。（3）确定未知数的初始值。在竖直航空摄影且地面控制点大体对称分布的情况下，可按如下方法确定初始值：

Z S0 = H = m*f

X S0 = 1/n*(X1+X2+…+X n)

Y S0 = 1/n*(Y1+Y2+…+Y n)

q0 = w0 = 0

k0可在航迹图上找出，或根据控制点坐标通过坐标正反变换求出。

（4）利用角元素近似值计算方向余弦值，组成旋转矩阵R。

（5）利用未知数的近似值按共线方程条件式计算控制点像点坐标的近似值x0，y0。（6）逐点计算误差方程式的系数和常数项，组成误差方程式。

（7）解求外方位元素改正数，并与相应的近似值求和，得到外方位元素新的近似值。

（8）规定限差限制迭代次数，检查计算是否收敛。

2.计算机框图

四程序实现

1.定义各种参数名并赋值。

由于题目是近似垂直观测，角元素q，w，k的初值我统一赋了0。从最终结果也可以看出三个值都非常接近0。

2.读取坐标

新建名为“坐标.txt”的文本文件，按顺序输入每个点的影像坐标和地面坐标。使用matlab读取文件语句textread将x，y，X，Y，Z分别存储到对应矩阵中。

3.从计算R阵开始，直到外方位元素改正，需要进行数次迭代循环。结束循环的条件为线元素最终改正值绝对值小于limit1（0.001），角元素最终改正值绝对值小于limit1（0.001），两者必须同时成立。

4.进入循环后，先按照书上公式，用程序中外方位元素赋的初值计算旋转矩阵R。

5.按照公式计算共线方程。共线方程左侧应为x-x0/y-y0，但本题中x0、y0都为0，故略去一步，直接记录左侧数值为x/y。存储在矩阵的顺序与c矩阵相同，方便二者作差求出L阵。

6.利用一个i（1~4）的循环，计算出系数阵A。在此之前需要计算Zbar。

7.计算常数项L，进而推导出外方位元素改正值dx。

8.改正外方位元素大小，判断是否满足循环条件。若不满足则输出结果。

9.循环结束后，分别输出外方位元素值、外方位元素误差和R阵，保存在txt文件中。

五结果分析

1.外方位元素结果如下；迭代次数为5次。

结果与参考答案一致，迭代次数也符合要求，说明外方位元素初值选取的比较好。本张像片近似垂直摄影。

2.外方位元素误差都比较小，同时角元素精度远高于线元素精度。

3.旋转矩阵R各元素大小四舍五入后与参考答案一致。

六实习体会

七代码展示

%单幅影像后方交会

%线元素10^-3,角元素10^-6

close all; clear all; clc;

disp('后方交会求解')%屏幕输出

[x,y,X,Y,Z]=textread('坐标.txt','%f %f %f %f %f');%读取文件

x=x/1000;y=y/1000;%修正单位

fk=153.24/1000;%主距修正单位

m=50000;%比例尺系数

limit1=0.001;limit2=0.000001;%改正值精度

Xs=sum(X)/4;Ys=sum(Y)/4;Zs=fk*m+sum(Z)/4;%待定参数初始值

q=0;w=0;k=0;%外方位角元素

R=zeros(3,3);%旋转矩阵

dx=ones(6,1);%改正值

H=m*fk;%航高

round=0;%迭代次数

c=[x(1);y(1);x(2);y(2);x(3);y(3);x(4);y(4)];%方便L计算

PI=3.1415926;

%循环条件：线元素改正值<=limit1且角元素改正值<=limit2

while(abs(dx(1,1))>limit1||abs(dx(2,1))>limit1||abs(dx(3,1))>limit1||a bs(dx(4,1))>limit2||abs(dx(5,1))>limit2||abs(dx(6,1))>limit2)

%迭代计算旋转矩阵R

R(1,1)=cos(q)*cos(k)-sin(q)*sin(w)*sin(k);

R(1,2)=-cos(q)*sin(k)-sin(q)*sin(w)*cos(k);

R(1,3)=-sin(q)*cos(w);

R(2,1)=cos(w)*sin(k);

R(2,2)=cos(w)*cos(k);

R(2,3)=-sin(w);

R(3,1)=sin(q)*cos(k)+cos(q)*sin(w)*sin(k);

R(3,2)=-sin(q)*sin(k)+cos(q)*sin(w)*cos(k);

R(3,3)=cos(q)*cos(w);

xy_0=zeros(8,1);%共线方程。奇数为x-x0,偶数为y-y0

xy_0(1,1)=-fk*(R(1,1)*(X(1)-Xs)+R(2,1)*(Y(1)-Ys)+R(3,1)*(Z(1)-Zs))/(R( 1,3)*(X(1)-Xs)+R(2,3)*(Y(1)-Ys)+R(3,3)*(Z(1)-Zs));

xy_0(2,1)=-fk*(R(1,2)*(X(1)-Xs)+R(2,2)*(Y(1)-Ys)+R(3,2)*(Z(1)-Zs))/(R( 1,3)*(X(1)-Xs)+R(2,3)*(Y(1)-Ys)+R(3,3)*(Z(1)-Zs));

xy_0(3,1)=-fk*(R(1,1)*(X(2)-Xs)+R(2,1)*(Y(2)-Ys)+R(3,1)*(Z(2)-Zs))/(R( 1,3)*(X(2)-Xs)+R(2,3)*(Y(2)-Ys)+R(3,3)*(Z(2)-Zs));

xy_0(4,1)=-fk*(R(1,2)*(X(2)-Xs)+R(2,2)*(Y(2)-Ys)+R(3,2)*(Z(2)-Zs))/(R( 1,3)*(X(2)-Xs)+R(2,3)*(Y(2)-Ys)+R(3,3)*(Z(2)-Zs));

xy_0(5,1)=-fk*(R(1,1)*(X(3)-Xs)+R(2,1)*(Y(3)-Ys)+R(3,1)*(Z(3)-Zs))/(R( 1,3)*(X(3)-Xs)+R(2,3)*(Y(3)-Ys)+R(3,3)*(Z(3)-Zs));

xy_0(6,1)=-fk*(R(1,2)*(X(3)-Xs)+R(2,2)*(Y(3)-Ys)+R(3,2)*(Z(3)-Zs))/(R( 1,3)*(X(3)-Xs)+R(2,3)*(Y(3)-Ys)+R(3,3)*(Z(3)-Zs));

xy_0(7,1)=-fk*(R(1,1)*(X(4)-Xs)+R(2,1)*(Y(4)-Ys)+R(3,1)*(Z(4)-Zs))/(R( 1,3)*(X(4)-Xs)+R(2,3)*(Y(4)-Ys)+R(3,3)*(Z(4)-Zs));

xy_0(8,1)=-fk*(R(1,2)*(X(4)-Xs)+R(2,2)*(Y(4)-Ys)+R(3,2)*(Z(4)-Zs))/(R( 1,3)*(X(4)-Xs)+R(2,3)*(Y(4)-Ys)+R(3,3)*(Z(4)-Zs));

A=zeros(8,6);%系数阵

for i=1:4

Xbar=R(1,1)*(X(i)-Xs)+R(2,1)*(Y(i)-Ys)+R(3,1)*(Z(i)-Zs);

Ybar=R(1,2)*(X(i)-Xs)+R(2,2)*(Y(i)-Ys)+R(3,2)*(Z(i)-Zs);

Zbar=R(1,3)*(X(i)-Xs)+R(2,3)*(Y(i)-Ys)+R(3,3)*(Z(i)-Zs);

A(2*i-1,1)=1/Zbar*(R(1,1)*fk+R(1,3)*xy_0(2*i-1,1));

A(2*i-1,2)=1/Zbar*(R(2,1)*fk+R(2,3)*xy_0(2*i-1,1));

A(2*i-1,3)=1/Zbar*(R(3,1)*fk+R(3,3)*xy_0(2*i-1,1));

A(2*i,1)=1/Zbar*(R(1,2)*fk+R(1,3)*xy_0(2*i,1));

A(2*i,2)=1/Zbar*(R(2,2)*fk+R(2,3)*xy_0(2*i,1));

A(2*i,3)=1/Zbar*(R(3,2)*fk+R(3,3)*xy_0(2*i,1));

A(2*i-1,4)=xy_0(2*i,1).*sin(w)-(xy_0(2*i-1,1)./fk.*(xy_0(2*i-1,1).*cos (k)-xy_0(2*i,1).*sin(k))+fk*cos(k))*cos(w);

A(2*i-1,5)=-fk*sin(k)-xy_0(2*i-1,1)./fk.*(xy_0(2*i-1,1).*sin(k)+xy_0(2 *i,1).*cos(k));

A(2*i-1,6)=xy_0(2*i,1);

A(2*i,4)=-xy_0(2*i-1,1).*sin(w)-(xy_0(2*i,1)./fk.*(xy_0(2*i-1,1).*cos( k)-xy_0(2*i,1).*sin(k))-fk*sin(k))*cos(w);

A(2*i,5)=-fk*cos(k)-xy_0(2*i,1)./fk.*(xy_0(2*i-1,1).*sin(k)+xy_0(2*i,1 ).*cos(k));

A(2*i,6)=-xy_0(i*2-1,1);

end

round=round+1;

L=c-xy_0;

dx=((A')*A)\(A')*L;

%改正外方位元素

Xs=Xs+dx(1,1);Ys=Ys+dx(2,1);Zs=Zs+dx(3,1);

q=q+dx(4,1);w=w+dx(5,1);k=k+dx(6,1);

end

%写入文件

fid1=fopen('后方交会.txt','w');%输出外方位元素

fprintf(fid1,'Xs = %f\r\n',Xs);

fprintf(fid1,'Ys = %f\r\n',Ys);

fprintf(fid1,'Zs = %f\r\n\r\n',Zs);

fprintf(fid1,'q = %f\r\n',q);

fprintf(fid1,'w = %f\r\n',w);

fprintf(fid1,'k = %f\r\n\r\n',k);

fprintf(fid1,'round = %d\r\n',round);

fclose(fid1);

fid2=fopen('R阵.txt','w');%输出矩阵R

fprintf(fid2,'\r\n');

for m=1:3

for n=1:3

fprintf(fid2,'%f ',R(m,n));

end

fprintf(fid2,'\r\n');

end

fclose(fid2);

fid3=fopen('精度评定.txt','w');%输出中误差

fprintf(fid3,'dXs = ±%f（mm）\r\n',1000*abs(dx(1,1)));

fprintf(fid3,'dYs = ±%f（mm）\r\n',1000*abs(dx(2,1)));

fprintf(fid3,'dZs = ±%f（mm）\r\n',1000*abs(dx(3,1)));

fprintf(fid3,'\r\ndq = ±%f（″）\r\n',abs(dx(4,1))*180*3600/PI);

fprintf(fid3,'dw = ±%f（″）\r\n',abs(dx(5,1))*180*3600/PI); fprintf(fid3,'dk = ±%f（″）\r\n',abs(dx(6,1))*180*3600/PI); fclose(fid3);