单片空间后方交会（最详细原理、代码、结果分析）

大家好，欢迎来到IT知识分享网。

        使用数据：控制点物方空间坐标、控制点像平面量测坐标

        误差方程式：共线条件方程

        误差方程观测值：像点量测坐标

        误差方程未知数（改正值），本文我使用最全面的改正方程：外方位元素、内方位元素、畸变系数（可根据需求减少未知数，删除未知数改正值和对应的系数即可）

        程序流程：打开控制点物方坐标文件、控制点像点坐标文件读取数据，并根据合适的控制点，按照共线条件方程的空间后方交会误差方程形式，将像点坐标观测值作为观测值，将外方位元素改正值、内方位元素改正值、畸变系数改正值作为未知数，并设定合适的起算数据进行迭代解算。计算相片的内外方位元素和畸变系数（k1，k2，p1，p2）。并计算统计精度，保证中误差和残差在一定限度内。

        统计计算精度：单位权中误差、未知数中误差、像点观测值残差。

一、单片空间后方交会原理

        使用控制点物方坐标文件、控制点像点量测坐标文件，读入数据并保存在二维点、三维点数据结构中。然后将对应点的像平面坐标和物方坐标连接，保存成一条数据，便于后续使用。根据共线条件方程误差方程的空间后方交会形式，将像点量测值作为观测值，将内外方位元素、畸变系数（k1，k2，p1，p2）的改正值作为未知数，设定一个未知数的起算估计值进行迭代，求解最终的相片内外方位元素和畸变系数（k1，k2，p1，p2）。

      1.1误差方程构建

        定义旋转矩阵为：

        基于共线条件方程：

线性化后得到误差方程的一般形式，其中未知数为：

误差方程一般式为：

        （其中At为外方位元素、BcXc为控制点、BuXu为待定点、DadXad为畸变系数的系数及改正值，L为观测值误差项）

而后方交会的误差方程式将内外方位元素和畸变系数作为未知数：

        1.2误差方程系数填充

其中：

畸变模型选用：

各系数的定义如下：

        Dad即使用畸变模型对应未知数改正值前的系数即可。

系数求解中使用到的为像空间坐标系下的物方点坐标，故首先需要将控制点的物方坐标从物方坐标系转换到像空间坐标系S-xyz下：

观测值误差项，与使用畸变模型求解：

        1.3初始值（起算值）设置

        1.4逐点计算，合并系数矩阵

            将所有的像点都计算得到一个误差方程，计算每一个像点对应的误差方程未知数系数和误差项，并将所有的系数阵放在一起求解，系数合并为大系数阵A，大小为（2*num_of_points）*13（13为未知数/改正数的个数，结合实际求解未知数设置），观测值误差阵合并为L，L矩阵大小为（2*num_of_points）*1。

      1.5迭代计算求最优解

        最后利用最小二乘平差原理，计算未知数改正值：

        将所有的未知数加上计算得到的改正数，得到新的值，再次代回第一步迭代计算，直到每一项的改正值小于某一个阈值为止（我设置的是0.0001）。输出最后得到的未知数平差计算结果。

        1.6精度统计

二、单片空间后方交会代码

#include <opencv2/highgui.hpp> #include <opencv2/core.hpp> //#include <iostream> #include <math.h> #include<fstream> #include<sstream> #include<vector> #include<iostream> #include<iomanip> using namespace std; using namespace cv; // 控制场平面坐标系 // ^ X // | // | // |——————> Y(左手系和右手系的转换） #define Pi 3. #define ROWS 4160 #define COLS 6240 #define PIXSIZE 0.00376 //--------------------------结构体申明-------------------------------------- struct ControlPoint { int flag; double X; double Y; double Z; }; struct imgPoint { int flag; double x; double y; }; struct merge_point { int flag; double X; double Y; double Z; double x; double y; }; //方位元素和畸变系数 struct orienParam { //内方位元素 double x0 = 0; double y0 = 0; double f = 0; //畸变系数 double k1 = 0; double k2 = 0; double p1 = 0; double p2 = 0; //外方位元素 double Phi = 0; double Omega = 0; double Kappa = 0; double Xs = 0; double Ys = 0; double Zs = 0; }; //------------------坐标文件读取与坐标系变换--------------------- // 标志点量测坐标定义 标志点像平面坐标定义 // --------> x(pixel) ^ y(mm) -(y-COLS/2) // | | // | | // | -----------------> x (x-ROWS/2) // | | // y | // void Pixel2Img_with_co_change(vector<imgPoint>& imgPoints) { for (int i = 0; i < imgPoints.size(); i++) { imgPoints[i].x = (imgPoints[i].x - COLS / 2) * PIXSIZE; imgPoints[i].y = -1.0 * (imgPoints[i].y - ROWS / 2) * PIXSIZE; } } void readctrl_ptData(char* file, vector<ControlPoint>& ControlPoints) { ifstream inFile(file); if (!inFile) { cout << "文件读取失败，请检查文件路径" << endl; return; } string line; string firstLine; getline(inFile, firstLine); while (getline(inFile, line)) { ControlPoint ctrl_pt; istringstream iss(line); iss >> ctrl_pt.flag; //坐标系转换 iss >> ctrl_pt.Z; iss >> ctrl_pt.X; iss >> ctrl_pt.Y; ctrl_pt.Z = -1.0 * ctrl_pt.Z; /* Z Y | X | | / => |______>X |/ / ————————>Y /Z */ ControlPoints.push_back(ctrl_pt); } inFile.close(); } void readImgData(char* file, vector<imgPoint>& imgPoints) { ifstream inFile(file); if (!inFile) { cout << "文件读取失败，请检查文件路径" << endl; return; } string line; while (getline(inFile, line)) { imgPoint img_pt; istringstream iss(line); iss >> img_pt.flag; iss >> img_pt.x; iss >> img_pt.y; imgPoints.push_back(img_pt); } inFile.close(); } void Merge_Point(vector<imgPoint>& imgPoints, vector<ControlPoint>& ControlPoints, vector<merge_point>& merge_points) { merge_point pair; for (int i = 0; i < imgPoints.size(); i++) { for (int j = 0; j < ControlPoints.size(); j++) { if (imgPoints[i].flag == ControlPoints[j].flag) { pair.flag = imgPoints[i].flag; pair.X = ControlPoints[j].X; pair.Y = ControlPoints[j].Y; pair.Z = ControlPoints[j].Z; pair.x = imgPoints[i].x; pair.y = imgPoints[i].y; merge_points.push_back(pair); break; } } } } //-----------------误差方程建立------------------------ void cal_Coefficient(Mat& A, Mat& l, vector<merge_point>& merge_points, orienParam& orien) { double x0 = orien.x0; double y0 = orien.y0; double f = orien.f; double k1 = orien.k1; double k2 = orien.k2; double p1 = orien.p1; double p2 = orien.p2; double Phi = orien.Phi; double Omega = orien.Omega; double Kappa = orien.Kappa; double Xs = orien.Xs; double Ys = orien.Ys; double Zs = orien.Zs; //计算旋转矩阵 Mat R = Mat::zeros(3, 3, CV_64FC1); R.at<double>(0, 0) = cos(Phi) * cos(Kappa) - sin(Phi) * sin(Omega) * sin(Kappa); R.at<double>(0, 1) = -1.0 * cos(Phi) * sin(Kappa) - sin(Phi) * sin(Omega) * cos(Kappa); R.at<double>(0, 2) = -1.0 * sin(Phi) * cos(Omega); R.at<double>(1, 0) = cos(Omega) * sin(Kappa); R.at<double>(1, 1) = cos(Omega) * cos(Kappa); R.at<double>(1, 2) = -1.0 * sin(Omega); R.at<double>(2, 0) = sin(Phi) * cos(Kappa) + cos(Phi) * sin(Omega) * sin(Kappa); R.at<double>(2, 1) = -1.0 * sin(Phi) * sin(Kappa) + cos(Phi) * sin(Omega) * cos(Kappa); R.at<double>(2, 2) = cos(Phi) * cos(Omega); for (int i = 0; i < merge_points.size(); i++) { double X = merge_points[i].X; double Y = merge_points[i].Y; double Z = merge_points[i].Z; double x = merge_points[i].x; double y = merge_points[i].y; double r_2 = pow(x - x0, 2) + pow(y - y0, 2); double delta_x = (x - x0) * (k1 * r_2 + k2 * r_2 * r_2) + p1 * (r_2 + pow((x - x0), 2)) + 2 * p2 * (x - x0) * (y - y0); double delta_y = (y - y0) * (k1 * r_2 + k2 * r_2 * r_2) + p2 * (r_2 + pow((y - y0), 2)) + 2 * p1 * (x - x0) * (y - y0); double Xbar = R.at<double>(0, 0) * (X - Xs) + R.at<double>(1, 0) * (Y - Ys) + R.at<double>(2, 0) * (Z - Zs); double Ybar = R.at<double>(0, 1) * (X - Xs) + R.at<double>(1, 1) * (Y - Ys) + R.at<double>(2, 1) * (Z - Zs); double Zbar = R.at<double>(0, 2) * (X - Xs) + R.at<double>(1, 2) * (Y - Ys) + R.at<double>(2, 2) * (Z - Zs); //外方位元素系数->线元素 double a11 = (R.at<double>(0, 0) * f + R.at<double>(0, 2) * (x - x0 + delta_x)) / Zbar; double a12 = (R.at<double>(1, 0) * f + R.at<double>(1, 2) * (x - x0 + delta_x)) / Zbar; double a13 = (R.at<double>(2, 0) * f + R.at<double>(2, 2) * (x - x0 + delta_x)) / Zbar; double a21 = (R.at<double>(0, 1) * f + R.at<double>(0, 2) * (y - y0 + delta_y)) / Zbar; double a22 = (R.at<double>(1, 1) * f + R.at<double>(1, 2) * (y - y0 + delta_y)) / Zbar; double a23 = (R.at<double>(2, 1) * f + R.at<double>(2, 2) * (y - y0 + delta_y)) / Zbar; //外方位元素系数->角元素 double a14 = (y - y0 + delta_y) * sin(Omega) - (((x - x0 + delta_x) / f) * ((x - x0 + delta_x) * cos(Kappa) - (y - y0 + delta_y) * sin(Kappa)) + f * cos(Kappa)) * cos(Omega); double a15 = -f * sin(Kappa) - ((x - x0 + delta_x) / f) * ((x - x0 + delta_x) * sin(Kappa) + (y - y0 + delta_y) * cos(Kappa)); double a16 = y - y0 + delta_y; double a24 = -1.0 * (x - x0 + delta_x) * sin(Omega) - (((y - y0 + delta_y) / f) * ((x - x0 + delta_x) * cos(Kappa) - (y - y0 + delta_y) * sin(Kappa)) - f * sin(Kappa)) * cos(Omega); double a25 = -f * cos(Kappa) - ((y - y0 + delta_y) / f) * ((x - x0 + delta_x) * sin(Kappa) + (y - y0 + delta_y) * cos(Kappa)); double a26 = -1.0 * (x - x0 + delta_x); //内方位元素系数 double a17 = (x - x0 + delta_x) / f; double a18 = 1; double a19 = 0; double a27 = (y - y0 + delta_y) / f; double a28 = 0; double a29 = 1; //畸变系数 double a110 = -1.0 * (x - x0 + delta_x) * r_2; //k1 double a111 = -1.0 * (x - x0 + delta_x) * r_2 * r_2; //k2 double a112 = -1.0 * (2 * pow((x - x0 + delta_x), 2) + r_2); //p1 double a113 = -2.0 * (y - y0 + delta_y) * (x - x0 + delta_x); //p2 double a210 = -1.0 * (y - y0 + delta_y) * r_2; //k1 double a211 = -1.0 * (y - y0 + delta_y) * r_2 * r_2; //k2 double a212 = -2.0 * (y - y0 + delta_y) * (x - x0 + delta_x); //p1 double a213 = -1.0 * (2 * pow((y - y0 + delta_y), 2) + r_2); //p2 //组合为系数阵 A.at<double>(2 * i, 0) = a11; A.at<double>(2 * i + 1, 0) = a21; A.at<double>(2 * i, 1) = a12; A.at<double>(2 * i + 1, 1) = a22; A.at<double>(2 * i, 2) = a13; A.at<double>(2 * i + 1, 2) = a23; A.at<double>(2 * i, 3) = a14; A.at<double>(2 * i + 1, 3) = a24; A.at<double>(2 * i, 4) = a15; A.at<double>(2 * i + 1, 4) = a25; A.at<double>(2 * i, 5) = a16; A.at<double>(2 * i + 1, 5) = a26; A.at<double>(2 * i, 6) = a17; A.at<double>(2 * i + 1, 6) = a27; A.at<double>(2 * i, 7) = a18; A.at<double>(2 * i + 1, 7) = a28; A.at<double>(2 * i, 8) = a19; A.at<double>(2 * i + 1, 8) = a29; A.at<double>(2 * i, 9) = a110; A.at<double>(2 * i + 1, 9) = a210; A.at<double>(2 * i, 10) = a111; A.at<double>(2 * i + 1, 10) = a211; A.at<double>(2 * i, 11) = a112; A.at<double>(2 * i + 1, 11) = a212; A.at<double>(2 * i, 12) = a113; A.at<double>(2 * i + 1, 12) = a213; //观测值 l.at<double>(2 * i, 0) = x - (x0 - f * Xbar / Zbar - delta_x); l.at<double>(2 * i + 1, 0) = y - (y0 - f * Ybar / Zbar - delta_y); } } //---------------------误差方程迭代求解-------------------- void calcu_left(char ctrl_path[],char left_path[]) { vector<ControlPoint> ControlPoints; vector<imgPoint> left_imgPoints; readctrl_ptData(ctrl_path, ControlPoints); readImgData(left_path, left_imgPoints); //定义两张相片的内外方位元素 orienParam left_orien; left_orien.x0 = 0; left_orien.y0 = 0; left_orien.f = 27; left_orien.Xs = 1500; left_orien.Ys = -200; left_orien.Zs = -1000; left_orien.Kappa = 0; left_orien.Omega = 0; left_orien.Phi = 0.3; left_orien.k1 = 0; left_orien.k2 = 0; left_orien.p1 = 0; left_orien.p2 = 0; //将像素坐标(pixel)转换为图像坐标(mm) Pixel2Img_with_co_change(left_imgPoints); //生成点对 vector<merge_point> left_merge_points; Merge_Point(left_imgPoints, ControlPoints, left_merge_points); //构建左片系数阵A和常数项l Mat A(2 * left_merge_points.size(), 13, CV_64FC1); Mat l(2 * left_merge_points.size(), 1, CV_64FC1); Mat X(13, 1, CV_64FC1); int i = 0; while (true) { i++; printf("%d", i); //填充系数阵 cal_Coefficient(A, l, left_merge_points, left_orien); //cout << "L = " << l << endl; //最小二乘 X = (A.t() * A).inv() * A.t() * l; //更新外方位元素 left_orien.Xs += X.at<double>(0, 0); left_orien.Ys += X.at<double>(1, 0); left_orien.Zs += X.at<double>(2, 0); left_orien.Phi += X.at<double>(3, 0); left_orien.Omega += X.at<double>(4, 0); left_orien.Kappa += X.at<double>(5, 0); left_orien.f += X.at<double>(6, 0); left_orien.x0 += X.at<double>(7, 0); left_orien.y0 += X.at<double>(8, 0); left_orien.k1 += X.at<double>(9, 0); left_orien.k2 += X.at<double>(10, 0); left_orien.p1 += X.at<double>(11, 0); left_orien.p2 += X.at<double>(12, 0); if (abs(X.at<double>(0, 0)) < 0.0001 && abs(X.at<double>(1, 0)) < 0.0001 && abs(X.at<double>(2, 0)) < 0.0001 && abs(X.at<double>(3, 0)) < 0.0001 && abs(X.at<double>(4, 0)) < 0.0001 && abs(X.at<double>(5, 0)) < 0.0001 && abs(X.at<double>(6, 0)) < 0.0001 && abs(X.at<double>(7, 0)) < 0.0001 && abs(X.at<double>(8, 0)) < 0.0001 && abs(X.at<double>(9, 0)) < 0.0001 && abs(X.at<double>(10, 0)) < 0.0001 && abs(X.at<double>(11, 0)) < 0.0001 && abs(X.at<double>(12, 0)) < 0.0001) break; } //输出定位参数left_orien cout << "内外方位元素和畸变系数：" << endl; cout << "Xs: " << left_orien.Xs << endl; cout << "Ys: " << left_orien.Ys << endl; cout << "Zs: " << left_orien.Zs << endl; cout << "Phi: " << left_orien.Phi << endl; cout << "Omega: " << left_orien.Omega << endl; cout << "Kappa: " << left_orien.Kappa << endl; cout << "f: " << left_orien.f << endl; cout << "x0: " << left_orien.x0 << endl; cout << "y0: " << left_orien.y0 << endl; cout << "k1: " << left_orien.k1 << endl; cout << "k2: " << left_orien.k2 << endl; cout << "p1: " << left_orien.p1 << endl; cout << "p2: " << left_orien.p2 << endl; cout << "----------------------------------" << endl; Mat V = A * X - l; Mat vtv = V.t() * V; double sigma0 = sqrt(vtv.at<double>(0, 0) / (2 * left_merge_points.size() - 13));//单位权中误差 cout << "单位权中误差: " << endl << sigma0 << "mm" << " or " << sigma0 / PIXSIZE << "pixel" << endl; cout << "----------------------------------" << endl; //精度统计 --> 各未知数的中误差 Mat Qxx = (A.t() * A).inv(); vector<string> str = { "Xs", "Ys", "Zs", "Phi", "Omega", "Kappa", "f", "x0", "y0", "k1", "k2", "p1", "p2" }; cout << "未知数中误差：" << endl; for (int i = 0; i < 13; i++) { double sigmaXi = sqrt(Qxx.at<double>(i, i)) * sigma0; cout << str[i] << ": " << sigmaXi << endl; } cout << "----------------------------------" << endl; //精度统计 --> 单位权中误差 cout << "像点观测值残差/pixel" << endl; cout << "id" << " vx " << " vy" << endl; for (int i = 0; i < left_merge_points.size(); i++) { cout << left_merge_points[i].flag << " " << setiosflags(ios::right) << fixed << setprecision(5) << V.at<double>(2 * i, 0) / PIXSIZE << " " << setiosflags(ios::right) << fixed << setprecision(5) << V.at<double>(2 * i + 1, 0) / PIXSIZE << endl; } } void calcu_right(char ctrl_path[], char right_path[]) { vector<ControlPoint> ControlPoints; vector<imgPoint> right_imgPoints; readctrl_ptData(ctrl_path, ControlPoints); readImgData(right_path, right_imgPoints); //定义两张相片的内外方位元素 orienParam right_orien; right_orien.x0 = 0; right_orien.y0 = 0; right_orien.f = 27; right_orien.Xs = 4500; right_orien.Ys = 200; right_orien.Zs = -1000; right_orien.Kappa = 0; right_orien.Omega = 0; right_orien.Phi = -0.2678; right_orien.k1 = 0; right_orien.k2 = 0; right_orien.p1 = 0; right_orien.p2 = 0; //将像素坐标(pixel)转换为图像坐标(mm) Pixel2Img_with_co_change(right_imgPoints); //生成点对 vector<merge_point> right_merge_points; Merge_Point(right_imgPoints, ControlPoints, right_merge_points); //构建左片系数阵A和常数项l Mat A(2 * right_merge_points.size(), 13, CV_64FC1); Mat l(2 * right_merge_points.size(), 1, CV_64FC1); Mat X(13, 1, CV_64FC1); int i = 0; while (true) { i++; printf("%d", i); //填充系数阵 cal_Coefficient(A, l, right_merge_points, right_orien); //cout << "L = " << l << endl; //最小二乘 X = (A.t() * A).inv() * A.t() * l; //更新外方位元素 right_orien.Xs += X.at<double>(0, 0); right_orien.Ys += X.at<double>(1, 0); right_orien.Zs += X.at<double>(2, 0); right_orien.Phi += X.at<double>(3, 0); right_orien.Omega += X.at<double>(4, 0); right_orien.Kappa += X.at<double>(5, 0); right_orien.f += X.at<double>(6, 0); right_orien.x0 += X.at<double>(7, 0); right_orien.y0 += X.at<double>(8, 0); right_orien.k1 += X.at<double>(9, 0); right_orien.k2 += X.at<double>(10, 0); right_orien.p1 += X.at<double>(11, 0); right_orien.p2 += X.at<double>(12, 0); if (abs(X.at<double>(0, 0)) < 0.0001 && abs(X.at<double>(1, 0)) < 0.0001 && abs(X.at<double>(2, 0)) < 0.0001 && abs(X.at<double>(3, 0)) < 0.0001 && abs(X.at<double>(4, 0)) < 0.0001 && abs(X.at<double>(5, 0)) < 0.0001 && abs(X.at<double>(6, 0)) < 0.0001 && abs(X.at<double>(7, 0)) < 0.0001 && abs(X.at<double>(8, 0)) < 0.0001 && abs(X.at<double>(9, 0)) < 0.0001 && abs(X.at<double>(10, 0)) < 0.0001 && abs(X.at<double>(11, 0)) < 0.0001 && abs(X.at<double>(12, 0)) < 0.0001) break; } //输出定位参数right_orien cout << "内外方位元素和畸变系数：" << endl; cout << "Xs: " << right_orien.Xs << endl; cout << "Ys: " << right_orien.Ys << endl; cout << "Zs: " << right_orien.Zs << endl; cout << "Phi: " << right_orien.Phi << endl; cout << "Omega: " << right_orien.Omega << endl; cout << "Kappa: " << right_orien.Kappa << endl; cout << "f: " << right_orien.f << endl; cout << "x0: " << right_orien.x0 << endl; cout << "y0: " << right_orien.y0 << endl; cout << "k1: " << right_orien.k1 << endl; cout << "k2: " << right_orien.k2 << endl; cout << "p1: " << right_orien.p1 << endl; cout << "p2: " << right_orien.p2 << endl; cout << "----------------------------------" << endl; Mat V = A * X - l; Mat vtv = V.t() * V; double sigma0 = sqrt(vtv.at<double>(0, 0) / (2 * right_merge_points.size() - 13));//单位权中误差 cout << "单位权中误差: " << endl << sigma0 << "mm" << " or " << sigma0 / PIXSIZE << "pixel" << endl; cout << "----------------------------------" << endl; //精度统计 --> 各未知数的中误差 Mat Qxx = (A.t() * A).inv(); vector<string> str = { "Xs", "Ys", "Zs", "Phi", "Omega", "Kappa", "f", "x0", "y0", "k1", "k2", "p1", "p2" }; cout << "未知数中误差：" << endl; for (int i = 0; i < 13; i++) { double sigmaXi = sqrt(Qxx.at<double>(i, i)) * sigma0; cout << str[i] << ": " << sigmaXi << endl; } cout << "----------------------------------" << endl; //精度统计 --> 单位权中误差 cout << "像点观测值残差/pixel" << endl; cout << "id" << " vx " << " vy" << endl; for (int i = 0; i < right_merge_points.size(); i++) { cout << right_merge_points[i].flag << " " << setiosflags(ios::right) << fixed << setprecision(5) << V.at<double>(2 * i, 0) / PIXSIZE << " " << setiosflags(ios::right) << fixed << setprecision(5) << V.at<double>(2 * i + 1, 0) / PIXSIZE << endl; } } //------------------------主函数------------------------ int main() { //读取控制点数据和左右片标志点的像素坐标 char ctrl_ptfile[] = "C:/Users/LENOVO/Desktop/近景摄影测量/2024实习_近景控制场_.txt"; //char left_file[] = "C:/Users/LENOVO/Desktop/近景摄影测量/coordinates.txt"; //char right_file[] = "C:/Users/LENOVO/Desktop/近景摄影测量/coordinates_right.txt"; char left_file[] = "C:/Users/LENOVO/Desktop/近景摄影测量/datapt/points_imgL.txt"; char right_file[] = "C:/Users/LENOVO/Desktop/近景摄影测量/datapt/points_imgR.txt"; printf("-----------左片-----------------\n"); calcu_left(ctrl_ptfile,left_file); printf("-----------右片-----------------\n"); calcu_right(ctrl_ptfile, right_file); return 0; } // 运行程序: Ctrl + F5 或调试 >“开始执行(不调试)”菜单 // 调试程序: F5 或调试 >“开始调试”菜单 // 入门使用技巧: // 1. 使用解决方案资源管理器窗口添加/管理文件 // 2. 使用团队资源管理器窗口连接到源代码管理 // 3. 使用输出窗口查看生成输出和其他消息 // 4. 使用错误列表窗口查看错误 // 5. 转到“项目”>“添加新项”以创建新的代码文件，或转到“项目”>“添加现有项”以将现有代码文件添加到项目 // 6. 将来，若要再次打开此项目，请转到“文件”>“打开”>“项目”并选择 .sln 文件 

三、单片空间后方交会结果与分析

        在这一步我设立了两组数据运行计算，分别利用81(左）/63（右）对控制点和9（左）/19（右）对控制点进行控制计算。

        3.1未知数求解结果

表1 未知数求解结果

        3.2精度统计-中误差

表2 中误差

        3.3精度统计-像点量测残差

        第一组81左片/63右片结果（单位/pixel）：

表3 像点量测残差

        第二组9左片/19右片结果（单位/pixel）：

表4 像点量测残差

四、误差剔除

        观察表1未知数求解结果，可以发现，使用较少的控制点和较多的控制点得到的结果相差较大，外方位元素可相差数十毫米，但整体解求结果无明显错误。

        观察表2可知，利用较少的控制点得到的结果误差较明显，例如Zs中误差，使用较少控制点得到的中误差达到了23.6963mm，而使用较多控制点得到的中误差减小到了2.6982mm。使用较多控制点计算得到的结果虽然误差明显下降且误差较小，外方位元素中误差可达到2~3mm以内，像点量测单位权中误差控制在了亚像素级别，但是像点量测单位权中误差却比使用较少控制点得到的误差大，我认为是较多的控制点中有一些量测误差过大的点。

        观察表3、4可以看出，两组实验得到的像点观测值残差都能基本保持在一个像素以内，而且正如上述观点，使用较多的控制点计算时，有许多像点量测误差过大的点，存在明显的量测错误，需要剔除或修改，否则影响整体精度。例如左片的：432、332、335、373、143号点，右片的：335、135、136、133号点，他们的量测误差达到了2~3个像元，尤其是335号点达到了数十个像元。

        我通过剔除与重新量测这些点来解决这个问题,剔除这些点后，单位权中误差减小到了和使用较少控制点相近的误差值，且所有未知数的误差值均有所下降，达到较为理想的结果。

表5 剔除粗差前后中误差

        最终得到的未知数结算结果为（其他未知数变化微小，在上述表中已经展示过，此处不再重复，仅展示外方位线元素）：

表6 剔除误差前后未知数结果（单位/mm）

五、误差来源分析

        经过误差剔除与重新测量（手动平移对准点）后的结果达到了较为理想的值，我将粗差较大的点记录下来，返回到像点量测程序中去再次框选查看，发现这些粗差点有一些规律：

• 位于影像边缘或倾斜程度较大，呈现椭圆形或变形。
• 标志点受光照射，有反光。
• 像点量测程序定位点存在明显偏差，需要手动调整

        所以像点量测程序其实也有待优化，可以使用其他的算法对抗这些存在的问题。

六、误差影响分析

        若手动加入一个很大的粗差，究竟会对整体求解造成哪些影响？我将后方交会程序中的左片的433号点的x坐标从41.38改到1000，手动加入了一个非常大的粗差，运行代码，发现会造成以下问题：

1. 迭代次数变多：从18次增加到了47次收敛
2. 未知数计算值明显变化：

表7 加入粗差前后未知数值

     3.统计误差变大

表8 加入粗差前后中误差

1. 像点量测残差变大，且带动周围点残差变大

表9 加入粗差前后像点量测残差

七、参考文献

冯文灏.近景摄影测量[M].武汉：武汉大学出版社，2002，116-159

欢迎指正！祝学业有成、学习顺利、开心每一天