数字摄影测量实习报告

1、.课程编号： 课程性质：必修数字摄影测量实习学院： 测绘学院 专业： 测绘工程 地点： 校 内 班级： 201116110 组号： 姓名： 学号： 教师： 徐 芳 李 妍 申丽丽 2014年4月10日 至 2014年4月30日*;目录一、 实习目的1二、 VirtuoZo实习过程11、 数据准备以及初步处理11） 新建测区12） 初步处理13） 进一步处理24） 产品生成22、 DEM编辑23、 DEM拼接34、 图像镶嵌45、 IGS数字化测图46、 空中三角测量51） 数据准备及预处理52） 交互式编辑63） 三角测量74） DEM拼接和正射影像制作7三、 GeoWay实习过程81、 交互

2、式编辑81） 新建工程并设置相关参数82） 生成内图廓.83） 交互式编辑92、 空三加密101） 数据准备及检查102） 控制网平差113） 精确匹配并查看匹配结果114） 刺点113、 立体测图121） 新建方案和数据准备122） 初步处理123） 批处理134） 地物采集135） DEM编辑136） TIN网编辑147） 质量检查15四、 编程171、 空间后方交会172、 特征匹配213、 同名核线提取25五、 心得体会28一、 实习目的通过对数字摄影测量软件的操作，熟悉数字摄影测量系统，掌握数字摄影测量软件的作业流程，加深对数字摄影测量过程的理解。通过编程实习，加深对数字摄影测量学相

3、关方法原理的理解，提高编程动手能力。二、 VirtuoZo实习过程1、数据准备以及初步处理1) 新建测区à设置测区参数，新建模型à设置模型参数，设置相机参数à设置控制点数据à导入影像数据，至此，准备工作完成； 图1 设置测区参数 图2 设置模型参数 图3 导入控制点 图4 导入相机参数图5 导入影像2) 内定向à相对定向à绝对定向，至此影像定向完成；图6 内定向图7 相对定向相对定向做完之后系统会自动进行绝对定向，我们可以到主界面查看绝对定向的结果：下图是模型157-156的绝对定向结果图：图8 157-156绝对定向结果3) 然后点

4、击处理à核线重采样à影像匹配，做完这些就可以准备生成产品；4) 然后点击产品，可以生成需要的产品，如DEM，DOM，等高线等；2、DEM编辑由于生成的DEM含有房屋，树木等地物，所以并不是真正意义上的数字高程模型，而应该是DTM（数字表面模型）。如下图：图9 数字表面模型图中的红色方框所包围的可以很容易看出来是房屋而不是地形。所以需要对上上述DEM进行编辑，生成正真意义上的DEM。点击“高级àDEM编辑”，进入DEM编辑界面。图10 编辑后的DEM图中的红色部分就是我们编辑的范围，经过编辑后加载DEM显示立体图，可以看出其中的房屋高度已经贴合地面了。3、DEM拼接

5、由于每一个DEM只是两张影像形成的，一测区有多个DEM，因此需要经过拼接才能形成整个测区的DEM。点击“高级->DEM拼接”弹出DEM拼接界面，加载需要拼接的DEM，用鼠标划定拼接范围，接下来就可以执行拼接。 图11图12 DEM拼接4、图像镶嵌生成整个测区的DEM后，我们就可以生成整个测区的正射影像了。依次点击“高级à正射影像制作”弹出制作界面，选定DEM以及相应的影像数据，然后执行，结果如图图13 拼接后的正射影像5、IGS数字化测图点击“测图àIGS数字化测图”进入到数字化测图界面，然后依次“新建VZV矢量文件à装载影像数据”，然后通过红绿眼睛对测区进

6、行矢量化。图14 IGS数字化测图图15 矢量图形6、空中三角测量1) 数据准备和图像预处理打开VirtuozoAAT，依次“新建测区设置参数-引入影像”图16 输入影像注意第二条航带的选项要选上相机旋转。然后进行“内定向à航带间偏移量”图17 航带间偏移量在图中选取航带间的同名点，系统会自动计算航带间偏移量设置好航带间偏移量之后，开始提取连接点，最后运行自动挑点程序。2) 交互式编辑打开交互式编辑界面，一次在每张图像上选点，保证每个蓝框周围至少有3个点就行。 图18 选点 图19 使用PATB平差 然后加入控制点，只要在测区中加入三个控制点，然后使用PATB进行平差，就可以使用“预

7、测控制点”预测出其余的控制点的位置，然后依次选出控制点坐标。 图20 控制点 图21 被挑出的点可能一次选择的点精度不够会有误差超限的点，显示在如图21中带有星号，然后对带有星号的点从新定位，这样直到没有星号的点。至此，交互式编辑完成。3) 三角测量打开VirtuozoNT，打开在AAT中建立的测区，自动创建测区，之就可以使用批处理功能进行自动化处理。至此，三角测量完成。图22 批处理界面4) DEM拼接和正射影像制作按照前面的拼接步骤进行DEM拼接和正射影像的制作，结果如下： 图23 DEM 图24 正射影像3、 GeoWay实习流程1、交互式编辑1）新建工程，设置相关参数（PS：对照表添加

8、 LAYER ENDSEC，然后将SYM符号库文件拷贝到C:WINDOWSGeoSymbol） 图1 参数设置界面2）生成内图廓线 图2 生成图廓线3）交互式地物编辑1、加载栅格图层。点击栅格图层，加载tfw数据，加入栅格图层，然后生成内图廓，格网线 图3 加载栅格图层2、勾选跨层输入编辑，设置激活图层，然后点击相应地物à设为输入地物类。3、点、单线、双线编辑4、点、线处理。悬挂点和伪节点处理，工具检查悬挂伪节点进行伪节点检测。批处理消除悬挂点和伪节点能够消除大部分错误，不能消除的部分进行手工编辑。 图4 节点编辑后的图5、面处理。新建房屋对照表，在拓扑构建时启用房屋对照表，然后拓扑

9、弧段预处理，最后进行拓扑构建。对于拓扑构建中被遗漏的要手工处理后再次构建图5 拓扑构建后的图2、GEOWAY空三加密1）新建工程，数据准备并检查依次点击“新建工程à选择数据所在的目录初始构网”然后点击“生成成果à生成单色点云”用mashlab打开.ply文件，查看点云和构网结果 图6 点云图 图7 网形图2）控制网平差“常规处理à光束法平差à删除残差大于域值的点à自由网平差”，重复上述步骤直到最小残差小于阈值60（PS：每次删除的点不能太多，以防删除过多正确的点）3）精确匹配并查看匹配结果点击“常规处理à精确匹配”，进一步匹配图8 图

10、94）刺点常规处理光束法平差à常规处理à绝对定向（GPS）à刺点图10 刺点界面5）空三结果图113、立体测图（无人机）1）新建方案和数据准备影像管理à加载影模型管理à新建模型控制点编辑à导入控制点像点导入à导入像点外方为元素导入à导入外方位元素2） 初步处理内定向相对定向绝对定向核线重采样图12 相对定向3）批处理图134） 地物采集新建加工工程，然后依次进行高程点、道路、房屋、比高点、等高线的采集图14 地物采集结果5）DEM编辑1、点线处理 采集特征线处理选择特征线处理2、面处理自适应内插，量测点内插局部置平

11、、水域处理 图15 DEM编辑前 图16 DEM编辑后6）TIN网处理图17 编辑前 图18 编辑后图19 编辑后 7）质量检查1、设置质检方案制作方案加入工作方案 图20 检查报告 图21 方案设置 2、新建检查任务，指定目录到数据所在目录。然后自动检查错误并查看错误，然后对矢量图进行修改图22 编辑前 图23 编辑后 图24 编辑前 图25 编辑后四、编程内容1、单像空间后方交会#include "hfjhdialog.h"#include <opencv.hpp>#include <vector>#include <string>#

12、include <fstream>#include <iostream>#include <QFileDialog>using namespace std;HFJHDialog:HFJHDialog(QWidget *parent) : QDialog(parent), ui(new Ui:HFJHDialog) ui->setupUi(this);char *title= "空间后方交会"QString qtitle = QString:fromLocal8Bit(title);setWindowTitle(qtitle);ui-

13、>textBrowser->setText("-HFJH-"); /焦距 单位pixel fx=4547.99665; fy=4547.87373; /内方位元素 pixel x0=47.48571; y0=12.02756; /畸变参数 k1=-5.00793e-009; k2=1.83462e-016; p1=-2.24419e-008; p2=1.76820e-008; /初始值 Xs0=4300.0; Ys0=-150.0; Zs0=-300.0;/单位mm phi0=0.0; omega0=0.0; kappa0=0.0;/单位radHFJHDialo

14、g:HFJHDialog() delete ui;/读取像点坐标int HFJHDialog:readxy(string path) int index = 1; /path ="F:我的课程数字摄影测量影像数据后方交会像点坐标.txt" ifstream fil; fil.open(path); if(!fil) char *log = "打开像点文件失败<<<"ui->textBrowser->append(QString:fromLocal8Bit(log); return -1; xy tmp; while(!fil

15、.eof() fil>>tmp.index>>tmp.x>>tmp.y;ui->textBrowser->append(QString("Point%1: x:%2 y:%3").arg(index).arg(tmp.x).arg(tmp.y);index+; pic.push_back(tmp); fil.close(); return 0;/读取物方点坐标int HFJHDialog:readXYZ(string path) int index = 1; ifstream fil; fil.open(path); if(!

16、fil) char *log = "打开控制点文件失败<<<"ui->textBrowser->append(QString:fromLocal8Bit(log); return -1; XYZ tmp; string head; for(int i = 0; i < 4; +i) fil>>head; while(!fil.eof() fil>>tmp.index>>tmp.X>>tmp.Y>>tmp.Z;ui->textBrowser->append(QStrin

17、g("Point%1: X:%2 Y:%3 Z:%3").arg(index).arg(tmp.X).arg(tmp.Y).arg(tmp.Z); index+; obj.push_back(tmp); fil.close(); return 0;/读取相机参数文件int HFJHDialog:readcamera(string path)int index = 1; /path ="F:我的课程数字摄影测量影像数据后方交会像点坐标.txt" ifstream fil; fil.open(path); if(!fil) char *log = "

18、;打开相机文件失败<<<"ui->textBrowser->append(QString:fromLocal8Bit(log); return -1; fil.close(); return 0;/空间后方交会算法主体函数void HFJHDialog:HFJH() /结果值 double Xs=Xs0,Ys=Ys0,Zs=Zs0;/定义 double phi=phi0,omega=omega0,kappa=kappa0;/定义 double x,y,X,Y,Z;/定义X,Y,Z double a1,a2,a3,b1,b2,b3,c1,c2,c3; do

19、uble X_,Y_,Z_; double a11,a12,a13,a14,a15,a16; double a21,a22,a23,a24,a25,a26; double lx,ly; double r_r, dx, dy; cv:Mat A, L, dX; /, A0, L0, A1; /迭代终止阈值 double e = 1.0 / 206265;/单位权中误差double m0; int index = 1; /求取主点坐标 double W = 5344.0; double H = 4008.0; x0 = W / 2 + x0; y0 = H / 2 - y0; do for (si

20、ze_t i = 0; i != pic.size(); +i) x = pici.x; y = pici.y; X = 0; Y = 0; Z = 0; for (size_t j = 0; j < obj.size(); j+) if (objj.index = pici.index) X = objj.X; Y = objj.Y; Z = objj.Z; break; if (abs(x) < 1) char *arg1 = "像点"char *arg2 = "没有找到对应的控制点"ui->textBrowser->appe

21、nd(QString("%1%2%3").arg(QString:fromLocal8Bit(arg1).arg(pici.index).arg(QString:fromLocal8Bit(arg2); continue; /求观测值的改正 r_r = (x - x0)*(x - x0) + (y - y0)*(y - y0); dx = (x - x0)*(k1*r_r + k2*r_r*r_r) + p1*(r_r + 2 * (x - x0)*(x - x0) - 2 * p2*(x - x0)*(y - y0); dy = -(y - y0)*(k1*r_r + k

22、2*r_r*r_r) + p2*(r_r + 2 * (y - y0)*(y - y0) - 2 * p1*(x - x0)*(y - y0); /经过畸变改正后的像点坐标观测值 x = (x - x0) - dx; y = -(y - y0) - dy; /求矩阵R各元素的值 a1 = cos(phi)*cos(kappa) - sin(phi)*sin(omega)*sin(kappa); a2 = -cos(phi)*sin(kappa) - sin(phi)*sin(omega)*cos(kappa); a3 = -sin(phi)*cos(omega); b1 = cos(omega

23、)*sin(kappa); b2 = cos(omega)*cos(kappa); b3 = -sin(omega); c1 = sin(phi)*cos(kappa) + cos(phi)*sin(omega)*sin(kappa); c2 = -sin(phi)*sin(kappa) + cos(phi)*sin(omega)*cos(kappa); c3 = cos(phi)*cos(omega); /求X_，Y_，Z_ X_ = a1*(X - Xs) + b1*(Y - Ys) + c1*(Z - Zs); Y_ = a2*(X - Xs) + b2*(Y - Ys) + c2*(Z

24、 - Zs); Z_ = a3*(X - Xs) + b3*(Y - Ys) + c3*(Z - Zs); /求x，y的近似值 double x_ = -fx*X_ / Z_; double y_ = -fy*Y_ / Z_; /求矩阵L各元素的值 lx = x - x_; ly = y - y_; /求矩阵A各元素的值 a11 = 1 / Z_*(a1*fx + a3*x); a12 = 1 / Z_*(b1*fx + b3*x); a13 = 1 / Z_*(c1*fx + c3*x); a21 = 1 / Z_*(a2*fy + a3*y); a22 = 1 / Z_*(b2*fy +

25、b3*y); a23 = 1 / Z_*(c2*fy + c3*y); a14 = y*sin(omega) - (x / fx*(x*cos(kappa) - y*sin(kappa) + fx*cos(kappa)*cos(omega); a15 = -fx*sin(kappa) - x / fx*(x*sin(kappa) + y*cos(kappa); a16 = y; a24 = -x*sin(omega) - (y / fy*(x*cos(kappa) - y*sin(kappa) - fy*sin(kappa)*cos(omega); a25 = -fy*cos(kappa) -

26、 y / fy*(x*sin(kappa) + y*cos(kappa); a26 = -x; /合成矩阵 double a12 = a11, a12, a13, a14, a15, a16, a21, a22, a23, a24, a25, a26 ; double l2 = lx, ly ; cv:Mat Ai(2, 6, CV_64FC1, a); cv:Mat Li(2, 1, CV_64FC1, l); A.push_back(Ai); L.push_back(Li); /计算 dX = (A.t() * A).inv() * A.t() * L; Xs+=dX.at<doub

27、le>(0,0); Ys+=dX.at<double>(1,0); Zs+=dX.at<double>(2,0); phi+=dX.at<double>(3,0); omega+=dX.at<double>(4,0); kappa+=dX.at<double>(5,0);ui->textBrowser->append(QString("#%1: %2 %3 %4 %5 %6 %7").arg(index).arg(dX.at<double>(0).arg(dX.at<double

28、>(1).arg(dX.at<double>(2).arg(dX.at<double>(3).arg(dX.at<double>(4).arg(dX.at<double>(5); index+;cv:Mat M0 = (A*dX - L).t()*(A*dX - L);m0 = sqrt(M0.at<double>(0) / (2 * pic.size() - 6); A.release(); L.release(); if (index > 100) break; while(fabs(dX.at<double&g

29、t;(0)>0.001|fabs(dX.at<double>(3)>e);ui->textBrowser->append(QString("Result:");ui->textBrowser->append(QString("Xs:%1 Ys:%2 Zs:%3 phi:%4 omega:%5 kappa:%6").arg(Xs).arg(Ys).arg(Zs).arg(phi).arg(omega).arg(kappa);char *log = "单位权中误差："ui->textB

30、rowser->append(QString("%1%2").arg(QString:fromLocal8Bit(log).arg(m0);void HFJHDialog:on_pushButton_clicked()QString qpath = QFileDialog:getOpenFileName(this);string path = qpath.toLocal8Bit();if(readxy(path) = -1)char *log ="读取像点坐标错误XXXXXX"ui->textBrowser->append(QStrin

31、g:fromLocal8Bit(log);return;char *log = "读取像点坐标结束OOOOOO"ui->textBrowser->append(QString:fromLocal8Bit(log);ui->lineEdit->setText(qpath);void HFJHDialog:on_pushButton_2_clicked()QString qpath = QFileDialog:getOpenFileName(this);string path = qpath.toLocal8Bit();if(readXYZ(path)

32、= -1)char *log = "读取控制点坐标错误XXXXXX"ui->textBrowser->append(QString:fromLocal8Bit(log);return;char *log = "读取控制点坐标结束OOOOOO"ui->textBrowser->append(QString:fromLocal8Bit(log);ui->lineEdit_2->setText(qpath);void HFJHDialog:on_pushButton_4_clicked()HFJH();void HFJHDi

33、alog:on_pushButton_5_clicked()HFJHDialog:close();图1 控件后方交会2、特征匹配#include "match.h"#include <opencv.hpp>#include <iostream>#include <fstream>using namespace std;using namespace cv;/构造函数Mmatch:Mmatch(int length, int flag)Mlength = length;Mflag = flag;match();/构建子窗体void Mmat

34、ch:findsubwin(Mat &src, Mat &dst, int i, int j, int length)dst = Mat:zeros(length, length, CV_8UC1);int xl = i - length / 2;int yl = j - length / 2;int xr = i + length / 2;int yr = j + length / 2;int xdl = 0;int ydl = 0;int xdr = length - 1;int ydr = length - 1;if (xl<0)xdl = 0 - xl;xl =

35、0;if (yl < 0)ydl = 0 - yl;yl = 0;if (xr >= src.cols)xdr = (dst.cols -1)- (xr - (src.cols-1);xr = src.cols - 1;if (yr >= src.rows)ydr = (dst.rows-1) - (yr - (src.rows-1);yr = src.rows - 1;Mat sROI = src(Range(yl+1, yr+1), Range(xl+1, xr+1);sROI.copyTo(dst(Range(ydl+1, ydr+1), Range(xdl+1, xd

36、r+1);/匹配算法主体函数int Mmatch:match()string path1 = "E:我的课程数字摄影测量影像数据匹配1BAK.bmp"string path2 = "E:我的课程数字摄影测量影像数据匹配2BAK.bmp"string path = "E:我的课程数字摄影测量影像数据匹配1BAK_points.txt"ifstream fil;vector<point> pointssrc;fil.open(path);if (!fil)cout << "read points erro

37、r" << endl;return -1;int index;point pt;while (!fil.eof()fil >> index >> pt.x >> pt.y;pointssrc.push_back(pt);Mat src1 = imread(path1,CV_LOAD_IMAGE_GRAYSCALE);Mat src2 = imread(path2,CV_LOAD_IMAGE_GRAYSCALE);if (src1.empty() | src2.empty()cout << "read image e

38、rror" << endl;return -1;vector<point> pointsdst;Mat subwinsrc, subwindst;for (int i = 0; i <pointssrc.size(); +i)point tpoint;findsubwin(src1, subwinsrc, pointssrci.x, pointssrci.y, Mlength);findsamename(src2, subwinsrc, tpoint);pointsdst.push_back(tpoint);Mat dst1, dst2;cvtColor

39、(src1, dst1, CV_GRAY2RGB);cvtColor(src2, dst2, CV_GRAY2RGB);CvFont font;cvInitFont(&font,CV_FONT_HERSHEY_PLAIN,1.0,1.0);for (int i = 0; i <pointssrc.size(); +i)char cind5;sprintf(cind,"%d",i+1);line(dst1, Point(pointssrci.x-3, pointssrci.y-3),Point(pointssrci.x+3, pointssrci.y+3), C

40、V_RGB(255, 0, 0);line(dst1, Point(pointssrci.x-3, pointssrci.y+3),Point(pointssrci.x+3, pointssrci.y-3), CV_RGB(255, 0, 0);putText(dst1,cind,Point(pointssrci.x+4,pointssrci.y+5),CV_FONT_HERSHEY_PLAIN,1.0,CV_RGB(255,0,0);line(dst2, Point(pointsdsti.x-3, pointsdsti.y-3), Point(pointsdsti.x+3, pointsds

41、ti.y+3), CV_RGB(255, 0, 0);line(dst2, Point(pointsdsti.x-3, pointsdsti.y+3), Point(pointsdsti.x+3, pointsdsti.y-3), CV_RGB(255, 0, 0);putText(dst2,cind,Point(pointsdsti.x+4,pointsdsti.y+5),CV_FONT_HERSHEY_PLAIN,1.0,CV_RGB(255,0,0);imshow("src", dst1);imshow("dst", dst2);waitKey()

42、;return 0;/寻找同名点函数void Mmatch:findsamename(Mat &src, Mat &subwin, point &dst)Mat subwind;int length = subwin.cols;int helf_length = length / 2;int flag = Mflag;float relate = 0;for (int i = 0; i < src.cols; +i)for (int j = 0; j < src.rows; +j)findsubwin(src, subwind, i, j, length);

43、if (i = helf_length && j = helf_length)relate = compare(subwin, subwind, flag);dst.x = i;dst.y = j;float tmp = compare(subwin, subwind, flag);if (tmp > relate)relate = tmp;dst.x = i;dst.y = j;/特征向量比较函数float Mmatch:compare(Mat &subwin, Mat &subwind,int flag)float sum = 0.0;float av

44、g1 = 0.0, avg2 = 0.0;float sdv1 = 0.0, sdv2 = 0.0;int nr = subwin.rows;int nc = subwin.cols;if(subwin.isContinuous()&&subwind.isContinuous()nc = nc*nr;nr = 1;switch (flag)case 1:/相关函数法for (int i = 0; i < nr; +i)uchar *subwin_ptr = subwin.ptr<uchar>(i);uchar *subwind_ptr= subwind.ptr

45、<uchar>(i);for (int j = 0; j < nc; j+)sum += (float)(*(subwin_ptr+j)* (*(subwind_ptr+j);return sum;case 2:/协方差函数法for (int i = 0; i < nr; +i)uchar *subwin_ptr = subwin.ptr<uchar>(i);uchar *subwind_ptr= subwind.ptr<uchar>(i);for (int j = 0; j < nc; j+)avg1 += *(subwin_ptr+j)

46、;avg2 += *(subwind_ptr+j);avg1 = avg1 / (subwin.rows*subwin.cols);avg2 = avg2 / (subwind.rows*subwind.cols);for (int i = 0; i < nr; +i)uchar *subwin_ptr = subwin.ptr<uchar>(i);uchar *subwind_ptr= subwind.ptr<uchar>(i);for (int j = 0; j < nc; j+)sum += (*(subwin_ptr+j)-avg1)*(*(subw

47、ind_ptr+j) - avg2);return sum;case 3:/相关系数法for (int i = 0; i < nr; +i)uchar *subwin_ptr = subwin.ptr<uchar>(i);uchar *subwind_ptr= subwind.ptr<uchar>(i);for (int j = 0; j < nc; j+)avg1 += *(subwin_ptr+j);avg2 += *(subwind_ptr+j);avg1 = avg1 / (subwin.rows*subwin.cols);avg2 = avg2 /

48、 (subwind.rows*subwind.cols);for (int i = 0; i < nr; +i)uchar *subwin_ptr = subwin.ptr<uchar>(i);uchar *subwind_ptr= subwind.ptr<uchar>(i);for (int j = 0; j < nc; j+)sum += (*(subwin_ptr+j)-avg1)*(*(subwind_ptr+j) - avg2);for (int i = 0; i < nr; +i)uchar *subwin_ptr = subwin.ptr

49、<uchar>(i);uchar *subwind_ptr= subwind.ptr<uchar>(i);for (int j = 0; j < nc; j+)sdv1 += (*(subwin_ptr+j)-avg1)*(*(subwin_ptr+j) - avg1);sdv2 += (*(subwind_ptr+j)-avg2)*(*(subwind_ptr+j) - avg2);sum = sum / sqrtf(sdv1*sdv2);return sum;case 4:/差平方和法for (int i = 0; i < nr; +i)uchar *subwin_ptr = subwin.ptr<uchar>(i);uchar *subwind_ptr= subwind.ptr<uchar>(i);for (int j = 0; j < nc; j+)sum += (float)(*(subwin_ptr+j)-*(subwind_ptr+j)*(*(s