摄影测量单像空间后方交会
10_GIS_2班
已知摄影机主距f=153.24mm，四对点的像点坐标与相应的地面坐标列入下表： 点号 像点坐标 地面坐标 x(mm) y(mm) x(m) y(m) z(m) 1 -86.15 -68.99 36589.41 25273.32 2195.17 2 -53.40 82.21 37631.08 31324.51 728.69 3 -14.78 -76.63 39100.97 24934.98 2386.50 4 10.46 64.43 40426.54 30319.81 757.31 求六个方位元素及其中误差 运用C++编写程序，具体步骤如下图：
首先将原始数据存储到.txt文本下，供后期编程过程中调用
程序代码如下：
// 空间后方交会.cpp : Defines the entry point for the console application. //
#include #include \"mtxfunction.h\" //矩阵求逆运算、转置、相乘 using namespace std; /*************************************** 声明结构类型GCP，用来保存控制点坐标 ***************************************/ struct GCP { double x; //影像控制点x坐标 double y; //影像控制点y坐标 矿业工程学院_10级GIS_215 - 1 - 摄影测量—空间后方交会 double Xt; //地面控制点Xt坐标 double Yt; //地面控制点Yt坐标 double Zt; //地面控制点Zt坐标 }; /***************************************/ // 声明读取文本数据函数 int Inputdata(int &GCPN, GCP *data, double &m, double &f, double &x0, double &y0); /***************************************/ // int main() { system(\"title 单像空间后方交会\"); //声明变量 int GCPN=100; //控制点对数GCPN double m =50000, f = 153.24, x0 = 0, y0 = 0; //比例尺m，主距f，外方位元素x0,y0 GCP data[4]; //控制点坐标矩阵 double Xs = 0, Ys = 0, Zs = 0, fai = 0, omiga = 0, k = 0; //内个外方位元素 double mtxX[6]; //X矩阵 X=[ΔXs,ΔYs,ΔZs,Δψ,Δω,Δκ] mtxX[3] = 1; double mtxR[9]; //旋转矩阵R double mtxresult1[6][6], mtxresult2[6]; //1保存ATA及其逆矩阵 和2保存ATL，相乘得到X结果 int i, n = 0; double sum = 0, m0; /***************************************/ // 读取数据 gotokaishi: ; //goto标记 if(Inputdata(GCPN, data, m, f, x0, y0)) { cout << \"读取“原始数据.txt”出现错误！\\n\"; cout << \"请将“原始数据.txt”置于程序所在目录中，并检查数据格式。\"; cout << \"\\n数据格式如下所示\\n(第一行为控制点个数和m，第二行为f,x0,y0，其余行为控制点坐标x,y,Xt,Yt,Zt)\\n\"; cout << \"4\0.0\\n\"; cout << \"153.24\0.0\0.0\\n\"; cout << \"-86.15\-68.99\36589.41\25273.32\2195.17\\n\"; cout << \"-53.40\ 82.21\37631.08\31324.51\728.69\\n\"; cout << \"-14.78\-76.63\39100.97\24934.98\2386.50\\n\"; cout << \" 10.46\ 64.43\40426.54\30319.81\757.31\\n\"; cout << \"是否重新读取文件？(Y or N)：\"; char ch; 矿业工程学院_10级GIS_215 - 2 - 摄影测量—空间后方交会 cin >> ch; if('Y' == ch || 'y' == ch) { system(\"cls\"); goto gotokaishi; } cout << \"程序即将退出，谢谢使用！\\n\" ; system(\"PAUSE\"); return 0; //读取数据不正确，返回1 } //计算Xs0,Ys0,Zs0 for(int i = 0; i < GCPN; ++i) { Xs += data[i].Xt / 4; Ys += data[i].Yt / 4; } Zs = f * m; array double mtxA[2 * 50][6]; //系数矩阵A 4个a11~a16，a21~a26 double mtxAT[6][2 * 50]; //A的转置 double mtxL[2 * 50]; //L矩阵 L=[lx,ly] cout << \"主距f=\" << f << \"比例尺m=\" << m << endl; cout << \"像点坐标\\\\" << \"地面点坐标\\n\" << \"x\\\" << \"y\\\" << \"Xt\\\" << \"Yt\\\" << \"Zt\\n\"; for(int i = 0; i < GCPN; ++i) { cout << data[i].x << \" \" << data[i].y << \" \" << data[i].Xt << \" \" \" \" << data[i].Yt << \" \" << data[i].Zt << endl; } /***************************************/ // 迭代计算 /***************************************/ while((mtxX[3] > 0.00001 || mtxX[4] > 0.00001 || mtxX[5] > 0.00001) && n < 100) { 矿业工程学院_10级GIS_215 - 3 - 摄影测量—空间后方交会 /*******每次迭代钱将ATA和ATL清空*********/ for(int i = 0; i < 6; ++i) { for(int j = 0; j < 6; ++j) { mtxresult1[i][j] = 0; } mtxresult2[i] = 0; } ++n;//统计迭代次数 cout << \"\\n******************************************\" \"***********************************\"; printf(\"\\n第%d次迭代\\n\ cout << \"Xs=\" << Xs << \"\Ys=\" << Ys << \"\Zs=\" << Zs << \"\\nψ=\" << fai << \"\ω=\" << omiga << \"\κ=\" << k << endl; ////////////****************旋转矩阵R************************/ mtxR[0] = cos(fai) * cos(k) - sin(fai) * sin(omiga) * sin(k); mtxR[1] = -cos(fai) * sin(k) - sin(fai) * sin(omiga) * cos(k); mtxR[2] = -sin(fai) * cos(omiga); mtxR[3] = cos(omiga) * sin(k); mtxR[4] = cos(omiga) * cos(k); mtxR[5] = -sin(omiga); mtxR[6] = sin(fai) * cos(k) + cos(fai) * sin(omiga) * sin(k); mtxR[7] = -sin(fai) * sin(k) + cos(fai) * sin(omiga) * cos(k); mtxR[8] = cos(fai) * cos(omiga); printf(\"-------------旋转矩阵R:-----------\\n\"); cout << mtxR[0] << '\' << mtxR[1] << '\' << mtxR[2] << endl; cout << mtxR[3] << '\' << mtxR[4] << '\' << mtxR[5] << endl; cout << mtxR[6] << '\' << mtxR[7] << '\' << mtxR[8] << endl << endl; ////////////**********************像点坐标计算值************************/ //cout<<\"-------------归算得到控制点坐标------------\\n\"; //cout<<\"Xtg\\"<<\" \" <<\"Ytg\\"<<\" \"<<\"Ztg\\\"< /******************归算控制点坐标************************/ Xtg[i] = mtxR[0] * (data[i].Xt - Xs) + mtxR[3] * (data[i].Yt - Ys) + mtxR[6] * 矿业工程学院_10级GIS_215 - 4 - 摄影测量—空间后方交会 (data[i].Zt - Zs); Ytg[i] = mtxR[1] * (data[i].Xt - Xs) + mtxR[4] * (data[i].Yt - Ys) + mtxR[7] * (data[i].Zt - Zs); Ztg[i] = mtxR[2] * (data[i].Xt - Xs) + mtxR[5] * (data[i].Yt - Ys) + mtxR[8] * (data[i].Zt - Zs); double x_x0, y_y0; //x-x0=-f(X/Z)见书38页（2-5-8） x_x0 = -f * Xtg[i] / Ztg[i]; y_y0 = -f * Ytg[i] / Ztg[i]; /**************计算L矩阵*********************/ mtxL[2 * i] = data[i].x - (x_x0); mtxL[2 * i + 1] = data[i].y - (y_y0); /*************计算A矩阵**********************/ mtxA[2 * i][0] = (mtxR[0] * f + mtxR[2] * (x_x0)) / Ztg[i]; mtxA[2 * i][1] = (mtxR[3] * f + mtxR[5] * (x_x0)) / Ztg[i]; mtxA[2 * i][2] = (mtxR[6] * f + mtxR[8] * (x_x0)) / Ztg[i]; mtxA[2 * i][3] = (y_y0) * sin(omiga) - (((x_x0) / f) * ((x_x0) * cos(k) - (y_y0) * sin(k)) + f * cos(k)) * cos(omiga); mtxA[2 * i][4] = -f * sin(k) - ((x_x0) / f) * ((x_x0) * sin(k) + (y_y0) * cos(k)); mtxA[2 * i][5] = (y_y0); mtxA[2 * i + 1][0] = (mtxR[1] * f + mtxR[2] * (y_y0)) / Ztg[i]; mtxA[2 * i + 1][1] = (mtxR[4] * f + mtxR[5] * (y_y0)) / Ztg[i]; mtxA[2 * i + 1][2] = (mtxR[7] * f + mtxR[8] * (y_y0)) / Ztg[i]; mtxA[2 * i + 1][3] = -(x_x0) * sin(omiga) - (((y_y0) / f) * ((x_x0) * cos(k) - (y_y0) * sin(k)) - f * sin(k)) * cos(omiga); mtxA[2 * i + 1][4] = -f * cos(k) - ((y_y0) / f) * ((x_x0) * sin(k) + (y_y0) * cos(k)); mtxA[2 * i + 1][5] = -(x_x0); } /***************输出矩阵A********************/ /* cout<<\"------------矩阵A-------------\\n\"; for(int i=0;i<8;++i) {for(int j=0;j<6;++j) {cout< Mxt_ZhuanZhi(&mtxA[0][0], &mtxAT[0][0], 2 * GCPN, 6); // 求A与AT的乘积，保存到mtxresult1矩阵 Mtx_XiangChen(&mtxAT[0][0], &mtxA[0][0], &mtxresult1[0][0], 6, 2 * GCPN, 6); // 求ATA的逆矩阵 Mtx_QiuNi(&mtxresult1[0][0], 6); 矿业工程学院_10级GIS_215 - 5 - 摄影测量—空间后方交会 // 求AT与L的乘积，结果保存到mtxresult2矩阵 Mtx_XiangChen(&mtxAT[0][0], mtxL, &mtxresult2[0], 6, 2 * GCPN, 1); // 求ATA的逆矩阵与ATL的乘积。即X=(ATA)-1 ATL Mtx_XiangChen(*mtxresult1, &mtxresult2[0], &mtxX[0], 6, 6, 1); //Xs=Xs+ΔXs mtxX即X矩阵 Ys，Zs同理 Xs += mtxX[0]; Ys += mtxX[1]; Zs += mtxX[2]; //ψ=ψ+Δψ mtxX即X矩阵 Κ，ω同理 fai += mtxX[3]; omiga += mtxX[4]; k += mtxX[5]; ///////////////输出每次迭代计算的改正数////////////////// cout << \"----------改正数--------\\n\" ; for(int i = 0; i < 6; ++i) { printf(\" mtxX[%d]=%lf\\n\ } } /*----------------旋转矩阵R-----------------------*/ mtxR[0] = cos(fai) * cos(k) - sin(fai) * sin(omiga) * sin(k); mtxR[1] = -cos(fai) * sin(k) - sin(fai) * sin(omiga) * cos(k); mtxR[2] = -sin(fai) * cos(omiga); mtxR[3] = cos(omiga) * sin(k); mtxR[4] = cos(omiga) * cos(k); mtxR[5] = -sin(omiga); mtxR[6] = sin(fai) * cos(k) + cos(fai) * sin(omiga) * sin(k); mtxR[7] = -sin(fai) * sin(k) + cos(fai) * sin(omiga) * cos(k); mtxR[8] = cos(fai) * cos(omiga); //输出结果到文件 FILE *fp_out; if(!(fp_out = fopen(\"空间后方交会计算结果.txt\错误，结果未能保存到文件\\n\"; /******************输出到文件******************/ fprintf(fp_out, \"\\n----------单像空间后方交会计算程序----------\"); fprintf(fp_out, \"\\n--------------解算完成---------------\"); fprintf(fp_out, \"\\n迭代次数为:%d\\n\ fprintf(fp_out, \"\\n像片外方位元素的解\\n\"); fprintf(fp_out, \"航向顷角ψ=%lf\\\n\ 矿业工程学院_10级GIS_215 - 6 - 摄影测量—空间后方交会 fprintf(fp_out, \"旁向倾角ω=%lf\\\n\ fprintf(fp_out, \"像片旋角κ=%lf\\\n\ fprintf(fp_out, \"\Xs=\%lf\\n\Ys=\%lf\\n\Zs=\%lf\\n\ fprintf(fp_out, \"\\n旋转矩阵R:\\n\"); fprintf(fp_out, \"%lf\\%lf\\%lf\\\n\ fprintf(fp_out, \"%lf\\%lf\\%lf\\\n\ fprintf(fp_out, \"%lf\\%lf\\%lf\\\n\ /******************输出到屏幕******************/ cout << \"\\n\\n----------单像空间后方交会计算程序 ----------\"; cout << \"\\n-----------------解算完成 ------------------\"; cout << \"\\n迭代次数为:\" << n << endl; printf(\"\\n像片外方位元素的解\\n\"); cout << \"航向顷角ψ=\\" << fai << endl; cout << \"旁向倾角ω=\\" << omiga << endl; cout << \"像片旋角κ=\\" << k << endl; printf(\"\Xs=\%lf\\n\Ys=\%lf\\n\Zs=\%lf\\n\ cout << endl; printf(\"旋转矩阵R:\\n\"); cout << mtxR[0] << '\' << mtxR[1] << '\' << mtxR[2] << endl; cout << mtxR[3] << '\' << mtxR[4] << '\' << mtxR[5] << endl; cout << mtxR[6] << '\' << mtxR[7] << '\' << mtxR[8] << endl << endl; /**************计算单位权中误差m0*******************/ for(i = 0; i < 2 * GCPN; i++) sum += pow(mtxL[i], 2); m0 = sqrt(sum / (2 * GCPN - 6)); /******************输出到文件******************/ fprintf(fp_out, \"\\n单位权中误差m0=%lf\\n\ fprintf(fp_out, \"\\n六个外方位元素精度:\\n\"); fprintf(fp_out, \"Xs的精度为:\%lf米\\n\ fprintf(fp_out, \"Ys的精度为:\%lf米\\n\ fprintf(fp_out, \"Zs的精度为:\%lf米\\n\ fprintf(fp_out, \" ψ的精度为:\%lf秒\\n\ fprintf(fp_out, \" ω的精度为:\%lf秒\\n\ fprintf(fp_out, \" κ的精度为:\%lf秒\\n\ fprintf(fp_out, \"\\n\\n----------解算完成----------\" \"-------------------------\"); /******************输出到屏幕******************/ cout << \"单位权中误差m0=\" << m0 << endl; cout << \"六个外方位元素精度:\" << endl; cout << \"Xs的精度为:\\" << m0 *sqrt(mtxresult1[0][0]) << \"米\" << endl; cout << \"Ys的精度为:\\" << m0 *sqrt(mtxresult1[1][1]) << \"米\" << endl; cout << \"Zs的精度为:\\" << m0 *sqrt(mtxresult1[2][2]) << \"米\" << endl; 矿业工程学院_10级GIS_215 - 7 - 摄影测量—空间后方交会 cout << \"ψ的精度为:\\" << m0 *sqrt(mtxresult1[3][3]) << \"秒\" << endl; cout << \"ω的精度为:\\" << m0 *sqrt(mtxresult1[4][4]) << \"秒\" << endl; cout << \"κ的精度为:\\" << m0 *sqrt(mtxresult1[5][5]) << \"秒\" << endl; cout << \"\\n\\n----------解算完成----------\\n\"; fclose(fp_out); fp_out = NULL; MessageBox(NULL, \"解算完成，点击确定打开结果保存文件！\空间后方交会\ system(\"空间后方交会计算结果.txt\"); //system(\"PAUSE\"); return 0; } /****************************************/ // GCPN是控制点对数，data用来保存控制点数据，m是比例尺，f,x0,y0是内方位元素 // 从程序所在文件夹读取“原始数据.txt\" int Inputdata(int &GCPN, GCP *data, double &m, double &f, double &x0, double &y0) { FILE *fp_data; if(!(fp_data = fopen(\"D:\\\\原始数据.txt\ { return 1; //如果文件读取失败，返回1 } cout << \"读取到的数据为\\n\"; fscanf(fp_data, \"%d%lf\ cout << \"控制点个数：\" << GCPN << endl; cout << \"航摄比例尺：\" << m << endl; if(GCPN < 4) { return 2; //控制点低于4个（不足），返回2 } fscanf(fp_data, \"%lf%lf%lf\ f /= 1000; //毫米单位换算为米 cout << \"主距f=\" << f << \"\x0=\" << x0 << \"\y0=\" << y0 << endl; if(m < 0 || f < 0) { return 3; //如果m<0或f<0，则返回3，即文本文件格式不对 } // data = new GCP [GCPN]; //分配控制点对数个GCP类型空间 // cout<<\"像点坐标\\\\"<<\"地面点坐标\\n\"<<\"x\\\"<<\"y\\\"<<\"Xt\\\"<<\"Yt\\\"<<\"Zt\\n\"; for(int i = 0; i < GCPN; ++i) { if(fscanf(fp_data, \"%lf%lf%lf%lf%lf\ 矿业工程学院_10级GIS_215 - 8 - 摄影测量—空间后方交会 &data[i].Zt) != 5) { return 3; //每对控制点5个数据不足，返回3 } data[i].x /= 1000; data[i].y /= 1000; //换算为米 // cout< /****************************************************/ // 如果比例尺为0，那么就用控制点计算比例尺 if(m == 0) { double image_d, ground_d, temp1, temp2, temp3; temp1 = pow(data[0].x - data[1].x, 2) + pow(data[0].y - data[1].y, 2); temp2 = pow(data[2].y - data[3].y, 2) + pow(data[2].y - data[3].y, 2); temp3 = pow(data[1].y - data[2].y, 2) + pow(data[1].y - data[2].y, 2); image_d = sqrt(temp1) + sqrt(temp2) + sqrt(temp3); temp1 = pow(data[0].Xt - data[1].Xt, 2) + pow(data[0].Yt - data[1].Yt, 2); temp2 = pow(data[2].Xt - data[3].Xt, 2) + pow(data[2].Yt - data[3].Yt, 2); temp3 = pow(data[1].Xt - data[2].Xt, 2) + pow(data[1].Yt - data[2].Yt, 2); ground_d = sqrt(temp1) + sqrt(temp2) + sqrt(temp3); m = 1 / (image_d / ground_d); cout << \"计算得到比例尺：\" << m << endl; } fclose(fp_data); fp_data = NULL; return 0; //一切顺利，返回0 } 结果截图： 矿业工程学院_10级GIS_215 - 9 - 摄影测量—空间后方交会 矿业工程学院_10级GIS_215 - 10 - 因篇幅问题不能全部显示，请点此查看更多更全内容