共线方程求解外方位元素--单片空间后方交会

单片空间后方交会是利用三个不在一条直线的已知点计算相片外方位元素的方法。

获取已知数据，包括：
- x[]、y[]：改正后的控制点的像方坐标
- X[]、Y[]、Z[]：控制点的物方坐标
- f: 相机焦距
  在空间后方交会中，通常采用像片四个角上（或者更多控制点），同时采用最小二乘原理进行平差计算。

确定未知数的初始值
在竖直摄影情况下，三个角元素初值都为0，三个线元素初始值：

/************************************************************************/
/*  外方位元素的初值    */
/************************************************************************/
for (int i = 0; i < 4; i++)
{Xs += X[i];Ys += Y[i];Zs += Z[i];}
Xs = Xs / 4;
Ys = Ys / 4;
Zs = Zs / 4 + H;
a1 = 0.0f;
a2 = 0.0f;
a3 = 0.0f;

进入循环，角度改正数小于指定值(一般为0.1)退出循环
计算旋转矩阵R
```
computeRMatrix(a1, a2, a3, mR);
```

逐点计算像点坐标近似值

/**
*计算像平面的近似坐标
*/for (int i = 0; i < 4; i++){_x[i] = -f*(R[0][0] * (X[i] - Xs) + R[1][0] * (Y[i] - Ys) + R[2][0] * (Z[i] - Zs))/ (R[0][2] * (X[i] - Xs) + R[1][2] * (Y[i] - Ys) + R[2][2] * (Z[i] - Zs));_y[i] = -f*(R[0][1] * (X[i] - Xs) + R[1][1] * (Y[i] - Ys) + R[2][1] * (Z[i] - Zs))/ (R[0][2] * (X[i] - Xs) + R[1][2] * (Y[i] - Ys) + R[2][2] * (Z[i] - Zs));}

组成误差方程式
每个点的误差方程式：

误差方程的系数矩阵L

/************************************************************************/
/* 误差方程的系数矩阵L                                                                  */
/************************************************************************/
    for (int i = 0; i < 4; i++){        L[2 * i] = x[i] - _x[i];
        L[2 * i + 1] = y[i] - _y[i];
    }

首先求解TZ：

float TZ = R[0][2] * (X[i / 2] - Xs) + R[1][2] * (Y[i / 2] - Ys) + R[2][2] * (Z[i / 2] - Zs);

偏导XS求解（YS、ZS类似）

偏导phi求解（偏导omega、偏导kapa类似）

根据平差原理，法方程为：

这里采用高斯消元法进行求解，也可以采用逆矩阵求解。

/************************************************************************/
/* 误差方程的系数矩阵A*/
/************************************************************************/
    for (int i = 0; i < 8; i = i + 2)
    {        float TZ = R[0][2] * (X[i / 2] - Xs) + R[1][2] * (Y[i / 2] - Ys) + R[2][2] * (Z[i / 2] - Zs);
        A[i][0] = 1 / TZ*(R[0][0] * f + R[0][2] * _x[i / 2]); mA.write(i, 0, A[i][0]);
        A[i][1] = 1 / TZ*(R[1][0] * f + R[1][2] * _x[i / 2]); mA.write(i, 1, A[i][1]);
        A[i][2] = 1 / TZ*(R[2][0] * f + R[2][2] * _x[i / 2]); mA.write(i, 2, A[i][2]);        A[i][3] = _y[i / 2] * sin(omega) - (_x[i / 2] / f*(_x[i / 2] * cos(kapa) - _y[i / 2] * sin(kapa)) + f*cos(kapa))*cos(omega); mA.write(i, 3, A[i][3]);
        A[i][4] = -f*sin(kapa) - _x[i / 2] / f*(_x[i / 2] * sin(kapa) + _y[i / 2] * cos(kapa));; mA.write(i, 4, A[i][4]);
        A[i][5] = _y[i / 2]; mA.write(i, 5, A[i][5]);
    }

更新角度值

    phi += Dx.read(3, 0);omega += Dx.read(4, 0);kapa += Dx.read(5, 0);

工程地址：共线方程求解外方位元素–单片空间后方交会

main.cpp

#include "_Matrix.h"
#include <iostream>
#include <cmath>
#include <fstream>
using namespace std;_Matrix_Calc matixCalc;
_Matrix mA(8, 6), mAT(6, 8), mATA(6, 6), mATL(6, 1), mL(8, 1), Dx(6, 1), mR(3, 3);
float f = 0.15324f;
float m = 0;
float x[] = { -0.08615f, -0.05340f, -0.01478f, 0.01046f };
float y[] = { -0.06899f, 0.08221f, -0.07663f, 0.06443f };
float X[] = { 36589.41f, 37631.08f, 39100.97f, 40426.54f };
float Y[] = { 25273.32f, 31324.51f, 24934.98f, 30319.81f };
float Z[] = { 2195.17f, 728.69f, 2386.50f, 757.31f };
float Xs, Ys, Zs, dXs, dYs, dZs;
float phi, omega, kapa, dphi, domega, dkapa;
float H;
float R[3][3];
float _x[4], _y[4];
float A[8][6];
float L[8];void computeRMatrix(float&phi, float&omega, float&kapa, _Matrix&m){R[0][0] = cos(phi)*cos(kapa) - sin(phi)*sin(omega)*sin(kapa); m.write(0, 0, R[0][0]);R[0][1] = -cos(phi)*sin(kapa) - sin(phi)*sin(omega)*cos(kapa); m.write(0, 1, R[0][1]);R[0][2] = -sin(phi)*cos(omega); m.write(0, 2, R[0][2]);R[1][0] = cos(omega)*sin(kapa); m.write(1, 0, R[1][0]);R[1][1] = cos(omega)*cos(kapa); m.write(1, 1, R[1][1]);R[1][2] = -sin(omega); m.write(1, 2, R[1][2]);R[2][0] = sin(phi)*cos(kapa) + cos(phi)*sin(omega)*sin(kapa); m.write(2, 0, R[2][0]);R[2][1] = -sin(phi)*sin(kapa) + cos(phi)*sin(omega)*cos(kapa); m.write(2, 1, R[2][1]);R[2][2] = cos(phi)*cos(omega); m.write(2, 2, R[2][2]);
}int main(){mA.init_matrix();mAT.init_matrix();mATA.init_matrix();mATL.init_matrix();mL.init_matrix();Dx.init_matrix();mR.init_matrix();float temp1 = 0, temp2 = 0;for (int i = 0; i < 3; i++){temp1 += sqrt(pow(x[i] - x[i + 1], 2) + pow(y[i] - y[i + 1], 2));temp2 += sqrt(pow(X[i] - X[i + 1], 2) + pow(Y[i] - Y[i + 1], 2));}m = temp2 / temp1;H = m*f;for (int i = 0; i < 4; i++){Xs += X[i];Ys += Y[i];Zs += Z[i];}Xs = Xs / 4;Ys = Ys / 4;Zs = Zs / 4 + H;phi = 0.0f;omega = 0.0f;kapa = 0.0f;int num = 0;while (num == 0 || abs(Dx.read(3, 0)) * 206265 > 0.1 || abs(Dx.read(4, 0)) * 206265 > 0.1 || abs(Dx.read(5, 0)) * 206265 > 0.1){computeRMatrix(phi, omega, kapa, mR);for (int i = 0; i < 4; i++){_x[i] = -f*(R[0][0] * (X[i] - Xs) + R[1][0] * (Y[i] - Ys) + R[2][0] * (Z[i] - Zs))/ (R[0][2] * (X[i] - Xs) + R[1][2] * (Y[i] - Ys) + R[2][2] * (Z[i] - Zs));_y[i] = -f*(R[0][1] * (X[i] - Xs) + R[1][1] * (Y[i] - Ys) + R[2][1] * (Z[i] - Zs))/ (R[0][2] * (X[i] - Xs) + R[1][2] * (Y[i] - Ys) + R[2][2] * (Z[i] - Zs));}for (int i = 0; i < 4; i++){L[2 * i] = x[i] - _x[i];L[2 * i + 1] = y[i] - _y[i];}mL.arr = L;for (int i = 0; i < 8; i = i + 2){float TZ = R[0][2] * (X[i / 2] - Xs) + R[1][2] * (Y[i / 2] - Ys) + R[2][2] * (Z[i / 2] - Zs);A[i][0] = 1 / TZ*(R[0][0] * f + R[0][2] * _x[i / 2]); mA.write(i, 0, A[i][0]);A[i][1] = 1 / TZ*(R[1][0] * f + R[1][2] * _x[i / 2]); mA.write(i, 1, A[i][1]);A[i][2] = 1 / TZ*(R[2][0] * f + R[2][2] * _x[i / 2]); mA.write(i, 2, A[i][2]);A[i][3] = _y[i / 2] * sin(omega) - (_x[i / 2] / f*(_x[i / 2] * cos(kapa) - _y[i / 2] * sin(kapa)) + f*cos(kapa))*cos(omega); mA.write(i, 3, A[i][3]);A[i][4] = -f*sin(kapa) - _x[i / 2] / f*(_x[i / 2] * sin(kapa) + _y[i / 2] * cos(kapa));; mA.write(i, 4, A[i][4]);A[i][5] = _y[i / 2]; mA.write(i, 5, A[i][5]);}for (int i = 1; i < 8; i = i + 2){float TZ = R[0][2] * (X[(i - 1) / 2] - Xs) + R[1][2] * (Y[(i - 1) / 2] - Ys) + R[2][2] * (Z[(i - 1) / 2] - Zs);A[i][0] = 1 / TZ*(R[0][1] * f + R[0][2] * _y[(i - 1) / 2]); mA.write(i, 0, A[i][0]);A[i][1] = 1 / TZ*(R[1][1] * f + R[1][2] * _y[(i - 1) / 2]); mA.write(i, 1, A[i][1]);A[i][2] = 1 / TZ*(R[2][1] * f + R[2][2] * _y[(i - 1) / 2]); mA.write(i, 2, A[i][2]);A[i][3] = _x[(i - 1) / 2] * sin(omega) - (_y[(i - 1) / 2] / f*(_x[(i - 1) / 2] * cos(kapa) - _y[(i - 1) / 2] * sin(kapa)) - f*sin(kapa))*cos(omega); mA.write(i, 3, A[i][3]);A[i][4] = -f*cos(kapa) - _y[(i - 1) / 2] / f*(_x[(i - 1) / 2] * sin(kapa) + _y[(i - 1) / 2] * cos(kapa)); mA.write(i, 4, A[i][4]);A[i][5] = -_x[(i - 1) / 2]; mA.write(i, 5, A[i][5]);}matixCalc.transpos(&mA, &mAT);matixCalc.multiply(&mAT, &mA, &mATA);matixCalc.multiply(&mAT, &mL, &mATL);matixCalc.Gauss_Col(&mATA, &mATL, &Dx);Xs += Dx.read(0, 0);Ys += Dx.read(1, 0);Zs += Dx.read(2, 0);phi += Dx.read(3, 0);omega += Dx.read(4, 0);kapa += Dx.read(5, 0);num++;}printf("一共迭代了%d次\n", num);printf("外方位元素\n");printf("%f  %f  %f \n", Xs, Ys, Zs);printf("%f  %f  %f \n", phi, omega, kapa);printf("旋转矩阵\n");mR.print();return 0;
}

共线方程求解外方位元素--单片空间后方交会相关推荐

空间前方交会（利用相机外方位元素和像点坐标进行解算）
目录一.前言二.空间前方交会 1. 前方交会的概念 2. 基本公式三.代码展示四.小结一.前言在摄影测量过程中,得到相机的外方位元素以及地面控制点对应的像点坐标之后,如何解算地面控制点 ...
Matlab解算空间后方交会外方位元素
目录 1 问题描述 2 实现部分参考文献本问题为武汉大学摄影测量学教材课后习题,现在用MATLAB实现,供大家学习参考. 1 问题描述已知摄像机主距f=153.24mm,四对点的像点坐标 ...
立体像对空间前方交会（利用外方位元素交会出地面点三维坐标）
要通过外方位元素进行前方交会首先需要的已知量有: 然后计算左右片在地辅坐标系中旋转矩阵的方向余弦再计算基线分量计算像点的像空间辅助坐标计算投影系数计算地面点的左像辅系坐标计算地面点的地面坐标 ...
摄影测量外方位元素代码
航摄相片的外方位元素表示的是摄影摄影瞬间相片上的点对于地面上的点之间的关系的一些参数,在测绘工作中,如果求出了一张航拍相片的外方位元素,那么就可以根据像素点的坐标计算出对应的地面点的坐标,而这些解算过 ...
航飞原始影像外方位元素_【技术】无人机倾斜摄影建模技术在虚拟现实中的应用...
(如有侵权,请联系删除) 摘要针对于虚拟现实平台中构建三维场景的费时费力问题,基于无人机倾斜摄影建模技术构建三维模型,利用 3DS Max 建模软件进行模型优化,并结合 Unity 3D 引擎构建 ...
外参矩阵（旋转矩阵+平移向量）以及外方位元素的关系
外参包括旋转矩阵R3×3.平移向量T3×1,它们共同描述了如何把点从世界坐标系转换到摄像机坐标系,旋转矩阵描述了世界坐标系的坐标轴相对于摄像机坐标轴的方向, 平移向量描述了在摄像机坐标系下空间原点的位 ...
航飞原始影像外方位元素_影响无人机航测精度的因素都有哪些？
‍[摘要]本文通过对1∶500 无人机航测法成图过程中误差产生的来源进行分析,研究‍提高成图精度的关键技术,经过试验建立1∶500 无人机航测法高精度成图技术路线和工艺流程,并给出实际生产项目中的具体 ...
单片空间后方交会程序设计（代码共享）
单片空间后方交会程序设计 1 目的用程序设计语言(VC或者VB)编写一个完整的单片空间后方交会程序,通过对提供的试验数据进行计算,输出像片的外方位元素并评定精度.本实验的目的在于让学生深入理解单片空 ...
航飞原始影像外方位元素_浅谈大型倾斜航摄仪（飞思）的数据处理流程
引言:虽然,现阶段无人机倾斜作业盛行,但对于大面积.特殊敏感区域上空,能及时有效获取外业倾斜数据的问题上依旧少不了大型倾斜航摄仪的身影.对于此类型数据处理,M3D同样适用.本文以飞思航摄仪数据进行处理 ...
用C#编写摄影测量后方交会求外方位元素
程序源码程序源码下载地址 https://download.csdn.net/download/u011713916/11743497 实验原理 **详细题目: ** 编程思想: 具体代码我使用了M ...

共线方程求解外方位元素--单片空间后方交会

共线方程求解外方位元素--单片空间后方交会相关推荐

最新文章

热门文章