单张影像空间后方交会程序

1、摄影测量实习报告单张影像空间后方交会程序实习时间 2013.5.20-2013.5.24学生班级10测绘（1）班学生姓名 房新明学生学号1072143138所在院系矿业工程学院指导老师张会战邵亚琴一.实习目的1 .深入理解单张影像空间后方交会的原理，体会在有多余观测情况下，用 最小二乘平差方法编程实现解求影像外方位元素的过程。2 .利用visual c+编写一个完整的单张影像空间后方交会程序，通过对提 供的试验数据进行计算，输出像片的外方位元素并进行评定精度。3.通过编写程序实现单张影像空间后方交会计算，掌握非线性方程线性化 的过程、相应数据读入与存储的方法以及迭代计算的特点，巩固各类基础课程

2、及 计算机课程的学习内容，培养上机调试程序的动手能力，通过对实验结果的分析， 增强综合运用所学知识解决专业实际问题的能力。4、掌握空间后方交会的定义和实现算法(1)定义：空间后方交会是以单幅影像为基础，从该影像所覆盖地面范围内 若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程， 解求该影像在航空摄影时刻的外方位元素 xs, ys, zs,小，k。(2)算法：由于每一对像方和物方共腕点可列出2个方程，因此若有3个已知地面坐标的控制点，则可列出 6个方程，解求6个外方位元素的改正数 xs, ys, azs, 小， k。实际应用中为了提高解算精度，常有多余观 测方程，通常是在影像

3、的四个角上选取 4个或均匀地选择更多的地面控制点，因 而要用最小二乘平差方法进行计算。5、了解空间后方交会的基本过程(1)获取已知数据。从摄影资料中查取影像比例尺1/m,平均摄影距离(航空摄影的航高)、内方位元素x0, y0, f;获取控制点的空间坐标xt, yt, zt(2)量测控制点的像点坐标并进行必要的影像坐标系统误差改正，得到像点 坐标。(3)确定未知数的初始值。单像空间后方交会必须给出待定参数的初始值， 在竖直航空摄影且地面控制点大体对称分布的情况下，xs0和ys0为均值，zs0为航高，小、k的初值都设为00或者k的初值可在航迹图上找出或根据控 制点坐标通过坐标正反变换求出。(4)计

4、算旋转矩阵ro利用角元素近似值计算方向余弦值，组成 r阵。(5)逐点计算像点坐标的近似值。利用未知数的近似值按共线条件式计算控 制点像点坐标的近似值(x), (y)。(6)逐点计算误差方程式的系数和常数项，组成误差方程式。(7)计算法方程的系数矩阵ata与常数项atl,组成法方程式。(8)解求外方位元素。根据法方程，解求外方位元素改正数，并与相应的近 似值求和，得到外方位元素新的近似值。(9)检查计算是否收敛。将所求得的外方位元素的改正数与规定的限差比较， 通常对小，k的改正数4小， k给予限差，通常为0.1,当3个 改正数均小于0.1时，迭代结束。否则用新的近似值重复(4)(8)步骤的计 算

5、，直到满足要求为止。二、实习要求1、认真复习单张影像空间后方交会的原理及有关内容；2、编写程序，计算左片号23和右片号24的外方位元素并进行精度评定；3、验证数据部分外方位元素和旋转矩阵 r已经给出，请验证自己编写的程序是 否正确；三、实习环境1、硬件环境：windows操作系统2、软件环境：vc+ 四、实习原理1、以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面 坐标和相应点的像坐标量测值出发，根据共线条件方程，求解该影象在航空摄影时刻的像片外方位元素 xs, ys , zs，/，k共线条件方程如下： x-x0=-f*a1(x-xs)+b1(y-ys)+c1(z-zs)/a3(

6、x-xs)+b3(y-ys)+c3(z-zs) y-y0=-f*a2(x-xs)+b2(y-ys)+c2(z-zs)/a3(x-xs)+b3(y-ys)+c3(z-zs)x,y为像点的像平面坐标；xs,ys，伪影像的外方位元素；xs, ys ,zs为摄站点的物方空间坐标；x,y,z为物方点的物方空间坐标；旋转矩阵为r ;由于此共线条件方程是非线性方程，先对其进行线性化，像点观测值一般视为等权，即p=i;矩阵形式：v=ax-l , p=i;通过间接平差，为提高精度，增加多余观测方程，根据最小二乘平差原理，可计 算出外方位元素的改正数。经过迭代计算，每次迭代用未知数的近似值与上次迭 代计算的改正数

7、之和作为新的近似值，重复计算，求出新的改正数，这样反复趋近，直到改正数小于某个限值为止。2、精度评定mi=m0*，qii 其中 m0=v【vm / (2n-6)3、公式推导一一 .干:0 干：0 千：一 一0fx(xs,ys,zs,q%m)=%(xs-xs)+e(ys-ys)+=hzs-zs) + xsyszs0 0 0二(：-0) ( - - 0) (1 -1 0) fx(x；,y：,z0j0,-0j0) cto干j 0mxs,yszsfixlxxs-xs) +0,fy. ys汗0(ys-ys0) 7zs(zs-z：)干0 c ：fc ：:fy (:-* y ( , - .0) y (.0)

8、 f(xo y0 z ： ,00):().：. ()() y( s, s，s，, )同理fy,xsdxs千ydys-ys千y dzs:zs干y0干：汗;0d d df；/的所以干x0dxs 二 dys，s干十.:z0-dzsffd :清生d.，小0emfy0 dxsxs以dys-ys干y0-zsdzs千y0：fy0yd- fy0 =01、基本原理共线条件线性方程式为:f dxs tdzs （ f2x- ）d4 - ?do + ydm - （ x - x计）=0fyxyy2dys 工dzs - - d -（ f 一 ）d - xd -（ y - y计）=0 zzff误差方程式为:f2= dxs

9、+工dzs -（ f + ）d - do + yd _（ x _ x计）=vx z zf f工 dys = dzs -2；d4-（f+y ）d。- xdm-（ y - y计）=vy计算改正值：dxs,dys,dzs,d* ,ds计计算改后的外方位元素：xs xs dxs ys ys dys zs zs dzs二 d . 二 d . 二 d六、 程 序 过 程 框 图开始七.实习数据1、模拟像片一对：左片23号 右片24号2、相片比例尺：1/300003、航摄机主距：f=150mm4、每张相片有四个控制点5、各片像点坐标及地面坐标片号占 八、号像点坐标（mm）地面坐标（m）xyxyz231-91

10、.596-74.859100000.0137500.011.003-94.23081.446100000.0142500.0036.00795.207-75.512106000.0137500.0042.00996.79783.077106000.0142500.0056.00244-102.695-79.618103000.0137500.0090.006-99.90481.754103000.0142500.0031.001086.890-77.540109000.0137500.007.001288.90476.257109000.0142500.005.00八.验证数据已知四对点的影像

11、坐标和地面坐标及航摄仪内方位元素f=153.42mm, x0=y0=0影像坐标地面坐标x(mm)y(mm)x (mm)y (mm)z (mm)1-86.15-68.9936589.4125273.322195.172-53.4082.2137631.0831324.51728.693-14.78-76.6339100.9724934.982386.50410.4664.4340426.5430319.81757.313、结算结果0.067530.003990.99772 -0.002210.001840.99990xs=39795.450.99771ys=27476.46 r= -0.0675

12、3zs=7572.69-0.00412九.实现程序输入文件形式如下:c+源程序如下：#include #include #include #include #include using namespace std;const int n=6;void inverse (double cnn);templatevoid transpose (t1*mat1,t2*mat2,int a,int b); templatevoid multi(t1*mat1,t2 * mat2,t2 * result,int a,int b,int c);templatevoid input (t*mat,int a

13、,int b);templatevoid output(t*mat,char*s,int a,int b);int main()ofstream outfile;cout.precision(5);double x0=0.0, y0=0.0; double fk=0.15324;/内方位元素double m=39689; /估算比例尺double b45=0.0,r33,xg61,at68,ata66,atl61;input (b,4,5);/从文件中读取控制点的影像坐标和地面坐标，存入数组bdouble xs=0.0, ys=0.0, zs=0.0,q=0.0,w=0.0,k=0.0;dou

14、ble x,y,z,l81,a86;/确定未知数的出始值for(int i=0;i4;i+)xs=xs+bi2;ys=ys+bi3;zs=zs+bi4;xs=xs/4; ys=ys/4; zs=zs/4+m*fk;int f=0;do/迭代计算f+;/组成旋转矩阵r00=cos(q)*cos(k)-sin(q)*sin(w)*sin(k);r01=-cos(q)*sin(k)-sin(q)*sin(w)*cos(k);r02=-sin(q)*cos(w);r10=cos(w)*sin(k);r11=cos(w)*cos(k);r12=-sin(w);r20=sin(q)*cos(k)+cos(

15、q)*sin(w)*sin(k);r21=-sin(q)*sin(k)+cos(q)*sin(w)*cos(k);r22=cos(q)*cos(w);/计算系数阵和常数项for(int i=0,k=0,j=0;i=6.0/206265.0|xg40=6.0/206265.0|xg50=6.0/2062 65.0);cout”迭代次数为：fendl;/精度评定double axg81,v81,vt18,vtv11,m0,d66;multi(a,xg ,axg,8,6,1);for( i=0;i8;i+)计算改正数vi0=axgi0-li0;transpose (/vt,1,8);multi(vt

16、,v,vtv,1,8,1);m0=vtv00/2;for(i=0;i6;i+)for(int j=0;j6;j+) dij=m0*ataij;/屏幕输出误差方程系数阵、常数项、改正数output(a,误差方程系数阵a为：,8,6);output(l,常数项 l 为：,8,1);output(xg,改正数为：,6,1);outfile.open(aim.txt,ios:app);/打开并添加 aim.txt 文件outfile.precision(10);/以文件的形式输出像片外方位元素、旋转矩阵、方差阵outfile 一、像片的外方位元素为：endlendl;outfilesetw(10)xs

17、=xssetw(10)ys=yssetw(10)zs=zsendl;outfilesetw(20)航向倾角为：qsetw(10)”旁向倾角为：wsetw(10)像片旋角为：kendl;outfileendl二、旋转矩阵 r 为:endlendl;for( i=0;i3;i+)for(int j=0;j3;j+)outfilesetw(25)rijsetw(25);outfileendl;outfileendl;outfilesetw(0)三、精度评定结果为：endl;outfile.precision(5);for(i=0;i6;i+)for(int j=0;j6;j+)outfilesetw

18、(14)dijsetw(14);outfileendl;outfile.close();return 0;templatevoid transpose(t1*mat1,t2*mat2,int a,int b) int i,j;for(i=0;ib;i+)for(j=0;ja;j+)mat2皿产mat1皿;return;templatevoid multi(t1*mat1,t2 * mat2,t2 * result,int a,int b,int c) int i,j,k;for(i=0;ia;i+)for(j=0;jc;j+)resultij=0;for(k=0;kb;k+)resultij+

19、=mat1ik*mat2kj;return;)template void input (t*mat,int a,int b) ifstream infile;infile.open(控制点坐标.txt);while(!infile.eof()for (int i=0;ia;i+)for(int j=0;jmatij;)infile.close();return;)templatevoid output(t*mat,char*s,int a,int b) coutsetw(15)sendl;for(int i=0;ia;i+)for(int j=0;jb;j+)coutsetw(13)matij

20、;coutendl;) return;)void inverse(double cnn) int i,j,h,k;double p;double qn12;for(i=0;in;i+)/构造高斯矩阵for(j=0;jn;j+)qij=cij;for(i=0;in;i+)for(j=n;j12;j+)if(i+6=j)qij=1;elseqij=0;for(h=k=0;kn-1;k+,h+)/消去对角线以下的数据for(i=k+1;in;i+)if(qih=0)continue;p=qkh/qih;for(j=0;j0;k-,h-) 消去对角线以上的数据 for(i=k-1;i=0;i-)if(

21、qih=0)continue;p=qkh/qih;for(j=0;j12;j+)q皿*=p;q皿产qk皿for(i=0;in;i+)/将对角线上数据化为1 p=1.o/qii；for(j=0;j12;j+)qij*=p;for(i=0;in;i+) / 提取逆矩阵for(j=0;jn;j+)cij=qij+6；程序的结果输出如下：（包括文本输出结果和荧屏输出中间数据） cju serswa ngweideilrtopsi：.副是程，生林a n haufa n gjiaohui.exe-2.8325e-0051.8876e-00fc1.6113e_005-0.203$8-0.025079-0.06899-1.8705e-00fc-2.841te-0051.28e-005-3,05096-0.181260.0b615-2.2250-005-q.18935q.04b1870.bh221-1.5562e-00fc-2.22e-00s-1.20046-0060.01s369-0.198830.0534-2.9484e-0051.993e-00b2.9731e-006-5.154970.00307-0.07663-1.937ba-qbb-2,?s27e-005lt729w-。瓯-b.02b278-q.19h620.b147