前方交会实验报告

实验一

1 实验任务

理解摄影测量中核心模型-共线方程作用，掌握航空影像中重要的点线面的透视关系以及物方与像方之间的解析关系，单幅影像上像点坐标与相应地面点坐标之间的关系。通过编程实现外方位元素的求解，提升编程能力。

2 理论模型与方法

单张像片的空间后方交会的基本思想：以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应的像坐标量测值处发，根据共线条件方程，解求该影像在航空摄影时刻的元素XS,YS,ZS,,,。 （1）共线方程

xx0fa1(XAXS)b1(YAYS)c1(ZAZS)a3(XAXS)b3(YAYS)c3(ZAZS)a(XXS)b2(YAYS)c2(ZAZS)yy0f2Aa3(XAXS)b3(YAYS)c3(ZAZS)（2）旋转矩阵R

a1coscossinsinsina2cossinsinsincosa3sincosb1cossinb2coscosb3sinc1sincoscossinsinc2sinsincossincosc3coscos（3）经过线性化，得到x，y的误差方程式

vxa11Xsa12Ysa13Zsa14a15a16(x)xvya21Xsa22Ysa23Zsa24a25a26(x)y矩阵形式如下： VAXL 系数方程

a12aA11a21a22

a13a23a14a24a15a25a16 a26

其中：

a11[a1fa3x]/Za12[b1fb3x]/Za13[c1fc3x]/Za21[a2fa3y]/Za22[b2fb3y]/Za23[c2fc3y]/Za14ysin{x/f[xcosysin]fcos}cosa15fsinx/f[xsinycos]a16ya24xsin[y/f(xcosysin)fsin]cosa25fcosy/f(xsinycos)a26xTL(x)x,(y)y

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x0fX/Zyy0fY/ZXXSYS

（4）由最小二乘间接平差原理可得：

ZS

X(ATA)1ATL

XSXSXSXSYSYSYSYSZSZSZSZS 0120120120120120123 程序设计

本地方仅列出核心代码： %确定初值 x0=0;

y0=0; f=153.24;

m=sqrt(((x(1)-x(2))^2+(y(1)-y(2))^2))/(sqrt(((X(1)-X(2))^2+(Y(1)-Y(2))^2))); Zs=f/m;%3个线性元素 Xs=mean(X); Ys=mean(Y);

aa=0; %3个外方位角元素初值 ww=0; kk=0;

cx=zeros(6,1);

p=0.1/2062.806247096363;%将0.1秒限差化为弧度 aa1=1; ww1=1; kk1=1; k=0;

while abs(aa-aa1)>p||abs(ww-ww1)>p||abs(kk-kk1)>p aa1=aa;%赋值 ww1=ww; kk1=kk;

%计算旋转矩阵

a1=cos(aa)*cos(kk)-sin(aa)*sin(ww)*sin(kk); a2=-cos(aa)*sin(kk)-sin(aa)*sin(ww)*cos(kk); a3=-sin(aa)*cos(ww); b1=cos(ww)*sin(kk); b2=cos(ww)*cos(kk); b3=-sin(ww);

c1=sin(aa)*cos(kk)+cos(aa)*sin(ww)*sin(kk); c2=-sin(aa)*sin(kk)+cos(aa)*sin(ww)*cos(kk); c3=cos(aa)*cos(ww);

R=[a1 a2 a3;b1 b2 b3;c1 c2 c3]; %计算误差方程系数 for i=1:d1

%计算近似值

XX=a1*(X(i)-Xs)+b1*(Y(i)-Ys)+c1*(Z(i)-Zs); YY=a2*(X(i)-Xs)+b2*(Y(i)-Ys)+c2*(Z(i)-Zs); ZZ=a3*(X(i)-Xs)+b3*(Y(i)-Ys)+c3*(Z(i)-Zs); a11=(1/ZZ)*(a1*f+a3*(x(i)-x0)); a12=(1/ZZ)*(b1*f+b3*(x(i)-x0)); a13=(1/ZZ)*(c1*f+c3*(x(i)-x0));

a14=(y(i)-y0)*sin(ww)-(((x(i)-x0)/f)*((x(i)-x0)*cos(kk)-(y(i)-y0)*sin(kk))+f*cos(kk))*cos(ww); a15=-f*sin(kk)-((x(i)-x0)/f)*((x(i)-x0)*sin(kk)+(y(i)-y0)*cos(kk)); a16=(y(i)-y0);

a21=(1/ZZ)*(a2*f+a3*(y(i)-y0));

a22=(1/ZZ)*(b2*f+b3*(y(i)-y0)); a23=(1/ZZ)*(c2*f+c3*(y(i)-y0));

a24=-(x(i)-x0)*sin(ww)-(((y(i)-y0)/f)*((x(i)-x0)*cos(kk)-(y(i)-y0)*sin(kk))-f*sin(kk))*cos(ww); a25=-f*cos(kk)-((y(i)-y0)/f)*((x(i)-x0)*sin(kk)+(y(i)-y0)*cos(kk)); a26=-(x(i)-x0);

A(2*i-1,:)=[a11 a12 a13 a14 a15 a16]; A(2*i,:)=[a21 a22 a23 a24 a25 a26]; l(2*i-1,:)=x(i)-(x0-f*XX/ZZ); l(2*i,:)=y(i)-(y0-f*YY/ZZ); end

cx1=inv(A'*A)*A'*l;%最小二乘平差 Xs=Xs+cx1(1); Ys=Ys+cx1(2); Zs=Zs+cx1(3); aa=aa+cx1(4); ww=ww+cx1(5); kk=kk+cx1(6);

k=k+1;%统计循环次数 end

V==A*cx1-l;

sigma=sqrt(V'*V/11);

4 结论与体会

本次实验使用了一种相对比较简单的MATLAB语言来编写外方位元素求解，总体来说，过程比较简单，但是在数据检验方面经常不不符，追其根源是在代码书写过程中错误太多，导致结果不对。相对于Ｃ语言而言，使用MATLAB编写代码较为简单。

实验二

1 实验任务

理解根据立体相对左右倆影像的内、外方位元素和同名像点影像坐标量测值来确定该点的物方空间坐标的过程。学会利用计算机语言实验立体相对的前方交会过程，求解地面点坐标。

2 理论模型与方法

利用点投影系数法实现一对相片的空间前方交会。

NX1BXNX2NY1BYNY2 NZ1BZNZ2 其中：

NBXZ2BZX2X1Z2Z1X1BZBZX1N'X1X1Z2Z1X1

N和N分别为左像点和右像点投影到地面的系数。在相对定向中的上式中元素都是根据像点坐和相对元素来计算，用左右影像的外方位元素来计算，则

BXXS2XS1BYYS2YS1BZZS2ZS1由左右影像的外方位元素1,1,1和2,2,2计算相应的正交矩阵R1和R2，则

X1x1X2x2 Y1R1y1，Y2R2y2

ZzZz1122任一地面点坐标为：

XXS1NX1XS1BXNX2YYS1NY1YS1BY1NY1NYBY 2ZZS1NZ1XS1BXNZ23 程序设计

%内方位参数 X10=0; Y10=0;

f=120.000000/1000;

%左相片外方位参数 Xs1=547350.2538; Ys1=4625202.1458; Zs1=18.1872; aa1=-0.002580; ww1=-0.003753; kk1=0.217876;

%右相片外方位参数 Xs2=547811.4390; Ys2=4625350.1678; Zs2=1886.4757; aa2=-0.008686; ww2=0.002983; kk2=0.212013; %左相片旋转矩阵

a11=cos(aa1)*cos(kk1)-sin(aa1)*sin(ww1)*sin(kk1); a12=-cos(aa1)*sin(kk1)-sin(aa1)*sin(ww1)*cos(kk1); a13=-sin(aa1)*cos(ww1); b11=cos(ww1)*sin(kk1); b12=cos(ww1)*cos(kk1); b13=-sin(ww1);

c11=sin(aa1)*cos(kk1)+cos(aa1)*sin(ww1)*sin(kk1); c12=-sin(aa1)*sin(kk1)+cos(aa1)*sin(ww1)*cos(kk1); c13=cos(aa1)*cos(ww1);

R1=[a11 a12 a13;b11 b12 b13;c11 c12 c13]; %右相片旋转矩阵

a21=cos(aa2)*cos(kk2)-sin(aa2)*sin(ww2)*sin(kk2); a22=-cos(aa2)*sin(kk2)-sin(aa2)*sin(ww2)*cos(kk2); a23=-sin(aa2)*cos(ww2); b21=cos(ww2)*sin(kk2); b22=cos(ww2)*cos(kk2); b23=-sin(ww2);

c21=sin(aa2)*cos(kk2)+cos(aa2)*sin(ww2)*sin(kk2); c22=-sin(aa2)*sin(kk2)+cos(aa2)*sin(ww2)*cos(kk2); c23=cos(aa2)*cos(ww2);

R2=[a21 a22 a23;b21 b22 b23;c21 c22 c23]; %以左相片投影系数为原点

%利用点投影系数法进行前方交会 Bx=(Xs2-Xs1); By=(Ys2-Ys1); Bz=(Zs2-Zs1);

%根据像坐标计算像空间辅助坐标系坐标 for i=1:d3 %左相片

A1=R1*[x1(i);y1(i);-f]; X1=A1(1,1); Y1=A1(2,1); Z1=A1(3,1); %右相片

A2=R2*[x2(i);y2(i);-f]; X2=A2(1,1); Y2=A2(2,1); Z2=A2(3,1); %点投影系数

N=(Bx*Z2-Bz*X2)/(X1*Z2-Z1*X2); N1=(Bx*Z1-Bz*X1)/(X1*Z2-Z1*X2); X(i)=Xs1+N*X1;

Y(i)=Ys1+0.5*(N*Y1+N1*Y2+By); Z(i)=Zs1+Bz+N1*Z2; end

4 结论与体会