摄影测量单像后方交会实验报告
2024-05-16 06:50:05
摄影测量单像后方交会实验报告
文章目录
- 摄影测量单像后方交会实验报告
- 一、理论基础
- 二、问题
- 三、源程序
- 四、说明
- 1.初值
- 2.迭代次数
- 3.结果
- 4.讨论
- 4.1关于输出结果的有效数字位数：
- 4.2矩阵求逆：
- 4.3优化：
一、理论基础
外方位元素是确定摄影瞬间像片在地面直角坐标系中空间位置和姿态的参数，根据6个外方位元素，可以恢复像片与目标间的位置关系，重建目标的模型，其重要性不言而喻。
单像后方交会，即根据像片覆盖范围内一定数量的分布合理的已知坐标的地面控制点，利用共线方程求解像片的外方位元素。
二、问题
三、源程序
#include<iostream>
#include<iomanip>
#include<string.h>
#include<math.h>using namespace std;void PrintArray(double* a, int b, int c); //输出矩阵
void TranspositionArray(double* a, double* aT, int b, int c); //转置矩阵
void MultiplyArray(double* a, double* b, double* c, int m, int n, int l); //矩阵相乘
void make2array(double**& a, int n); //矩阵扩充（进行初等变换）
void deletarray(double**& a, int n); //释放内存
int ConverseArray(double* ip, double* rp, int n); //矩阵求逆
void SubtratArray(double* a, double* b, double* c, int m, int n); //矩阵相减int main()
{//像点坐标double ICP[4][2] = { -86.15,-68.99,-53.40,82.21,-14.78,-76.63,10.46,64.43 };//地面点坐标double GCP[4][3] = { 36589.41, 25273.32,2195.17,37631.08,31324.51,728.69,39100.97,24934.98,2386.50,40426.54,30319.81,757.31 };//估算摄影比例尺double scale = (GCP[0][0] - GCP[1][0]) / (ICP[0][0] - ICP[1][0])*1000;double f = 0.15324;//统一长度单位为米for (int i = 0; i < 4; i++){for (int j = 0; j < 2; j++){ICP[i][j] /= 1000.0;}}//定义外方位元素double Xs = 0.0, Ys = 0.0, Zs = 0.0, FAI = 0.0, OMEGA = 0.0, KAPPA = 0.0;//外方位元素的改正数double delta[6] = { 0.0 };//定义旋转矩阵double Rotate[3][3] = { 0.0 };//确定Xs,YS,Zs初值for (int i = 0; i < 4; i++){Xs = (Xs + GCP[i][0]) / 4.0;Ys = (Ys + GCP[i][1]) / 4.0;}Zs = f * scale;//输出外方位元素的初始值cout << "外方位元素的初始值为：" << endl;cout <<"Xs初值："<< Xs << " " <<"Ys初值：" <<Ys << " " << "Zs初值："<<Zs << endl;cout << "FAI初值："<<FAI << " " << "OMEGA初值："<<OMEGA << " " <<"KAPPA初值："<< KAPPA << endl;//定义double x[4] = { 0 }, y[4] = { 0 }, A[8][6] = { 0 }, AT[6][8] = { 0 }, ATA[6][6] = { 0 }, ATAInv[6][6] = { 0 }, L[8] = { 0 }, ATAAT[6][8] = { 0 };double v[8] = { 0 }, m0 = 0, m[6] = { 0 }, AX[8] = { 0 }, vv = 0;int n = 0;//累计循环次数n//主循环while (1){//求旋转矩阵Rotate[0][0] = cos(FAI) * cos(KAPPA) - sin(FAI) * sin(OMEGA) * sin(KAPPA); //a1Rotate[0][1] = (-1.0) * cos(FAI) * sin(KAPPA) - sin(FAI) * sin(OMEGA) * cos(KAPPA); //a2Rotate[0][2] = (-1.0) * sin(FAI) * cos(OMEGA); //a3Rotate[1][0] = cos(OMEGA) * sin(KAPPA); //b1Rotate[1][1] = cos(OMEGA) * cos(KAPPA); //b2Rotate[1][2] = (-1.0) * sin(OMEGA); //b3Rotate[2][0] = sin(FAI) * cos(KAPPA) + cos(FAI) * sin(OMEGA) * sin(KAPPA); //c1Rotate[2][1] = (-1.0) * sin(FAI) * sin(KAPPA) + cos(FAI) * sin(OMEGA) * cos(KAPPA); //c2Rotate[2][2] = cos(FAI) * cos(OMEGA); //c3for (int i = 0; i < 4; i++){//共线方程x[i] = (-1.0) * f * (Rotate[0][0] * (GCP[i][0] - Xs) + Rotate[1][0] * (GCP[i][1] - Ys)+ Rotate[2][0] * (GCP[i][2] - Zs)) / (Rotate[0][2] * (GCP[i][0] - Xs)+ Rotate[1][2] * (GCP[i][1] - Ys) + Rotate[2][2] * (GCP[i][2] - Zs));y[i] = (-1.0) * f * (Rotate[0][1] * (GCP[i][0] - Xs) + Rotate[1][1] * (GCP[i][1] - Ys)+ Rotate[2][1] * (GCP[i][2] - Zs)) / (Rotate[0][2] * (GCP[i][0] - Xs)+ Rotate[1][2] * (GCP[i][1] - Ys) + Rotate[2][2] * (GCP[i][2] - Zs));//矩阵LL=像点坐标-共线方程计算值L[i * 2] = ICP[i][0] - x[i];L[i * 2 + 1] = ICP[i][1] - y[i];//偏导数矩阵A[i * 2][0] = (-1.0) * f / (Zs - GCP[i][2]); //a11A[i * 2][1] = 0.0; //a12A[i * 2][2] = (-1.0) * x[i] / (Zs - GCP[i][2]); //a13A[i * 2][3] = (-1.0) * f * (1 + (x[i] * x[i]) / (f * f)); //a14A[i * 2][4] = (-1.0) * x[i] * y[i] / f; //a15A[i * 2][5] = y[i]; //a16A[i * 2 + 1][0] = 0.0; //a21A[i * 2 + 1][1] = A[i * 2][0]; //a22A[i * 2 + 1][2] = (-1.0) * y[i] / (Zs - GCP[i][2]); //a23A[i * 2 + 1][3] = A[i * 2][4]; //a24A[i * 2 + 1][4] = (-1.0) * f * (1 + (y[i] * y[i]) / (f * f)); //a25A[i * 2 + 1][5] = (-1.0) * x[i]; //a26}//求矩阵A的转置TranspositionArray(*A, *AT, 8, 6);//矩阵A的转置与矩阵A相乘MultiplyArray(*AT, *A, *ATA, 6, 8, 6);//矩阵A的转置与矩阵A相乘的逆ConverseArray(*ATA, *ATAInv, 6);//矩阵A的转置与矩阵A相乘的逆与A的转职相乘MultiplyArray(*ATAInv, *AT, *ATAAT, 6, 6, 8);//求得误差方程的解MultiplyArray(*ATAAT, L, delta, 6, 8, 1);Xs += delta[0];Ys += delta[1];Zs += delta[2];FAI += delta[3];OMEGA += delta[4];KAPPA += delta[5];++n;if ((fabs(delta[0]) <1e-4 ) && (fabs(delta[1]) < 1e-4) && (fabs(delta[2]) < 1e-4)&&n<1000){break;}else if(n>=1000){cout << "1000次循环未解出来" << endl;}}cout << "总计" << n << "次循环" << endl;MultiplyArray(*A, delta, AX, 8, 6, 1);SubtratArray(AX, L, v, 8, 1);for (int i = 0; i < 8; i++){vv += v[i] * v[i];}m0 = sqrt(vv / 2);for (int i = 0; i < 6; i++){m[i] = m0 * sqrt(fabs(ATAInv[i][i]));}cout << "相片的外方位元素为：" << endl;cout << "Xs=" << setprecision(7) << Xs << endl;cout << "Ys=" << setprecision(7) << Ys << endl;cout << "Zs=" << setprecision(7) << Zs << endl;cout << "FAI=" << FAI << endl;cout << "OMEGA=" << OMEGA << endl;cout << "KAPPA=" << KAPPA << endl;cout << "旋转矩阵R为：" << endl;PrintArray(*Rotate, 3, 3);return 0;
}//输出矩阵
void PrintArray(double* a, int b, int c)
{ for (int i = 0; i < b; i++) {for (int j = 0; j < c; j++) {cout << a[i * c + j] << " ";}cout << endl;}
}//转置矩阵
void TranspositionArray(double* a, double* aT, int b, int c)
{for (int i = 0; i < b; i++) {for (int j = 0; j < c; j++) {aT[j * b + i] = a[i * c + j];}}
}//矩阵a，b相乘
void MultiplyArray(double* a, double* b, double* c, int m, int n, int l)
{for (int i = 0; i < m * l; i++) //矩阵初始化为0{c[i] = 0;}for (int i = 0; i < m; i++) {for (int j = 0; j < l; j++) {for (int k = 0; k < n; k++) {c[i * l + j] += a[i * n + k] * b[k * l + j];}}}
}//把n阶方阵扩充为n×2n矩阵并初始化为0
void make2array(double**& a, int n)
{int i, j;a = new double* [n];for (i = 0; i < n; i++){a[i] = new double[2 * n];}for (i = 0; i < n; i++){for (j = 0; j < 2 * n; j++){a[i][j] = 0;}}
}//释放内存
void deletarray(double**& a, int n)
{int i;for (i = 0; i < n; i++){delete[]a[i];}delete[]a;
}//利用初等行变换法进行矩阵求逆
int ConverseArray(double* ip, double* rp, int n)
{double** mat = NULL; //做行变换的矩阵 int i, j, r;double k, temp;//初始化mat为0,大小为n*2nmake2array(mat, n);for (i = 0; i < n; i++){for (j = 0; j < n; j++){mat[i][j] = ip[i * n + j];}mat[i][n + i] = 1; //初始化右侧的单位矩阵 }//做初等行变换化为上三角阵 for (i = 0; i < n; i++){if (mat[i][i] == 0) //第i行i列为0 {for (j = i + 1; j < n; j++){if (mat[j][i] != 0) //找一个非0行 {for (r = i; r < 2 * n; r++){temp = mat[j][r];mat[j][r] = mat[i][r];mat[i][r] = temp;}break; //跳出j循环 }}}if (mat[i][i] == 0) return 0; //行列式为0则返回 for (j = i + 1; j < n; j++){if (mat[j][i] != 0) //i行i列下方的j行i列元素不为0 {k = -mat[j][i] / mat[i][i]; //做行变换 for (r = i; r < 2 * n; r++)mat[j][r] = mat[j][r] + k * mat[i][r];}}}//化成单位矩阵 for (i = n - 1; i >= 0; i--){k = mat[i][i];for (r = i; r < 2 * n; r++)mat[i][r] = mat[i][r] / k;for (j = i - 1; j >= 0; j--){k = -mat[j][i];for (r = i; r < 2 * n; r++)mat[j][r] = mat[j][r] + k * mat[i][r];}}//将结果输出 for (i = 0; i < n; i++)for (j = 0; j < n; j++)rp[i * n + j] = mat[i][j + n];//mat矩阵n释放 deletarray(mat, n);return 1;
}//矩阵相减
void SubtratArray(double* a, double* b, double* c, int m, int n)
{int i, j;for (i = 0; i < m; i++){for (j = 0; j < n; j++){c[i * n + j] = a[i * n + j] - b[i * n + j];}}
}
四、说明
1.初值
代码如下（示例）：
double scale = (GCP[0][0] - GCP[1][0]) / (ICP[0][0] - ICP[1][0])*1000;for (int i = 0; i < 4; i++){Xs = (Xs + GCP[i][0]) / 4.0;Ys = (Ys + GCP[i][1]) / 4.0;}Zs = f * scale;
Xs、Ys的初值为四个物点坐标求均值，Zs的初值为f乘以比例尺。本题目中比例尺未给出，由已知数据估算：两个地面控制点的x坐标差除以对应像点的x坐标差。初值的设定关系着外方位元素改正数的收敛性。如果初值不合理，计算结果发散，则会超出设定的循环次数（本程序为1000次），程序终止。
2.迭代次数
if ((fabs(delta[0]) <1e-4 ) && (fabs(delta[1]) < 1e-4) && (fabs(delta[2]) < 1e-4)&&n<1000){break;}else if(n>=1000){cout << "1000次循环未解出来" << endl;}
迭代次数与计算结果精度有关。本程序结果精度为1e-3，如果要继续提高精度，可以适当减小Xs、Ys、Zs的改正数。
本程序的迭代次数为：
如果改为
if ((fabs(delta[0]) <1e-8 ) && (fabs(delta[1]) < 1e-8) && (fabs(delta[2]) < 1e-8)&&n<1000){break;}else if(n>=1000){cout << "1000次循环未解出来" << endl;}
迭代次数为：
因此随着精度要求的提高，迭代次数增加。
3.结果
本程序的运行结果如下：
可见结果与题目提示的结果较好地符合。
4.讨论
4.1关于输出结果的有效数字位数：
cout << "Xs=" << setprecision(7) << Xs << endl;cout << "Ys=" << setprecision(7) << Ys << endl;cout << "Zs=" << setprecision(7) << Zs << endl;
通过改变setprecision（）里的值，可以改变输出结果的有效数字位数。本程序为保留七位有效数字。
4.2矩阵求逆：
本程序的矩阵求逆采用初等行变换法进行求逆，也可采用伴随矩阵法求逆。
4.3优化：
本程序的矩阵运算全部采用循环的方法进行计算，可以引入有关矩阵运算库进行简化。
摄影测量单像后方交会实验报告相关推荐
- 计算机公共基础知识实验报告,MIPS单周期CPU实验报告总结.doc
<计算机组成原理实验> 实验报告 (实验二) 学 院 名 称 : 专业(班级) : 学 生 姓 名 : 学号 : 时间:2017年11月25日 成 绩 : 实 验 二 : 单周期 CPU设 ...
- 《数据结构》单链表操作实验报告
用于复习链表,自己写的有关链表的一些操作,包括后插前插创建链表,删除,插入,有序合并,查找,遍历,反转,查中间节点,删倒数节点,检测环结构的功能.可以用做C语言实训参考. #include<io ...
- 大学物理实验长度的测量实验报告_(完整精品)大学物理实验报告之长度基本测量.doc...
(完整精品)大学物理实验报告之长度基本测量.doc 大学物理实验报告 评 分姓名 学号 评 分 学院 班级 实验日期 2017 年5 月23 日 实验地点:实验楼B411室 实验名称 长度的基本测量 ...
- c语言综合实验报告与材料专业,实验报告汇总
种子萌发的实验报告 一.做实验1.材料工具(1)常见的种子(如:绿豆 黄豆)40粒.(2)有盖的罐头4个,小勺1个,餐巾纸8张,4张分别标有1.2.3.4的标签,胶水,清水.2．方法步骤(1)在第一个 ...
- 计算机视觉与摄影测量的不同2
计算机视觉和摄影测量两套理论分别在上个世纪60年代提出,经过计算机技术和数码相机技术数十年的迅猛发展,从一开始的相去甚远到今日的水乳交融,从理论到实践都发生了巨大的变化.直到今天,两门学科已经彼此紧密 ...
- 液位单闭环实验计算机控制,过程控制实验报告3(液位单闭环实验)
<过程控制实验报告3(液位单闭环实验)>由会员分享,可在线阅读,更多相关<过程控制实验报告3(液位单闭环实验)(2页珍藏版)>请在人人文库网上搜索. 1.班级:082班 座号: ...
- 单容水箱液位pid控制实验报告_实验二(单容水箱液位pid控制实验)实验报告电子版.doc...
实验二(单容水箱液位pid控制实验)实验报告电子版 电子科技大学中山学院学生实验报告 系别: 机电工程学院 专业: 自动化 课程名称:过程控制与自动化仪表 班级: 自动化 姓名: 学号: 组别: 实验 ...
- 有滞后单容对象MATLAB仿真,自动化生产线实训实验报告
<自动化生产线实训实验报告>由会员分享,可在线阅读,更多相关<自动化生产线实训实验报告(47页珍藏版)>请在人人文库网上搜索. 1.北京科技大学自动化生产线实训实验报告班 级: ...
- 哈工大 c语言测试与系统控制 ad,哈工大——c语言在测量与控制中应用实验报告.pdf...
哈工大--c语言在测量与控制中应用实验报告 Harbin Institute of Technology Harbin Institute of Technology C 语言在测量与控制中的 C 语 ...
- C++——《算法分析与设计》实验报告——单源最短路径问题
实验名称: 单源最短路径问题 实验地点: 实验目的: 1. 理解分支限界法的剪枝搜索策略: 2. 掌握分支限界法的算法柜架: 3. 掌握分支限界法的算法步骤: 4. 通过应用范例学习动态规划算 ...
最新文章
- 脚本1)启动jetty的脚本
- ‘torch.nn‘ has no attribute ‘SiLU‘
- 用 GDI 操作 EMF 文件[8]: 绘制图元文件时改变画笔与画刷
- python读excel乱码_Python读写excel练习_去除excel中乱码行，并添加列
- [转载]项目经理必备工具包：项目管理中的22个思维导图
- c++STL容器的Map和multimap
- 九九乘法表-九九乘法表数据输出
- MapReduce原理全剖析
- 论文小综 | 知识图谱表示学习中的零样本实体研究
- java spring mvc json ajax 优势_SpringMVC后台json数据前台ajax获取不到！！！急求解答！！！...
- phpmyadmin4.8.1远程文件包含漏洞
- SQLi LABS Less 27a 联合注入+布尔盲注+时间盲注
- 深度学习面试100题
- 圣诞节苹果服务器没有人维护2020,2020圣诞节真的推迟到1月8号吗
- bzoj4199品酒大会
- 笔记 ：AVS2背景建模
- BlackHat DEFCON现场秀：阿里安全专家演示“视频水印叠加”和“一分钟越狱iOS 11.4”...
- Http协议和Python调试过程
- vue3+vant开发微信公众号网页爬坑不完全指北
- PostgreSQL集群方案-Postgres-XL