一、前言

在摄影测量过程中，得到相机的外方位元素以及地面控制点对应的像点坐标之后，如何解算地面控制点坐标？
在研究过程中，可以利用空间误差处理中的平差知识以及计算机编程进行实现。

二、空间前方交会

1. 前方交会的概念

利用立体像对两张像片的内方位元素、同名像点坐标和像对的相对方位元素(或外方位元素)解算模型点坐标(或地面点坐标)的工作，称为空间前方交会。在摄影测量中主要有两种：
(1). 利用立体像对两张像片的相对方位元素，计算模型点的三位坐标。
(2). 利用立体像对两张像片的外方位元素，计算地面点的地面坐标。
而在本文章中，将重点介绍在得到相机的外方位元素以及地面控制点对应的像点坐标，利用平差的知识，求算出地面控制点的坐标。

2. 基本公式

由于共线方程大家都已知，就不在这里赘述。这里主要列出误差方程、误差系数以及旋转矩阵。

误差方程
v x 1 = ∂ x 1 ∂ X Δ X + ∂ x 1 ∂ Y Δ Y + ∂ x 1 ∂ Z Δ Z − l x 1 v y 1 = ∂ y 1 ∂ X Δ X + ∂ y 1 ∂ Y Δ Y + ∂ y 1 ∂ Z Δ Z − l y 1 v x 2 = ∂ x 2 ∂ X Δ X + ∂ x 2 ∂ Y Δ Y + ∂ y 2 ∂ Z Δ Z − l x 2 v y 2 = ∂ y 2 ∂ X Δ X + ∂ y 2 ∂ Y Δ Y + ∂ y 2 ∂ Z Δ Z − l y 2 ⋯ v x n = ∂ x n ∂ X Δ X + ∂ x n ∂ Y Δ Y + ∂ x n ∂ Z Δ Z − l x n v y n = ∂ y n ∂ X Δ X + ∂ y n ∂ Y Δ Y + ∂ y n ∂ Z Δ Z − l y n \begin{align} v_{x1} &= \frac{\partial x_1}{\partial X}\Delta X+\frac{\partial x_1}{\partial Y}\Delta Y+\frac{\partial x_1}{\partial Z}\Delta Z-l_{x1}\\ v_{y1} &= \frac{\partial y_1}{\partial X}\Delta X+\frac{\partial y_1}{\partial Y}\Delta Y+\frac{\partial y_1}{\partial Z}\Delta Z-l_{y1}\\ v_{x2} &= \frac{\partial x_2}{\partial X}\Delta X+\frac{\partial x_2}{\partial Y}\Delta Y+\frac{\partial y_2}{\partial Z}\Delta Z-l_{x2}\\ v_{y2} &= \frac{\partial y_2}{\partial X}\Delta X+\frac{\partial y_2}{\partial Y}\Delta Y+\frac{\partial y_2}{\partial Z}\Delta Z-l_{y2}\\ &\dotsb\nonumber \\ v_{xn} &= \frac{\partial x_n}{\partial X}\Delta X+\frac{\partial x_n}{\partial Y}\Delta Y+\frac{\partial x_n}{\partial Z}\Delta Z-l_{xn} \tag{2n-1}\\ v_{yn} &= \frac{\partial y_n}{\partial X}\Delta X+\frac{\partial y_n}{\partial Y}\Delta Y+\frac{\partial y_n}{\partial Z}\Delta Z-l_{yn} \tag{n} \end{align} vx1vy1vx2vy2vxnvyn=∂X∂x1ΔX+∂Y∂x1ΔY+∂Z∂x1ΔZ−lx1=∂X∂y1ΔX+∂Y∂y1ΔY+∂Z∂y1ΔZ−ly1=∂X∂x2ΔX+∂Y∂x2ΔY+∂Z∂y2ΔZ−lx2=∂X∂y2ΔX+∂Y∂y2ΔY+∂Z∂y2ΔZ−ly2⋯=∂X∂xnΔX+∂Y∂xnΔY+∂Z∂xnΔZ−lxn=∂X∂ynΔX+∂Y∂ynΔY+∂Z∂ynΔZ−lyn(2n-1)(n)
误差系数
∂ x ∂ X = − 1 Z ( a 1 f x + a 3 x ) ∂ x ∂ Y = − 1 Z ( b 1 f x + b 3 x ) ∂ x ∂ Z = − 1 Z ( c 1 f x + c 3 x ) ∂ y ∂ X = − 1 Z ( a 2 f y + a 3 y ) ∂ y ∂ Y = − 1 Z ( b 2 f y + b 3 y ) ∂ y ∂ Z = − 1 Z ( c 2 f y + c 3 y ) \begin{align} \frac{\partial x}{\partial X} &=-\frac{1}{Z}(a_1f_x+a_3x)\tag{1}\\ \frac{\partial x}{\partial Y} &=-\frac{1}{Z}(b_1f_x+b_3x)\tag{2}\\ \frac{\partial x}{\partial Z} &=-\frac{1}{Z}(c_1f_x+c_3x)\tag{3}\\ \frac{\partial y}{\partial X} &=-\frac{1}{Z}(a_2f_y+a_3y)\tag{4}\\ \frac{\partial y}{\partial Y} &=-\frac{1}{Z}(b_2f_y+b_3y)\tag{5}\\ \frac{\partial y}{\partial Z} &=-\frac{1}{Z}(c_2f_y+c_3y)\tag{6} \end{align} ∂X∂x∂Y∂x∂Z∂x∂X∂y∂Y∂y∂Z∂y=−Z1(a1fx+a3x)=−Z1(b1fx+b3x)=−Z1(c1fx+c3x)=−Z1(a2fy+a3y)=−Z1(b2fy+b3y)=−Z1(c2fy+c3y)(1)(2)(3)(4)(5)(6)
旋转矩阵
[ a 1 a 2 a 3 b 1 b 2 b 3 c 1 c 2 c 3 ] = [ c o s ϕ cos ⁡ κ − s i n ϕ s i n ω c o s κ − c o s ϕ s i n κ − s i n ϕ s i n ω c o s κ − s i n ϕ c o s ω c o s ϕ s i n κ c o s ω c o s κ − s i n ω s i n ϕ cos ⁡ κ + c o s ϕ s i n ω s i n κ − s i n ϕ s i n κ + c o s ϕ s i n ω c o s κ c o s ϕ c o s ω ] \begin{equation} \begin{bmatrix} a_1 & a_2 & a_3 \\ b_1 & b_2 & b_3 \\ c_1 & c_2 & c_3 \end{bmatrix} = \begin{bmatrix} cos\phi\cos\kappa -sin\phi sin\omega cos\kappa & -cos\phi sin\kappa -sin\phi sin\omega cos\kappa & -sin\phi cos\omega \\ cos\phi sin\kappa & cos\omega cos\kappa & -sin\omega\\ sin\phi\cos\kappa +cos\phi sin\omega sin\kappa & -sin\phi sin\kappa +cos\phi sin\omega cos\kappa & cos\phi cos\omega\tag{1} \end{bmatrix} \end{equation} a1b1c1a2b2c2a3b3c3 = cosϕcosκ−sinϕsinωcosκcosϕsinκsinϕcosκ+cosϕsinωsinκ−cosϕsinκ−sinϕsinωcosκcosωcosκ−sinϕsinκ+cosϕsinωcosκ−sinϕcosω−sinωcosϕcosω (1)

三、代码展示

#include <Eigen/Dense>
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <vector>
#include<Eigen/Dense>
#include<math.h>
#include <iomanip> /*
1.为了便于矩阵运算，采用了Eigen库来进行矩阵的转置求逆等运算
*/using namespace std;
using namespace Eigen;const double f = 105.200;struct Img_Parameter//定义相机外方位元素
{int ImageID;double Xs;double Ys;double Zs;double Phi;double Omega;double Kappa;
};struct Point3D  //定义控制点地面坐标
{string PointName;double X;double Y;double Z;
};struct image //影像坐标
{int ID;double x;double y;
};//函数声明
void Read_phtData(string filename, vector<Img_Parameter>& imgPara_data);//读取相机外方位元素void Read_ptsData(string filename, vector<Point3D>& pts_data, vector<image>& img_data, vector<int>& ImagePointsNumber);//读取地面控制点X、Y、Zvoid FindID(const vector<Img_Parameter>& imgPara_data, const vector<image>& img_data,//将像点坐标所对应相机的外方位元素取出并整理成子集vector<Img_Parameter>& subImgPara_data);                                  //subImgPara_data方便进行条件平差void condition_adjustment(const vector<Point3D>& pts_data, const vector<image>& img_data,//进行条件平差const vector<Img_Parameter>& imgPara_data, const vector<int>& ImagePointsNumber,vector<Point3D>& my_data);void WriteToFile(vector<Point3D>& my_data, const string& filename);//写文件到“output.txt”int main()
{string file1 = "Data.pht";string file2 = "Data.pts";string outputfile = "output.txt";vector<Img_Parameter> imgPara_data;//将相机的内方位元素读取进img_datavector<Point3D> pts_data,my_data;//将控制点的物方坐标存储进pts_datavector<image> img_data;//读取影像坐标vector<int> ImagePointsNumber;//将每个点对应的影像数存储Read_phtData(file1, imgPara_data);Read_ptsData(file2, pts_data, img_data, ImagePointsNumber);condition_adjustment(pts_data, img_data, imgPara_data, ImagePointsNumber, my_data);WriteToFile(my_data, outputfile);return 0;
}void Read_phtData(string filename, vector<Img_Parameter>& imgPara_data)
{ifstream infile1(filename);if (!infile1){cout << "Cannot open file: " << filename << endl;}string line;//用getline读取每一行int data_count = 0;while (getline(infile1, line)) {// 忽略以 '$' 开头的注释行istringstream count(line);if (line[0] == '$') {continue;}else{count >> data_count;break;}}for (int i = 0; i < data_count; i++){getline(infile1, line);       //从读取到Image个数之后下一行开始读istringstream iss(line);Img_Parameter image;iss >> image.ImageID >> image.Xs >> image.Ys >> image.Zs>> image.Phi >> image.Omega >> image.Kappa;imgPara_data.push_back(image);//将其存入数组中}
}void Read_ptsData(string filename, vector<Point3D>& pts_data, vector<image>& img_data, vector<int>& ImagePointsNumber)
{ifstream file(filename);if (!file){cout << "Cannot open file." << endl;}string line;//用来读取每一行设置的字符串int Point_Number = 0;while (getline(file, line)){if (line[0] == '$')continue;else{istringstream iss(line);iss >> Point_Number;break;}}for (int i = 0; i < Point_Number; i++){getline(file, line);istringstream is(line);Point3D pts;is >> pts.PointName>> pts.X >> pts.Y >> pts.Z;pts_data.push_back(pts);int img_pts_Num;//设置照片点数getline(file, line);istringstream img_pts_num(line);img_pts_num >> img_pts_Num;ImagePointsNumber.push_back(img_pts_Num);//将每个点对应的影像数存储for (int j = 0; j < img_pts_Num; j++){image img;getline(file, line);istringstream is_image(line);is_image >> img.ID >> img.x >> img.y;img_data.push_back(img);}}}void condition_adjustment( const vector<Point3D>& pts_data, const vector<image>& img_data,const vector<Img_Parameter>& imgPara_data, const vector<int>& ImagePointsNumber, vector<Point3D>& my_data)
{int i = 0;int start = 0;for (auto it1 = pts_data.begin(); it1 != pts_data.end(); ++it1,i++)//遍历循环X, Y, Z{int o = 1;vector<Img_Parameter> subImgPara_data;//用来取出对应image的内方位元素Point3D mydata;MatrixXd X(3,1);mydata.PointName = it1->PointName;mydata.X = 720000;mydata.Y = 4340000;mydata.Z = 1000;int size = ImagePointsNumber[i];vector<image> subimage(img_data.begin() + start, img_data.begin() + start + size);// cout << "subimage " << i+1 << ":" << start<<" " << size<<endl;start += size;FindID(imgPara_data, subimage, subImgPara_data);//找到对应的外方位元素,这里的即为图片所对应的外方位元素do {double* _X = new double[size];//设置动态double* _Y = new double[size];double* _Z = new double[size];double* col_equa_x = new double[size];double* col_equa_y = new double[size];MatrixXd A(size * 2, 3);MatrixXd l(size * 2, 1);int j = 0;for (auto it2 = subImgPara_data.begin(); it2 != subImgPara_data.end(); ++it2, j++){double R[3][3] = { 0.0 }; //旋转矩阵RR[0][0] = cos(it2->Phi) * cos(it2->Kappa) - sin(it2->Phi) * sin(it2->Omega) * sin(it2->Kappa);R[0][1] = -cos(it2->Phi) * sin(it2->Kappa) - sin(it2->Phi) * sin(it2->Omega) * cos(it2->Kappa);R[0][2] = -sin(it2->Phi) * cos(it2->Omega);R[1][0] = cos(it2->Omega) * sin(it2->Kappa);R[1][1] = cos(it2->Omega) * cos(it2->Kappa);R[1][2] = -sin(it2->Omega);R[2][0] = sin(it2->Phi) * cos(it2->Kappa) + cos(it2->Phi) * sin(it2->Omega) * sin(it2->Kappa);R[2][1] = -sin(it2->Phi) * sin(it2->Kappa) + cos(it2->Phi) * sin(it2->Omega) * cos(it2->Kappa);R[2][2] = cos(it2->Phi) * cos(it2->Omega);_X[j] = R[0][0] * (mydata.X - it2->Xs) + R[1][0] * (mydata.Y - it2->Ys)+ R[2][0] * (mydata.Z - it2->Zs);_Y[j] = R[0][1] * (mydata.X - it2->Xs) + R[1][1] * (mydata.Y - it2->Ys)+ R[2][1] * (mydata.Z - it2->Zs);_Z[j] = R[0][2] * (mydata.X - it2->Xs) + R[1][2] * (mydata.Y - it2->Ys)+ R[2][2] * (mydata.Z - it2->Zs);col_equa_x[j] = -f * _X[j] / _Z[j];col_equa_y[j] = -f * _Y[j] / _Z[j];A(j * 2, 0) = -(1 / _Z[j]) * (R[0][0] * f + R[0][2] * col_equa_x[j]);A(j * 2, 1) = -(1 / _Z[j]) * (R[1][0] * f + R[1][2] * col_equa_x[j]);A(j * 2, 2) = -(1 / _Z[j]) * (R[2][0] * f + R[2][2] * col_equa_x[j]);A(j * 2 + 1, 0) = -(1 / _Z[j]) * (R[0][1] * f + R[0][2] * col_equa_y[j]);A(j * 2 + 1, 1) = -(1 / _Z[j]) * (R[1][1] * f + R[1][2] * col_equa_y[j]);A(j * 2 + 1, 2) = -(1 / _Z[j]) * (R[2][1] * f + R[2][2] * col_equa_y[j]);l(j * 2, 0) = subimage[j].x + f * _X[j] / _Z[j];l(j * 2 + 1, 0) = subimage[j].y + f * _Y[j] / _Z[j];}X = (A.transpose() * A).inverse() * A.transpose() * l;mydata.X += X(0,0);mydata.Y += X(1,0);mydata.Z += X(2,0);// cout << "间接平差系数B:" << A.format(Eigen::IOFormat(Eigen::StreamPrecision, Eigen::DontAlignCols, ", ", ",\n", "", "")) << endl;//cout << "第" << i + 1 << "组控制点坐标:" << endl;//  cout << fixed << setprecision(6) << X(0,0) <<" "<<X(1,0)<<" "<<X(2,0)<< endl;//cout << o++ << endl;delete[]_X;delete[]_Y;delete[]_Z;delete[]col_equa_x;delete[]col_equa_y;} while (abs(mydata.X - it1->X) > 0.01 || abs(mydata.Y - it1->Y) > 0.01 || abs(mydata.Z - it1->Z) > 0.01);/* cout << "控制点" << it1->PointName << "坐标:" << endl;cout <<"X: " << setprecision(9) << mydata.X << endl;cout <<"Y: " << setprecision(10) << mydata.Y << endl;cout <<"Z: " << setprecision(7) << mydata.Z << endl;cout << endl;*/my_data.push_back(mydata);subImgPara_data.clear();}
}void FindID(const vector<Img_Parameter>& imgPara_data, const vector<image>& img_data,vector<Img_Parameter>& subImgPara_data)
{for (int i = 0; i < img_data.size(); i++){for (auto it1 = imgPara_data.begin(); it1 != imgPara_data.end(); ++it1){if (img_data[i].ID == it1->ImageID){subImgPara_data.push_back(*it1);}}}}void WriteToFile(vector<Point3D>& my_data, const string& filename)
{ofstream outfile;outfile.open(filename, ios::out || ios::trunc);if (!outfile.is_open()) {cerr << "Error: failed to open the file " << filename << endl;return;}for (auto it = my_data.begin(); it != my_data.end(); ++it){       outfile << "控制点" << it->PointName << "坐标:" << endl;outfile << "X: " << setprecision(9) << it->X << endl;outfile << "Y: " << setprecision(10) << it->Y << endl;outfile << "Z: " << setprecision(7) << it->Z << endl;outfile << endl;}
}

四、小结

此算法很好的完成了所需任务，将与实际测量的误差控制在0.01m以内，达到了很高的精度。同时本算法也具有健壮性和完备性，用相同数据格式的数据集便可以达到相同的效果。
但是由于编程简便性，并没有将焦距等信息用函数读取，而是当作已知量直接进行计算。这也是以后值得完善的地方。
在编程的过程中，由于没有注意到各个变量的单位变化，导致循环迭代一直发散并不收敛，经过多次调试演算，才发现这一错误。同时也让我在接下来的学习过程中更加注意这些细节，从而完善自己。

空间前方交会（利用相机外方位元素和像点坐标进行解算）相关推荐

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

空间前方交会（利用相机外方位元素和像点坐标进行解算）

目录