單向后方交會(huì)實(shí)驗(yàn)報(bào)告C++

1、實(shí)驗(yàn)報(bào)告 班級(jí)：測繪一班 學(xué)號(hào): 日期：2016.5.5目錄一、計(jì)算原理3二、算法流程4三、源程序5四、計(jì)算結(jié)果13五、結(jié)果分析13六、心得體會(huì)13一、計(jì)算原理 已知條件攝影機(jī)主距f=153.24mm，x0=0.01mm，y0=0.02mm, 像片比例尺為1:40000，有四對(duì)點(diǎn)的像點(diǎn)坐標(biāo)與相應(yīng)的地面坐標(biāo)如下表。 點(diǎn)號(hào)像點(diǎn)坐標(biāo)地面坐標(biāo)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以單像空間后方交會(huì)方法，求解該像片的外方位元素。二、算法流程（1）獲取已知數(shù)據(jù)。從航攝資料中差取平均航高與攝影機(jī)主距；獲取控制點(diǎn)的地面測量坐標(biāo)并轉(zhuǎn)換為地面攝影坐標(biāo)。（2）量測控制點(diǎn)的像點(diǎn)坐標(biāo)并作系統(tǒng)誤差改正。（3）確定未知數(shù)的初始值。在豎直攝影且地面控制點(diǎn)大體對(duì)稱分布的情況下，按如下方法確定初始值，即，（4） 用三個(gè)角元素的初始值按下式，計(jì)算各個(gè)方向余弦值，組成旋轉(zhuǎn)矩陣R （5） 逐點(diǎn)計(jì)算像點(diǎn)坐標(biāo)的近似值。利用未知數(shù)的近似值和控制點(diǎn)的地面坐標(biāo)；帶入共線方程式，逐點(diǎn)近似像點(diǎn)坐標(biāo)的近似值（x）、（y）。（6） 逐點(diǎn)計(jì)算誤差方程式的系數(shù)和常數(shù)

3、項(xiàng)，組成誤差方程式。（7） 計(jì)算法方程的系數(shù)矩陣和常數(shù)項(xiàng)，組成法方程式。 （8） 解法方程，求得外方位元素的改正數(shù)。（9） 用前次迭代取得的近似值，加本次迭代的改正數(shù)，計(jì)算外方位元素的新值。（10） 將求得的外方位元素改正數(shù)與規(guī)定的限差比較，若小于限差，則迭代結(jié)束。負(fù)責(zé)用新的近似值重復(fù)（4）-（9），直到滿足要求為止。 用共線方程進(jìn)行空間后方交會(huì)的程序框如圖所示。輸入原始數(shù)據(jù)像點(diǎn)坐標(biāo)計(jì)算，系統(tǒng)誤差正確定外方位因素初始值組成旋轉(zhuǎn)矩陣R逐點(diǎn)組成誤差方程式并法化所有像點(diǎn)完否是解法方程，求外方位元素改正數(shù)計(jì)算改正后的外方位元素外方位元素改正數(shù)是否小于限差是輸出計(jì)算成果，計(jì)算并結(jié)束結(jié)束并顯示錯(cuò)誤信息迭代

4、次數(shù)小于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、數(shù)double X06 = 0 ;/設(shè)定未知數(shù)(XS,YS,ZS, ，)初始值double f = 0.15324;/攝影機(jī)主距f=153.24mmdouble a = 1 / 40000.0;/像片比例尺為1:40000double R33 = 0 ;/初始化旋轉(zhuǎn)矩陣R double dayue_x8 = 0 ;/用于存放像點(diǎn)估計(jì)值double A86 = 0 ;/系數(shù)陣double AT68 = 0 ;/A的轉(zhuǎn)置矩陣double L8 = 0 ;/存放常數(shù)項(xiàng)const double pi = 3.1415926535897932;double Asum66 = 0 ;double jieg

7、uo26 = 0 ;double jieguo168 = 0 ;double sumXYZ3 = 0 ;cout.precision(5);cout << "已知像點(diǎn)坐標(biāo)為：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 << "已知地面四個(gè)點(diǎn)的坐標(biāo)為：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);/第一個(gè)像點(diǎn)的估計(jì)值，其坐標(biāo)位于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);/第二個(gè)像點(diǎn)的估計(jì)值,其坐標(biāo)位于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);/第三個(gè)像點(diǎn)的估計(jì)值,其坐標(biāo)位于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);/第四個(gè)像點(diǎn)的估計(jì)值,其坐標(biāo)位于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個(gè)像點(diǎn)估計(jì)值放在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;/初始化常數(shù)項(xiàng)for (i = 0; i<4; i+)L2 * i = xi0 - dayue_x2 * i;L2 * i + 1 = xi1 - dayue_x2 * i + 1;/A的轉(zhuǎn)置矩陣for (i = 0; i<8; i+)for (j = 0; j<6; j+)ATji = Aij;/實(shí)現(xiàn)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);/實(shí)現(xiàn)矩陣Asum66與AT68的相乘,結(jié)果存放在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;/實(shí)現(xiàn)result168與l8的相乘，得到結(jié)果放在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 << "進(jìn)行第" << m<< "次迭代帶得到Xs,Ys,Zs, ，改正數(shù)結(jié)果為：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滿足條件的結(jié)果為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+) /計(jì)算改正數(shù)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<<"評(píng)定完精度的結(jié)果是"<<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.