摄影测量实验报告前方交汇后方交汇.doc

摄影测量实验报告前方交汇后方交汇.doc

《摄影测量实验报告前方交汇后方交汇.doc》由会员分享，可在线阅读，更多相关《摄影测量实验报告前方交汇后方交汇.doc（12页珍藏版）》请在冰豆网上搜索。

摄影测量学

实验报告

学院：

地信院

班级：

测绘0904班

老师：

邹峥嵘

姓名：

张文佳

学号：

0405090921

2011年11月11日

空间后方交会——空间前方交会

程序编程实验

一．实验目的

1、要求掌握运用摄影测量中空间后方交会-空间前方交会求解地面点的空间位置的方法和原理。

2、学会运用空间后方交会的原理，根据所给控制点的地面摄影测量坐标系坐标以及相应的像平面坐标系中的坐标，利用计算机编程语言实现空间后方交会的计算，完成所给像对中两张像片各自的六个外方位元素的求解和精度评定。

3、根据空间后方交会所得的两张像片的内外方位元素，利用同名像点在左右像片上的坐标，利用计算机编程语言前方交会编程，求解其对应的地面点在摄影测量坐标系中的坐标，从而达到通过摄影测量量测地面地理数据的目的。

二．实验仪器

1、计算机

2、MATLAB计算机编程软件

三、实验数据

实验数据实验数据包含四个地面控制点（GCP）的地面摄影测量坐标及在左右像片中的像平面坐标。

此四对坐标运用最小二乘法求解左右像片的外方位元素，即完成了空间后方的过程。

另外还给出了5对地面点在左右像片中的像平面坐标和左右像片的内方位元素。

实验数据如下：

f=150.000mm，x0=0，y0=0

点号

左片

右片

地面摄影测量坐标

x

y

x

y

X

Y

Z

GCP1

16.012

79.963

-73.93

78.706

5083.205

5852.099

527.925

GCP2

88.56

81.134

-5.252

78.184

5780.02

5906.365

571.549

GCP3

13.362

-79.37

-79.122

-78.879

5210.879

4258.446

461.81

GCP4

82.24

-80.027

-9.887

-80.089

5909.264

4314.283

455.484

1

51.758

80.555

-39.953

78.463

2

14.618

-0.231

-76.006

0.036

3

49.88

-0.782

-42.201

-1.022

4

86.14

-1.346

-7.706

-2.112

5

48.035

-79.962

-44.438

-79.736

四、程序设计流程图

1、后方交会

输入GCP的像点坐标x，y

计算旋转矩阵R

计算像点在像空间坐标系中的近似值xj,yj,并组成误差方程的常数项,L=x-xj

计算误差方程的系数项组成系数矩阵A

组成法方程式

计算系数A’A常数项A’L

并求解外方位元素

计算c、w、k、Xs、Ys、Zs改正后的值

确定初始值c=w=k=0，Xs0，Ys0，Zs0

小于

跳出循环，完成迭代计算，精度评定

判断改正数是否小于限差?

大于

此过程完成空间后方交会求解像片的外方位元素，其中改正数小于限差（长度改正数小于0.01m,角度改正数小于0.0003，相当于1’的角度值）为止。

在这个过程中采用迭代计算的方法，是外方位元素逐渐收敛于理论值，每次迭代所得的改正数都应加到上一次的初始值之中。

2、前方交会

输入所需计算点的像平面坐标x,y

根据后方交会所得的旋转矩阵Ra,Rb计算像点在左、右像空间辅助坐标系中的坐标XaYaZa，XbYbZb

计算摄影基线的三个坐标分量BxByBz

计算个点在左右像片中的的投影系数Na,Nb

计算地面所求点在地面摄影测量坐标系中的坐标XYZ

计算完毕

七、实验原理公式

1、后方交会中运用的共线方程数学模型

3、前方交会与后方交会中均用到旋转矩阵进行的坐标转换

4、精度评定中均采用最小二乘准则进行平差计算

5.前方交会的转换

八．实验源程序

1.空间后方交会（以左片为例）

%已知地面摄影测量坐标

Xg=[5083.205,5780.02,5210.879,5909.264];

Yg=[5852.099,5906.365,4258.466,4314.283];

Zg=[527.925,571.549,461.81,455.484];

%对应地面坐标的左片像点坐标

x=[0.016012,0.08856,0.013362,0.08224];

y=[0.079963,0.081134,-0.07937,-0.080027];

%设置初始值

c=0;w=0;k=0;

f=0.15;

Dg=sqrt（（Xg

（1）-Xg

（2））^2+（Yg

（1）-Yg

（2））^2）;

Ds=sqrt（（x

（1）-x

（2））^2+（y

（1）-y

（2））^2）;

p=Dg/Ds;

Xs0=1/4*（Xg

（1）+Xg

（2）+Xg（3）+Xg（4））;

Ys0=1/4*（Yg

（1）+Yg

（2）+Yg（3）+Yg（4））;

Zs0=p*f+（Zg

（1）+Zg

（2）+Zg（3）+Zg（4））/4;

W=0%统计迭代计算次数

%完成迭代计算，检验改正是是否符合要求

while1

%计算旋转矩阵系数

a1=cos（c）*cos（k）-sin（c）*sin（w）*sin（k）;

a2=-cos（c）*sin（k）-sin（c）*sin（w）*cos（k）;

a3=-sin（c）*cos（w）;

b1=cos（w）*sin（k）;

b2=cos（w）*cos（k）;

b3=-sin（w）;

c1=sin（c）*cos（k）+cos（c）*sin（w）*sin（k）;

c2=-sin（c）*cos（k）+cos（c）*sin（w）*cos（k）;

c3=cos（c）*cos（w）;

R=[a1,a2,a3;b1,b2,b3;c1,c2,c3]

forn=1:

1:

4

X（n）=a1*（Xg（n）-Xs0）+b1*（Yg（n）-Ys0）+c1*（Zg（n）-Zs0）;

Y（n）=a2*（Xg（n）-Xs0）+b2*（Yg（n）-Ys0）+c2*（Zg（n）-Zs0）;

Z（n）=a3*（Xg（n）-Xs0）+b3*（Yg（n）-Ys0）+c3*（Zg（n）-Zs0）;

xj（n）=-f*X（n）/Z（n）;%计算近似点坐标

yj（n）=-f*Y（n）/Z（n）;

a11（n）=1/Z（n）*（a1*f+a3*x（n））;%矩阵系数

a12（n）=1/Z（n）*（b1*f+b3*x（n））;

a13（n）=1/Z（n）*（c1*f+c3*x（n））;

a21（n）=1/Z（n）*（a2*f+a3*y（n））;

a22（n）=1/Z（n）*（b2*f+b3*y（n））;

a23（n）=1/Z（n）*（c2*f+c3*y（n））;

a14（n）=y（n）*sin（w）-（x（n）/f*（x（n）*cos（k）-y（n）*sin（k））+f*cos（k））*cos（w）;

a15（n）=-f*sin（k）-x（n）/f*（x（n）*sin（k）+y（n）*cos（k））;

a16（n）=y（n）;

a24（n）=-x（n）*sin（w）-（y（n）/f*（x（n）*cos（k）-y（n）*sin（k））-f*sin（k））*cos（w）;

a25（n）=-f*cos（k）-y（n）/f*（x（n）*sin（k）+y（n）*cos（k））;

a26（n）=-x（n）;

end

A=[a11

（1）,a12

（1）,a13

（1）,a14

（1）,a15

（1）,a16

（1）;

a21

（1）,a22

（1）,a23

（1）,a24

（1）,a25

（1）,a26

（1）;

a11

（2）,a12

（2）,a13

（2）,a14

（2）,a15

（2）,a16

（2）;

a21

（2）,a22

（2）,a23

（2）,a24

（2）,a25

（2）,a26

（2）;

a11（3）,a12（3）,a13（3）,a14（3）,a15（3）,a16（3）;

a21（3）,a22（3）,a23（3）,a24（3）,a25（3）,a26（3）;

a11（4）,a12（4）,a13（4）,a14（4）,a15（4）,a16（4）;

a21（4）,a22（4）,a23（4）,a24（4）,a25（4）,a26（4）];

L=[x

（1）-xj

（1）;y

（1）-yj

（1）;x

（2）-xj

（2）;y

（2）-yj

（2）;x（3）-xj（3）;y（3）-yj（3）;x（4）-xj（4）;y（4）-yj（4）];

Xp=inv（（A'）*A）*（A'）*L

V=A*Xp-L

Xs0=Xs0+Xp（1,1）

Ys0=Ys0+Xp（2,1）

Zs0=Zs0+Xp（3,1）

c=c+Xp（4,1）;w=w+Xp（5,1）;k=k+Xp（6,1）;

W=W+1

%判断收敛条件if（（abs（Xp（1,1））<0.01&&abs（Xp（2,1））<0.01&&abs（Xp（3,1））<0.01&&abs（Xp（4,1））<0.0003&&abs（Xp（5,1））<0.0003&&abs（Xp（6,1））<0.0003））

break

end

end

%精度评定

Qxx=inv（（A'）*A）%外方元素协因素阵

mo=sqrt（（（V）'*V）/（2*8-6））%单位权中误差

mi=mo*sqrt（Qxx）%外方元素改正数中误差

2.前方交会

%后方交会中计算出的改正后外方元素

Xs1=4999.7;Ys1=5000.1;Zs1=2000.0;

Xs2=5896.3;Ys2=5087.9;Zs2=2029.9;

%量测像点左右片坐标

xa=[0.051758,0.014618,0.04988,0.08614,0.048035];

ya=[0.080555,-0.000231,-0.000782,-0.001346,-0.079962];

xb=[-0.039953,-0.076006,-0.042201,-0.007706,-0.044438];

yb=[0.078463,0.000036,-0.001022,-0.002112,-0.079736];

%由外方位线元素计算基线分量Bx,By,Bz

Bx=Xs2-Xs1;

By=Ys2-Ys1;

Bz=Zs2-Zs1;

%由外方位角元素计算像空间辅助坐标X1,Y1,Z1,X2,Y2,Z2

R1=[0.9955,-0.