单向后方交会实验报告C++

1、southwest jiaotong university班级，测绘一班 学号：日期：2016.5.5目录一、计算原理3二、算法流程4三、源程序5四、计算结果13五、结果分析13六、心得体会13一、计算原理已知条件摄影机主距f=153. 24mm, xo=o. 01mm, y0=0. 02mm,像片比例尺为算法流程1:40000,有四对点的像点坐标与相应的地面坐标如卜表。勺像点坐标地面坐标x (mm)y (mm)x(m)y(m)z(m)1-86. 15-68. 9936589.4125273. 322195. 172-53. 4082.2137631.0831324.51728. 693-14

2、. 78-76. 6339100.9724934. 982386.50410. 4664. 4340426.5430319.81757. 31以单像空间后方交会方法，求解该像片的外方位元索。(1) 获取已知数据。从航摄资料巾差取平均航高与摄影机主距；获取控制点的 地面测量坐.标并转换为地面摄影坐.标。(2) 量测控制点的像点坐标并作系统误差改正。(3) 确定未知数的初始值。在竖直摄影且地面控制点大体对称分布的情况卜*, 按如卜*方法确定初姑值，即n，z5° = inf + nn=0(4) 用三个角元素的初始值按下式，计算各个方向余弦值，组成旋转矩阵ral = cos (p cos /

3、c - sin (p sin co sin k a2= -cospsin zc-sin sin69coszc a3 = - sin cos 69/?, = cos 69 sin zc b2 = cos69coszcb3 = - sin o)c, = sin cos k + cos sin ?sin k c2 =-sin(p sin k + cos sin 69 cos k c3 = cos 炉 cos 仞(5) 逐点计算像点坐标的近似值。利用未知数的近似值和控制点的地而坐标； 带入共线方程式，逐点近似像点坐标的近似值(x)、(y)。(6) 逐点计算误差方程式的系数和常数项，组成误差方程式。(7

4、) 计算法方程的系数矩阵和常数项组成法方程式。(8) 解法方程，求得外方位元素的改正数6/xs、e/ys、d(p、do)、dk。(9) 用前次迭代取得的近似值，加木次迭代的改正数，计算外方位元素的新值。+dxks,ysk =ysk' +dysk,zks = z；-' +dzks(pk =(pk-' +d(pk，0)k =o)ks do)k ,kk =kkx +dkk(10) 将求得的外方位元素改正数与规定的限差比较,若小于限差，则迭代结朿。 负责用新的近似值重复(4) - (9),直到满足要求为止。用共线方程进行空间后方交会的程序框如阁所示。解法方程，求外方位元素改正数

5、并显示错误信息否v计算改正后的外方位元素; 外方位元素改¥数是否小于限差7 是小/否二.源程序/2014113214 徐福辉#include<iostream>#include<fstream>#include<iomanip>using namespace std;const int n = 6;void inverse(double cn n);tempiate<typename tl, typename t2>void transpose(t1*mat1, t2氺!nat2, int a, int b); template<

6、typename tl, typename t2>void multi (tl*matl, t2*mat2, t2 * result, int a, int b，int c);int main ()double x42 = -0.08616, -0.06897，-0. 05341, 0.08223, -0.01479, -0.07661， 0.01045, 0. 06445 ;double x4 3 = 36589.41, 25273.32, 2195.17, 37631.08, 31324.51, 728.69, 39100.97, 24934.98, 2386.50, 40426.5

7、4, 30319.81, 757.31 ;int i, j, m = l;/n为迭代次数double x06 = 0 ;/设定未知数(xs, ys, zs, v, o , k)初始值 double f = 0. 15324;/摄影机主距 f=153. 24mm double a = 1 / 40000. 0;/像片比例尺为 1:40000 double r33 = 0 ;/初始化旋转矩阵r double dayue_x8 = 0 ;/用于存放像点估计值double a8 6 = 0 ;/系数阵double at6 8 = 0 ;/a 的转置矩阵double l8 = 0 ;/存放常数项cons

8、t double pi = 3. 1415926535897932;double asum66 = 0 ;double jieguo26 = 0 ;double jieguol68 = 0 ;double sumxyz3 = 0 ;cout. precision(5);cout « "已知像点坐标为：n"for (i = 0; i<4; i+)for (j = 0; j<2; j+十)cout « fixed;if (j = 0)cout << x i + 1 < = setw(10) « xijelsecout&

9、lt; y< i + 1 << = << setw(lo) « xij<cout « "已知地面四个点的坐标为:n;for (i = 0; i4; i+)for (j = 0; j3; j+)if (j = 0)cout « /zx/zi + 1 « =xij< 'elseif (j = 1)cout << "y" « i + 1 « /z= /z « xij «elsecout "z + 1; cout =end

10、l;cout << endl;for (j = 0; j3; j+)for (i = 0; i<4; i+)sumxyzej += xij;for (i = 0; i<2; i+)x0i = suraxyzti / 4;/x0，y0 初始化x0i = 1 / a*f + sumxyz2 / 4.0;/z0 初始化 dor00 = cos(x03)*cos(x05) - sin(x03)*sin(x04)*sin(x05); r0l = _cos(x03)*sin(x05) - sin(x03)*sin(x04)*cos(x05);r02 = -sin(x03)*cos(

11、x04);rl0 = cos(x04)*sin(x05);rll = cos(x04)*cos(x05);rl2 = -sin(x04)；r20 = sin(x03)*cos(x05) + cos(x03)*sin(x04)*sin(xo5); r2l = -sin(x03)*sin(x05) + cos(xo3)*sin(xo4)*cos(xo5): r2 2 = cos(x03)*cos(x04);/第一个像点的估计值,其坐标位于x0 0, x0 1, x0 2dayue_x0 = -f*(roo * (x00 - x00) + rl0 * (x0l - x0l)+ r20 * (x0 2

12、 -x02) / (r02 * (x00 - x00) + rl2 * (x0l -x0l) + r22 * (x02 - x02);dayue_xl =-f*(ro 1 * (x00 - x00) + rl 1 * (x0l - x0l)+ r2l * (x0 2 - x02) / (r02 * (x00 - x00) + rl2 * (x0l -x0l) + r22 * (x02 - x02);/第二个像点的估计值,其坐标位于xl0,xll,xl2dayue_x2 = -f*(roo* (xl0- x00)+rl0 * (xll- x0l)+ r20 * (xl2 -x02) /(r02

13、*(xl0- x00) +r12 *(xll-x0l) + r22 * (xl2 - x02);dayue_x3 =-f*(ro 1* (xl0- x00)+r11 * (xll- x0l)+ r21 * (xl2 - x02) /(r02 *(xl0 x00) +r12 *(xll -x0l) + r22 * (xl2 - x02);/第三个像点的估计值,其坐标位于x2 0, x2 1, x2 2dayue_x4 =-f*(ro 0 * (x20 一 x00) +rl0 * (x2l - x0l)+ r20 * (x22 -x02) / (r0 2 * (x20 - x00) +r12 *

14、(x2 1 -x0l) + r22 * (x22 - x02);dayue_x5 = f*(ro 1 * (x20 - x00) + rl 1 * (x2l - x0l)+ r2l * (x2 2 -x02) / (r02 * (x2 0 - x00) +r12 * (x2l -x0l) + r22 * (x22 - x02);/第四个像点的估计值，其坐标位于x3 0，x3 1, x3 2dayue_x6 =-f*(ro 0* (x30- x00)+rl0 * (x3l- x0l)+ r20 * (x3 2-x02) /(r02 *(x30- x00) +r12 *(x3l-xol) + r2

15、2 * (x32 - x02);dayue_x7 =-f*(ro 1* (x30- x00)+r11 * (x3l- x0l)+ r21 * (x32-x02) /(r02 *(x30- x00) + rl 2 *(x3 1-x0l) + r22 * (x32 - x02); for (i = 0; i<4; i+)/第i个像点估计值放在dayue_x2*(i-l)a2*i0 = (r0 0*f+ r02 * dayue_x2 * i)/ (r02 *(xi0-x00)+rl2 * (xil -x0l) + r22 * (xi2- x02);a2*il = (rl0*f + rl2 *d

16、ayue_x2*i)/ (r02 *(xi0-x00)+rl2 * (xil -xo1) + r22 * (xi 2- x02);a2*i2 = (r2 0*f + r22 * dayue. x2* i)/ (r02 *(xi 0-x00)+rl2 * (xil -x0l) + r22 * (xi2- x02);a2 * i + 1 0 = (r0 1 * f + r0 2 * dayue_x2 * i + 1) / (r0 2* (xi0 -x00) +r12 * (xil -x0l) +r22 * (xi2 - x02);a2 * i + 11 = (rl 1 * f + rl 2 * d

17、ayue_x2 * i + 1) / (r0 2* (xi0 -xoo) +r12 * (xil -xol) +r22 * (xi 2 - x02);a2* i + 12 = (r2 1 * f + r2 2 * dayue_x2 * i + 1) / (r0 2* (xi0 -x00) +r12 * (xil -x0l) +r22 * (xi2 - x02);a2 * i3 = dayue_x2 * i + 1 * sin(x04) - (dayue_x2 * i / f*(dayue_x2 * i * cos(x05) - dayue_x2 * i + 1 * sin(x05) + f*c

18、os (x05)*cos(x04);a2 * i4 = -f氺sin(x05) - dayue_x2 * i / f*(dayue_x2 * i * sin(x05) + dayue_x2 * i + 1 * cos(x05);a2 * i5 = dayue_x2 * i + 1;a2 * i + 1 3 = -1 * dayue_x2 * i * sin(x04) - (dayue_x2 * i + 1 / f*(dayue x2 * i * cos(x05) - dayue_x2 * i + 1 * sin(x05)- f*sin(x05)*cos(x04);八2 * i + 1 4 =

19、-1 * f*cos(x05) - dayue_x2 * i + 1 / f*(dayue_x2* i * sin(x05) + dayue_x2 * i + 1 * cos(x05);a2 * i + 1 5 = -dayue x2 * i;/初始化常数项 for (i = 0; i<4; i+)l2 * i = xi0 - dayue_x2 * i;l2 * i + 1 = xil - davue_x2 * i + 1;/a的转®矩阵for (i = 0; i<8; i+)for (j = 0; j<6; j+)atji = aij;/实现a与at相乘 int

20、k = 0;for (i = 0; i<6; i+十) for (j = 0; j6; j+)asumij = 0;for (i = ; i<6; i+)for (k = 0; k<6; k+)for (j = 0; j<8; j+)asumik += atij * ajk;/得到at*a的逆矩陈存放在inverseasum6 6屮inverse(asum);/实现矩阵asum6 6与at6 8的相乘,结果存放在result 1 6 8中for(i =0; i<6;i+)for(j =0; j<8;j+)jieguolij=0for(i =0; i<6

21、;i+十)for(k =0; k<8;k+)for(j =0； j<6;j+)jieguolik += asumij * atjk;/实现result 1 6 8与18的相乘，得到结果放在result26中; for (i = 0; i6; i+)jieguo2i = 0; for (i = 0; i<6; i+) for (j = 0; j<8; j+)jieguo2i += jieguoli j * lj; for (i = 0; i<6; i+十)x0i = x0i + jieguo2i;ofstream f7("d:a. txt");f

22、7 << std:fixed;emit <进行笫< m次迭代带得到xs，ys，zs，v，co, k改正数结果为:for (i = 0; i<6; i+)cout « setw(12) « jieguo2i; f7 << setw(12) << jieguo2i;cout « cndl « endl;f7. closeo ;getcharo ;m+=l; while (abs(jieguo23 * 206265.0)>6 | abs(jieguo24 * 206265.0)>6 | abs(

23、jieguo25 * 206265. 0)>6);cout « "什满足条件的结果为n"cout « setw(12) << "xs"< setw(12) « ys" setw(12) « zs" < setw(12) << "p setw(12) << " g) "< setw(12) << " k " << endl;ofstream f7("d:a

24、. txt");f7 « std:fixed; cout. precision； for (i = 0; i6; i+)cout « setw(12) « x0i; f7 « setw(16) « x0i;f7. close();double xg61;for (i = 0; i<6; i+)xgi0 = jieguo2i;double axg81, v8l, vt18, vtv1 1, mo, d6 6;multi (a, xg, axg, 8, 6, 1);for (i = 0; i8; i+)/计算改正数vi0 = axg

25、i0 - li; transpose (v, vt, 1, 8); multi (vt, v, vtv, 1, 8, 1); m0 = vtv00 / 2; cout « endl;ofstream f6(d:what. txt)； cout"评记完精度的结果是"<<endl; for (i = 0; i<6; i+)for (int j = 0; j<6; j+)dij = m0*asumij; cout « setw(10) « dij; f6 « setw(10) << dij;cout &#

26、171; endl; f6 << endl;cout<"所得中误差为"<endl; for (i = 0; i<6; i+)cout « sqrt(dii) « endl; f6. close(); getcharo ; return 0;void inverse(double cnn)int i, j, h, k; double p; double qn12;for (i = 0; i<n; i+)for (j = 0; j<n; j+) qij = cij;for (i = 0; i<n; i+)for

27、 (j = n; j<12; j+)if (i + 6 = j)= 1；elseciij = 0;for (h = k = 0; k<n - 1; k十+，h十十) for (i = k + 1; in; i+)if (qi h = 0)continue;p = qkh / qih; for (j = 0; j12; j+)qi j p; qij -=for (h = k = n - 1; k0; k，h) for (i = k - 1; i= 0; i)if (qi h = 0)continue;p = qk h / qi h; for (j = 0; j<12; j+)q

28、i j *= p; qij -= qkj;for (i = 0; i<n; i+)p = 1.0 / qi i; for (j = 0; j<12; j+)qi j p;for (i = 0; i<n; i+十) for (j = 0; j<n; j+)ci j = qi j + 6;tempiate<typename tl, typename t2>void transpose(tlmatl, t2氺mat2, int a, int b)int i, j;for (i = 0; i<b; i+) for (j = 0; j<a; j+)mat2

29、j i = matli j; return;template<typename tl, typename t2>void multi(t1氺matl, t2 * mat2, t2 * result, int a, int b, int c)int i, j, k;for (i = 0; ia; i+)for (j = ; j<c; j+)resultfij = 0; for (k = 0; k<b; k+)resultij += matlik * mat2kj;return;程序截街专范®)sot2014113214 #mclude<iost ream&

30、gt; »mclude<fstr©am> #include<iojnanip> using namespace std; const int n = 6; void inverse(double cn template<typenanie tl, typ tcmp 1 ate<typename tl, typ int maino double x4 2 = -0 double x 4 3 = 36 int i, j, n = 1:/n为 double x06 = 0 ： double f = 0.15324;/ double a = 1

31、/ 40000. double r33 = 0 double dayue_x8 = double a86】= 0 double at68 = 0 double l8 = 0 ./, const double pi = 3. 1 double asun66 = dnuhlppir卜像占世标为.1 1 承rj'/ j xl =-0.08616yl=-0.0689?x2-0.05341y2=0.08223x3-0.01479ir3-0.076610.01045y4-0.06445已知地面四个点的坐标为：xl =36589.41000v1=25273.3200021=2195.17000x2

32、=37631.08000?2=31324.51000z2=728.69000x3 =39100.970003=24934.98000z3=2386.50000x4-40426.54000v430319.810b0z4-757.3100033 c:windowssystem32cmd.exe行第i次迭代帯得到iek改正数结果为：1406.38334-398.38130-53.00183-0.004420.00194行第2次迭代带得到xsjs.zs，冲，k改正数结果为： -47.80937-89.29912-20.905360.00b440.0b015m亍第3次迭代带得到xs.vs.zs，id,

33、«, k改正数结果为：0.281420.39487-0.06283-0.00002-0.00001-0.06757-0.000030.000的yjj.cpp x(全围)si c:windowssystem32cmd.exeb/2014113214 摘福辉#include<iostrcaa>#include<fstreanl>进行第3次迭代带得到xs.vs，zs，co, k改正数结果为:s#include<iomanip>0.281420.39487-0.06283-0.00002一0.000010.60000using nanespace std;const mt n = 6;void inverse(double c