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

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<iostream>#include<fstream>#include<iomanip>using namespace std;const int n = 6;void inverse(double cnn);template<typename T1, typename T2>void transpose(T1*mat1, T2*mat2, int a, int b);template<typename T1, typename T2>void multi(T

5、1*mat1, 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 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为迭代次

6、数double X06 = 0 ;/设定未知数(XS,YS,ZS, ，)初始值double f = 0.15324;/摄影机主距f=153.24mmdouble a = 1 / 40000.0;/像片比例尺为1:40000double 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 jieg

7、uo26 = 0 ;double jieguo168 = 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) << xij << ""

8、elsecout << "y" << i + 1 << "= " << setw(10) << xij << endl;cout << "已知地面四个点的坐标为：n"for (i = 0; i<4; i+)for (j = 0; j<3; j+)if (j = 0)cout << "X" << i + 1 << "= " << Xij <<

9、""elseif (j = 1)cout << "Y" << i + 1 << "= " << Xij << ""elsecout << "Z" << i + 1; cout << "= " cout << Xij << endl;cout << endl;for (j = 0; j<3; j+)for (i = 0; i<4; i+

10、)sumXYZj += Xij;for (i = 0; i<2; 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(X04);

11、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);dayue_x1

12、 = -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 = -f

13、*(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*(R01

14、 * (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 * (X

15、30 - X00) + R11 * (X31 - X01) + R21 * (X32 - X02) / (R02 * (X30 - X00) + R12 * (X31 - X01) + R22 * (X32 - X02);for (i = 0; i<4; 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 * i

16、) / (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 +

17、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 * s

18、in(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

19、* 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; i<4; i+)L2 * i = xi0 - dayue_x2 * i;L2 * i + 1 = xi1 - dayue_x2 * i + 1;/A的转置矩阵for (i = 0; i<8; i+)for (j = 0; j<6; j+)ATji = Aij;/实现A与AT相乘

20、int k = 0;for (i = 0; i<6; i+)for (j = 0; j<6; j+)Asumij = 0;for (i = 0; i<6; i+)for (k = 0; k<6; k+)for (j = 0; j<8; j+)Asumik += ATij * Ajk;/得到AT*A的逆矩阵存放在inverseAsum66中inverse(Asum);/实现矩阵Asum66与AT68的相乘,结果存放在result168中for (i = 0; i<6; i+)for (j = 0; j<8; j+)jieguo1ij = 0;for (i

21、 = 0; i<6; i+)for (k = 0; k<8; k+)for (j = 0; j<6; j+)jieguo1ik += Asumij * ATjk;/实现result168与l8的相乘，得到结果放在result26中；for (i = 0; i<6; i+)jieguo2i = 0;for (i = 0; i<6; i+)for (j = 0; j<8; j+)jieguo2i += jieguo1ij * Lj;for (i = 0; i<6; i+)X0i = X0i + jieguo2i;ofstream f7("d:A.

22、txt");f7 << std:fixed;cout << "进行第" << m<< "次迭代带得到Xs,Ys,Zs, ，改正数结果为：n"for (i = 0; i<6; i+)cout << setw(12) << jieguo2i;f7 << setw(12) << jieguo2i;cout << endl << endl;f7.close();getchar();m+=1; while (abs(jieguo2

23、3 * 206265.0)>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) << "" <<

24、setw(12) << "" << setw(12) << "" << endl;ofstream f7("d:A.txt");f7 << std:fixed;cout.precision(4);for (i = 0; i<6; i+)cout << setw(12) << X0i; f7 << setw(16) << X0i;f7.close();double XG61;for (i = 0; i<6; i+)XG

25、i0 = jieguo2i;double AXG81, V81, VT18, VTV11, m0, D66;multi(A, XG, AXG, 8, 6, 1);for (i = 0; i<8; 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

26、 = 0; i<6; i+)for (int j = 0; j<6; j+)Dij = m0*Asumij;cout << setw(10) << Dij;f6 << setw(10) << Dij;cout << endl;f6 << endl;cout<<"所得中误差为"<<endl;for (i = 0; i<6; i+)cout << sqrt(Dii) << endl;f6.close(); getchar();return 0

27、;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 (j = n; j<12; j+)if (i + 6 = j)qij = 1;elseqij = 0;for (h = k = 0; k<n - 1; k+, h+)for (i = k + 1; i<n; i+)if (qih = 0)continue;p = qkh / qih;for (j = 0; j<12; j+)qij *= p;qij -= qkj;for (h = k = n - 1; k>0; k-, h-) for (i = k - 1; i >= 0; i-)if (qih = 0)continue;p = qkh / qih;for (j = 0; j<12; j+)qij *= p;qij -= qkj;for (i = 0; i<n; i+)p = 1.