《单向后方交会实验报告C++.doc》由会员分享，可在线阅读，更多相关《单向后方交会实验报告C++.doc（13页珍藏版）》请在三一文库上搜索。
1、实验报告 班级：测绘一班 学号: 日期：2016.5.5目录一、计算原理3二、算法流程4三、源程序5四、计算结果13五、结果分析13六、心得体会13一、计算原理 已知条件摄影机主距f=153.24mm，x0=0.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.78-76.6339100.9724934.982386.50410.
2、4664.4340426.5430319.81757.31以单像空间后方交会方法，求解该像片的外方位元素。二、算法流程（1）获取已知数据。从航摄资料中差取平均航高与摄影机主距；获取控制点的地面测量坐标并转换为地面摄影坐标。（2）量测控制点的像点坐标并作系统误差改正。（3）确定未知数的初始值。在竖直摄影且地面控制点大体对称分布的情况下，按如下方法确定初始值，即，（4） 用三个角元素的初始值按下式，计算各个方向余弦值，组成旋转矩阵R （5） 逐点计算像点坐标的近似值。利用未知数的近似值和控制点的地面坐标；带入共线方程式，逐点近似像点坐标的近似值（x）、（y）。（6） 逐点计算误差方程式的系数和常数
3、项，组成误差方程式。（7） 计算法方程的系数矩阵和常数项，组成法方程式。 （8） 解法方程，求得外方位元素的改正数。（9） 用前次迭代取得的近似值，加本次迭代的改正数，计算外方位元素的新值。（10） 将求得的外方位元素改正数与规定的限差比较，若小于限差，则迭代结束。负责用新的近似值重复（4）-（9），直到满足要求为止。 用共线方程进行空间后方交会的程序框如图所示。输入原始数据像点坐标计算，系统误差正确定外方位因素初始值组成旋转矩阵R逐点组成误差方程式并法化所有像点完否是解法方程，求外方位元素改正数计算改正后的外方位元素外方位元素改正数是否小于限差是输出计算成果，计算并结束结束并显示错误信息迭代
4、次数小于n 是否 否2 源程序/2014113214 徐福辉#include#include#includeusing 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);int main()double x42 = -0.08616, -0.06897, -0.05341, 0.0822
5、3, -0.01479, -0.07661, 0.01045, 0.06445 ;double X43 = 36589.41, 25273.32, 2195.17, 37631.08, 31324.51, 728.69, 39100.97, 24934.98, 2386.50, 40426.54, 30319.81, 757.31 ;int i, j, m = 1;/n为迭代次数double X06 = 0 ;/设定未知数(XS,YS,ZS, ，)初始值double f = 0.15324;/摄影机主距f=153.24mmdouble a = 1 / 40000.0;/像片比例尺为1:4000
6、0double R33 = 0 ;/初始化旋转矩阵R double dayue_x8 = 0 ;/用于存放像点估计值double A86 = 0 ;/系数阵double AT68 = 0 ;/A的转置矩阵double L8 = 0 ;/存放常数项const double pi = 3.1415926535897932;double Asum66 = 0 ;double jieguo26 = 0 ;double jieguo168 = 0 ;double sumXYZ3 = 0 ;cout.precision(5);cout 已知像点坐标为：n;for (i = 0; i4; i+)for (j
7、 = 0; j2; j+)cout fixed;if (j = 0)cout x i + 1 = setw(10) xij ;elsecout y i + 1 = setw(10) xij endl;cout 已知地面四个点的坐标为：n;for (i = 0; i4; i+)for (j = 0; j3; j+)if (j = 0)cout X i + 1 = Xij ;elseif (j = 1)cout Y i + 1 = Xij ;elsecout Z i + 1; cout = ; cout Xij endl;cout endl;for (j = 0; j3; j+)for (i =
8、0; i4; i+)sumXYZj += Xij;for (i = 0; i2; i+)X0i = sumXYZi / 4;/X0,Y0初始化X0i = 1 / a*f + sumXYZ2 / 4.0;/Z0初始化doR00 = cos(X03)*cos(X05) - sin(X03)*sin(X04)*sin(X05);R01 = -cos(X03)*sin(X05) - sin(X03)*sin(X04)*cos(X05);R02 = -sin(X03)*cos(X04);R10 = cos(X04)*sin(X05);R11 = cos(X04)*cos(X05);R12 = -sin(
9、X04);R20 = sin(X03)*cos(X05) + cos(X03)*sin(X04)*sin(X05);R21 = -sin(X03)*sin(X05) + cos(X03)*sin(X04)*cos(X05);R22 = cos(X03)*cos(X04);/第一个像点的估计值，其坐标位于X00,X01,X02dayue_x0 = -f*(R00 * (X00 - X00) + R10 * (X01 - X01) + R20 * (X02 - X02) / (R02 * (X00 - X00) + R12 * (X01 - X01) + R22 * (X02 - X02);day
10、ue_x1 = -f*(R01 * (X00 - X00) + R11 * (X01 - X01) + R21 * (X02 - X02) / (R02 * (X00 - X00) + R12 * (X01 - X01) + R22 * (X02 - X02);/第二个像点的估计值,其坐标位于X10,X11,X12dayue_x2 = -f*(R00 * (X10 - X00) + R10 * (X11 - X01) + R20 * (X12 - X02) / (R02 * (X10 - X00) + R12 * (X11 - X01) + R22 * (X12 - X02);dayue_x3
11、 = -f*(R01 * (X10 - X00) + R11 * (X11 - X01) + R21 * (X12 - X02) / (R02 * (X10 - X00) + R12 * (X11 - X01) + R22 * (X12 - X02);/第三个像点的估计值,其坐标位于X20,X21,X22dayue_x4 = -f*(R00 * (X20 - X00) + R10 * (X21 - X01) + R20 * (X22 - X02) / (R02 * (X20 - X00) + R12 * (X21 - X01) + R22 * (X22 - X02);dayue_x5 = -f
12、*(R01 * (X20 - X00) + R11 * (X21 - X01) + R21 * (X22 - X02) / (R02 * (X20 - X00) + R12 * (X21 - X01) + R22 * (X22 - X02);/第四个像点的估计值,其坐标位于X30,X31,X32dayue_x6 = -f*(R00 * (X30 - X00) + R10 * (X31 - X01) + R20 * (X32 - X02) / (R02 * (X30 - X00) + R12 * (X31 - X01) + R22 * (X32 - X02);dayue_x7 = -f*(R01
13、 * (X30 - X00) + R11 * (X31 - X01) + R21 * (X32 - X02) / (R02 * (X30 - X00) + R12 * (X31 - X01) + R22 * (X32 - X02);for (i = 0; i4; i+)/第i个像点估计值放在dayue_x2*(i-1)A2 * i0 = (R00 * f + R02 * dayue_x2 * i) / (R02 * (Xi0 - X00) + R12 * (Xi1 - X01) + R22 * (Xi2 - X02);A2 * i1 = (R10 * f + R12 * dayue_x2 *
14、i) / (R02 * (Xi0 - X00) + R12 * (Xi1 - X01) + R22 * (Xi2 - X02);A2 * i2 = (R20 * f + R22 * dayue_x2 * i) / (R02 * (Xi0 - X00) + R12 * (Xi1 - X01) + R22 * (Xi2 - X02);A2 * i + 10 = (R01 * f + R02 * dayue_x2 * i + 1) / (R02 * (Xi0 - X00) + R12 * (Xi1 - X01) + R22 * (Xi2 - X02);A2 * i + 11 = (R11 * f +
15、 R12 * dayue_x2 * i + 1) / (R02 * (Xi0 - X00) + R12 * (Xi1 - X01) + R22 * (Xi2 - X02);A2 * i + 12 = (R21 * f + R22 * dayue_x2 * i + 1) / (R02 * (Xi0 - X00) + R12 * (Xi1 - X01) + R22 * (Xi2 - X02);A2 * i3 = dayue_x2 * i + 1 * sin(X04) - (dayue_x2 * i / f*(dayue_x2 * i * cos(X05) - dayue_x2 * i + 1 *
16、sin(X05) + f*cos(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 + 13 = -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);A2
17、 * i + 14 = -1 * f*cos(X05) - dayue_x2 * i + 1 / f*(dayue_x2 * i * sin(X05) + dayue_x2 * i + 1 * cos(X05);A2 * i + 15 = -dayue_x2 * i;/初始化常数项for (i = 0; i4; i+)L2 * i = xi0 - dayue_x2 * i;L2 * i + 1 = xi1 - dayue_x2 * i + 1;/A的转置矩阵for (i = 0; i8; i+)for (j = 0; j6; j+)ATji = Aij;/实现A与AT相乘int k = 0;f
18、or (i = 0; i6; i+)for (j = 0; j6; j+)Asumij = 0;for (i = 0; i6; i+)for (k = 0; k6; k+)for (j = 0; j8; j+)Asumik += ATij * Ajk;/得到AT*A的逆矩阵存放在inverseAsum66中inverse(Asum);/实现矩阵Asum66与AT68的相乘,结果存放在result168中for (i = 0; i6; i+)for (j = 0; j8; j+)jieguo1ij = 0;for (i = 0; i6; i+)for (k = 0; k8; k+)for (j
19、= 0; j6; j+)jieguo1ik += Asumij * ATjk;/实现result168与l8的相乘，得到结果放在result26中；for (i = 0; i6; i+)jieguo2i = 0;for (i = 0; i6; i+)for (j = 0; j8; j+)jieguo2i += jieguo1ij * Lj;for (i = 0; i6; i+)X0i = X0i + jieguo2i;ofstream f7(d:A.txt);f7 std:fixed;cout 进行第 m 次迭代带得到Xs,Ys,Zs, ，改正数结果为：n;for (i = 0; i6; i+
20、)cout setw(12) jieguo2i;f7 setw(12) jieguo2i;cout endl 6 | abs(jieguo24 * 206265.0)6 | abs(jieguo25 * 206265.0)6);cout n满足条件的结果为n;cout setw(12) Xs setw(12) Ys setw(12) Zs setw(12) setw(12) setw(12) endl;ofstream f7(d:A.txt);f7 std:fixed;cout.precision(4);for (i = 0; i6; i+)cout setw(12) X0i; f7 setw
21、(16) X0i;f7.close();double XG61;for (i = 0; i6; i+)XGi0 = jieguo2i;double AXG81, V81, VT18, VTV11, m0, D66;multi(A, XG, AXG, 8, 6, 1);for (i = 0; i8; i+) /计算改正数Vi0 = AXGi0 - 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
22、 = 0; i6; i+)for (int j = 0; j6; j+)Dij = m0*Asumij;cout setw(10) Dij;f6 setw(10) Dij;cout endl;f6 endl;cout所得中误差为endl;for (i = 0; i6; i+)cout sqrt(Dii) endl;f6.close(); getchar();return 0;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
23、 (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 (qih = 0)continue;p = qkh / qih;for (j = 0; j12; j+)qij *= p;qij -= qkj;for (i = 0; in;
24、i+)p = 1.0 / qii;for (j = 0; j12; j+)qij *= p;for (i = 0; in; i+) for (j = 0; jn; j+)cij = qij + 6;templatevoid transpose(T1*mat1, T2*mat2, int a, int b)int i, j;for (i = 0; ib; i+)for (j = 0; ja; j+)mat2ji = mat1ij;return;templatevoid multi(T1*mat1, T2 * mat2, T2 * result, int a, int b, int c)int i
25、, j, k;for (i = 0; ia; i+)for (j = 0; jc; j+)resultij = 0;for (k = 0; kb; k+)resultij += mat1ik * mat2kj;return; 程序截图3、 计算结果XsYsZs39795.855427475.86947572.5475-0.00400.00210.06764、 结果分析由计算结果可知:摄影中心在地面摄影测量坐标系中的坐标为(39795.452258.27476.462385. 7572.685988)m航向倾角为-0.003987弧度，旁向倾角为0.002114弧度，像片旋角为0.067578弧度。参数XsYsZs中误差1.26201.42390.55650.00020.00020.0001六、心得体会 关于本次实验，我的体会实在很多，从哪说起呢，我的编程能力原本是比较差的，说实话刚拿到这次课程设计的时候我觉得自己没法完成，但还好我并没有放弃，自己看了好多遍书，首先理解了实验原理，然后又复习了一下c+，当然在实验过程中也问了不少同学，才完成了这次实验，不过过程虽然艰难，但完成时的那种满足感真的是难以用语言表达，我觉得通过这次实验，我学到的或许比我想象的还多。