武汉大学测绘学院摄影测量集中实习.docx
《武汉大学测绘学院摄影测量集中实习.docx》由会员分享，可在线阅读，更多相关《武汉大学测绘学院摄影测量集中实习.docx（33页珍藏版）》请在冰点文库上搜索。
武汉大学测绘学院摄影测量集中实习
课程编号：
20162024123课程性质：
专业必修
数字摄影测量实习
实习报告
学院：
测绘学院
专业：
测绘工程
地点：
校
班级：
组号：
：
学号：
教师：
2017年6月19日至2017年7月8日
第1章空三报告
1.1NAT空中三角成果数据流程
1.设置测区基本参数
2.建立相机文件，导入老师给的文件
3.导入外业控制点
4.影像格式转换
5.自动转点
6.进行PATB平差计算
1.2空三加密需要注意
1.像素大小、焦距、像主点坐标、畸变参数，像幅尺寸（数码相机）
2.模型连接点平面及高程限差、地面定向点残差限差、地面检查点残差限差
3.控制点坐标、刺点片及点之记
1.3
空三流程图
1.4空三评查步骤
1.4.1需要平差的原因
1同一物体在不同视角下的图像有极大的不同，物体遮挡而丢失信息
2光照、场景中物体的几何形状和物理特性、物体和摄像机之间的空间关系
3从带畸变和噪声的、二维的、单一的灰度值中提取出尽可能不带畸变和噪声的、三维场景中的诸多因素
1.4.2粗差的影响
1观测值中存在粗差，则某一观测的误差对所有残差都有影响
2一个观测的残差包含着所有粗差的影响
1.4.3平差方法-循环粗差剔除
1先外围加4个点进行平差解算
像方的粗差点排除到一定程度后，开始加预测的控制点，不要一次性刺入所有预测控制点，建议从外向回字形一圈圈刺预测控制点，留区域中的3到4个控制点做检核点。
每刺入一轮都先平差，编辑得差不多再进行下一轮。
2继续调整控制点
关掉粗差探测，勾选验后方差，再平差一次，平差结束。
3上述为精度报告检查平差结果，最后需要到立体测图软件进行检查。
1.4.4自己的方案
在本次测区的15个点中，我分了三轮进行平差，每次选取4个点最后保留3个点作为检查点，难处在于有的点有灰尘遮盖造成精度不准，我的精度指标如下所示。
第2章Mapmatrix进行DOM,DEM生成
2.1外方位元素的确定
2.1.1定向：
确定与物点相对应的像点坐标
本次实习测量了四个框标，故采用仿射变换形式
X=a0+a1*x’+a2*y’;
Y=b0+b1*y’+b2*y’;
2.1.2扫描分辨率为什么不能错：
量测的像点坐标存在着从扫描坐标岛像点坐标的转换，如果没有调整合适的分辨率就会对影像定向产生影响。
2.1.3创建立体像对：
2.1.4相对定向：
1.单独像对相对定向
2.连续像对相对定向
目的是为了恢复摄影时两影像光束的相互关系，从而使同名光线对对相交。
5个元素
2.1.5绝对定向：
借助物空间坐标为已知的控制点来确定空间辅助坐标系与实际物空间坐标的变换关系称为立体模型的绝对定向。
2.1.6后方交会
后方交会：
如果知道每影像的6个外方位元素就可以知道确定被摄物体与航射影像的关系，利用影像围的控制点的空间坐标与影像坐标根据共线方程式，反求影像的外方位元素。
2.1.7后方交会绝对定向刺6个点能成功：
绝对定向有七个参数，空间相似变换的3个角元素，3个平移量，1个比例尺的缩放。
2个平高+1个高程控制点。
独立模型法区域空中三角测量的基本思想，把一个单元模型视为刚体，一个单元模型由两个或更多个立体像对构成，因为一个立体像对需要两影像，
所以如图所示，两影像的拼接必须要有三个同名像点作为公共点，满足X,Y,Z的传递。
连接过程中每个模型可以视为刚体，只能做平移，旋转，缩放。
两相片外方位元素12个，绝对定向的3个角元素，3个平移量，1个比例尺的缩放，相对定向5个元素。
通过相对定向＋绝对定向能够恢复两影像的外方位元素。
2.1.8核线、匹配采样
核线的性质：
同名像点必位于同名核线上。
通过摄影基线与地面所作的平面称为核面核面与影像面交线称为核线同名像点必定在同名核线上。
点击定义核线围，重采样即可完成。
2.1.9影像匹配
2.2DEM编辑生成
2.2.1选取新建的DEM成品自动生成
2.2.2DEM编辑
（对hammer编辑）
2.2.3DEM介绍
Hammer第一条航带的DEM编辑，数字高程模型（DigitalElevationModel），简称DEM，是通过有限的地形高程数据实现对地面地形的数字化模拟（即地形表面形态的数字化表达），它是用一组有序数值阵列形式表示地面高程的一种实体地面模型。
2.2.4问老师得到的对Dem编辑的理解
第一条航带编辑必须要把浮动到眼前，通过3D眼睛看到的明显具有高程的房屋建筑改化到地面高程中去，生成的是地物地貌模型，而不是具有人工建筑的地物模型。
将房屋等人工建筑周围的地形高程进行插，求出该店的实际高程。
方法：
匹配的点，没贴准地面需要编辑
插的点没贴准地面需要编辑
匹配点插方式：
不量测地面的点，用两边的点对中间点插
量测点插方式：
量测地面的点
2.2.5
Dem拼接裁剪
（Dem成图）
2.3DOM的生成编辑
2.3.1DOM概念
DOM-数字正射影像图是对航空（或航天）相片进行数字微分纠正和镶嵌，按一定图幅围裁剪生成的数字正射影像集。
它是同时具有地图几何精度和影像特征的图像。
2.3.2选取新建的DOM产品自动生成
调整扫描分辨率与参数
沿影像边缘生成，且原始影像单独生成DOM
2.3.3规格要求
✓数据能正确衔接
✓建筑等无错开、扭曲变形等现象
✓色调一致、合理
1YLD裁切后的影像（0.1米）（用DSM生成的拼接后的DOM不用单独提交，主要用于DOM建筑等人工地物的变形分析）
2hammer数据拼接后DOM0.5米（不用编辑，主要用于DOM色彩分析）
2.3.4生成72影像的DOM导入易拼图进行正射影像编辑
（影像镶嵌）
（待编辑的正射影像）
2.3.5EPT编辑工具使用总结
将72影像生成的DOM拼接导入EPT进行编辑编辑注意要不过水体，不过房屋，可以借助工具进行镶嵌线编辑，添加点、插入点、移动点、删除点、删除临时线和修补选区。
注意事项：
对于某些四点公用点的镶嵌线编辑
1.在与关键点相连的任意一条镶嵌线上单击左键添加点，开始绘制。
2.到第四条镶嵌线上，单击左键添加结束点。
3.最后单击右键，结束镶嵌线编辑。
4.其中绘制镶嵌线时，开始点和结束点必须捕捉到镶嵌线上，所绘制的镶嵌线必须跨越其他两条镶嵌线，否则无法结束此条镶嵌线的编辑。
2.3.6DOM建筑等人工地物变形对比分析
纠正前后对比图样
纠正前纠正后
进行镶嵌线编辑后要进行正射影像修补，用修补工具围起来变形的房子，右键即可完成修补工作。
2.3.7DOM色彩对比分析
DOM编辑对于某些地域进行了拼接线调整扭曲、变形、错位调整局部调色。
如下面的hammer测区最后生成的DOM图所示，航带衔接的航空图片有着明显的区分镶嵌线，如右图所示，明显看出本条航带有三影像。
（hammer生成的DOM影像）
2.3.8最后生成的产品
可见未经房屋过滤前生成的DEM图像，房屋屋顶以黄褐色标示，月亮岛小区平地多为黄色或者是浅黄色，地势较为低洼的湖泊池塘多用绿色来渲染，象征了高程的不同。
房屋过滤后的DEM影像较为单一，没有了大的凸起，进行了一层过滤。
（经过房屋过滤后的DEM）
（最终的Dom成果图）
第3章DLG的立测生成
3.1DLG的概念
DLG（DigitalLineGraphic，数字线划地图）是与现有线划基本一致的各地图要素的矢量数据集，且保存了各要素间的空间关系信息和相关的属性信息，是4D产品的一种。
3.2本次试验生成的DLG成果，导出为DXF格式
3.3
立测的精度检核
3.4DLG生成
1要求地物丰富
2培养立体感
3“符号化配置”成“MapMatrix2007国标旧版”
3.4.1经验总结
在绘制矢量线划图时，注意要绘制一个高程平面的地物，如在控制点附近可以参考控制点的高程坐标。
实验的每一步都是建立在上一步的基础上，DEM-DOM-DLG这样的顺序。
至于4D的另一个产品DRG数字栅格地图则是DOM,DEM派生出新的可视化信息。
绘制房子与绘制道路平面时要注意在不同的高程，否则立测看起来道路仿佛具有了高程，“浮在”了地面之上，不注意调整高程还会导致不同地物之间的覆盖，不能正确的表达地貌。
在这次实验认识到了摄影测量空三成果的基础性以及对于后面的成果有着重大影响，要求我们每一步都要有合适的精度作为生产标准，理解了平差精度的重要。
最大的收获是培养了立体感，当真真切切的带着立测眼睛看自己的工作时，感觉自己仿佛置身于正在测绘的无人机上一样，以一种高度来俯视测区，山丘的断崖，房屋的起伏都是真真切切的，航空摄影测量以后将会有着很好的前途，逐步替代低精度要求却外业工作很繁重的测绘任务。
实习又认识到了3R之间的异同处，无论是通过全站仪，GPS，无人机，都是为了生产出合格的dlg地图，都可以导入到cass里进行编辑，这都是殊途同归的。
第4章编程实习
4.1实验目的
在掌握共线方程原理的基础上，自己编写后方交会程序，在计算机上调试，输出计算结果：
像片的外方位元素。
4.2具体要求
已知值
像点坐标为扫描坐标，原点在像片的左上角，单位为pixel。
应将坐标原点平移至像片中心，作为像平面坐标系（x轴向右、y轴向上），代入共线方程。
像片高h=4008（pixel）y方向
像宽片w=5344（pixel）x方向
方位元素
x0=47.4857pixel
y0=12.0276pixel
f=4547.9352pixel
对共线方程进行线性化，列立误差方程式，系数应采用严密式，按照最小二乘原则进行解。
参考结果
Xs=3373.40082mm；
Ys=-141.55657mm；
Zs=92.77483mm；
Phi=-0.01205弧度；
Omega=0.09997弧度；
Kappa=-0.04101弧度。
4.3计算原理
如图所示，物点A和摄影中心S在地面摄影测量坐标系中的坐标依次是（X,Y,Z）、（XS,YS,ZS）；像点a在像空间坐标系中的坐标是（x,y,-f）。
那么由共线条件方程知，
其中ai，bi，ci是只含三个独立参数ψ，ω，κ的九个方向余弦。
在方程中共有六个未知参数XS,YS,ZS,ψ,ω,κ，所以有三个不在一条直线上的已知地面点坐标就可以求出像片的这六个外方位元素。
由于共线条件方程是非线性方程，为了便于迭代计算，需要把方程用泰勒级数展开，取一次项得到线型表达式，有
用新的符号表示各偏导数后为
其中（x）、（y）是函数近似值，
是外方位元素近似值的改正数，它们的系数为函数的偏导数。
为了便于推导，令
=
=
=
对于竖直摄影而言，像片的角方位元素都是小值，因而各系数的近似值为
为了提高精度和可靠性，通常需要测量四个或更多的地面控制点和对应的像点坐标，采用最小二乘平差方法解算。
此时像点坐标（x,y）作为观测值，加入相应的偶然误差改正数
，可列出每个点的误差方程式
用矩阵表示为
那么由
4.4算法解释
每一个物方坐标就有一个未知的x,y就可以列出两个误差方程，因为有6个外方位元素所以只有3个点就可以了，本次试验多余的数据作为平差多余计算参与平差解算。
4.5算法流程：
获取已知数据m,x0,y0,f,Xtp,Ytp,Ztp
量测控制点像点坐标x,y
确定未知数初值Xs0,Ys0,Zs0,0,0,0
组成误差方程式并法化
解求外方位元素改正数
检查迭代是否收敛
具体的流程程序框图如下：
确定外方位元素初始值
组成旋转矩阵R
逐点组成误差方程式并法化
所以像点完否
解法方程，求外方位元素改正数
计算改正后的外方位元素
外方位元素改正数是否小于限差
输出计算成果，计算结束
迭代次数小于n
结束并显示错误信息
是
否否
是
否
是
4.6MATLAB源代码
%aa=textread（'物方控制点坐标.txt'）;
fid=fopen（'CONTROL.txt'）;%打开数据总文件
B=textscan（fid,'%f%f%f%f'）;%把每一列的数据读入到读入到单元数组B中
C=[B{2}B{3}B{4}];%从单元数组B中提取每列数据赋值给矩阵C
n=max（size（C））;%N为矩阵的行总数也就是控制点位个数Xa
fid=fopen（'IAMGE.txt'）;%打开数据总文件
J=textscan（fid,'%f%f%f'）;%把每一列的数据读入到读入到单元数组B中
I=[J{2}J{3}];%从单元数组B中提取每列数据赋值给矩阵C%确定读入数据的坐标数目
%Cp为平均值
Cpx1=sum（C）;%控制点的X,Y,Z总和
Xap=Cpx1（1,1）/n;%控制点平均Xaverage
Yap=Cpx1（1,2）/n;
Zap=Cpx1（1,3）/n;
%要求待定参数Xs=3373.40082mm,Ys=-141.55657mm,Zs=92.77483mm,Phi,Omega,Kappa
%像方坐标，控制点，Xa,Ya,Za,delta
%设定外方位元素初始值
x0=47.4857;
y0=12.0276;
f=4547.9352;
Xs=3370;
Ys=-140;
Zs=92;
Phi=0;
Omega=0;
Kappa=0;
%RPhi=[cos（Phi）,0,-sin（Phi）;0,1,0;sin（Phi）,0,cos（Phi）];
%ROmega=[1,0,0;0,cos（Omega）,-sin（Omega）;0,sin（Omega）,cos（Omega）];
%RKappa=[cos（Kappa）,-sin（Kappa）,0;sin（Kappa）,cos（Kappa）,0;0,0,1];
%R=RPhi*ROmega*RKappa;
%a1=R（1,1）;
%a2=R（1,2）;
%a3=R（1,3）;
q=Phi;
w=Omega;
k=Kappa;
zy=0;
while1
forj=n:
-1:
1
i=n-j+1;
a
（1）=cos（q）*cos（k）-sin（q）*sin（w）*sin（k）;
a
（2）=-cos（q）*sin（k）-sin（q）*sin（w）*cos（k）;
a（3）=-sin（q）*cos（w）;
b
（1）=cos（w）*sin（k）;
b
（2）=cos（w）*cos（k）;
b（3）=-sin（w）;
c
（1）=sin（q）*cos（k）+cos（q）*sin（w）*sin（k）;
c
（2）=-sin（q）*sin（k）+cos（q）*sin（w）*cos（k）;
c（3）=cos（q）*cos（w）;
%以i开始迭代Control坐标
%开始迭代
Xbar（i）=a
（1）*（C（i,1）-Xs）+b
（1）*（C（i,2）-Ys）+c
（1）*（C（i,3）-Zs）;
Ybar（i）=a
（2）*（C（i,1）-Xs）+b
（2）*（C（i,2）-Ys）+c
（2）*（C（i,3）-Zs）;
Zbar（i）=a（3）*（C（i,1）-Xs）+b（3）*（C（i,2）-Ys）+c（3）*（C（i,3）-Zs）;
H（i）=-Zbar（i）;
%H赋予初值
%误差方程的系数aij改为xsij
xs（2*i-1,1）=-f/H（i）;
xs（2*i-1,2）=0;
xs（2*i-1,3）=-（I（i,1）-x0）/H（i）;
xs（2*i-1,4）=-（f+（I（i,1）-x0）*（I（i,1）-x0）/f）;
xs（2*i-1,5）=-（I（i,1）-x0）*（I（i,2）-y0）/f;
xs（2*i-1,6）=I（i,2）-y0;
xs（2*i,1）=0;
xs（2*i,2）=-f/H（i）;
xs（2*i,3）=-（I（i,2）-y0）/H（i）;
xs（2*i,4）=-（I（i,1）-x0）*（I（i,2）-y0）/f;
xs（2*i,5）=-（f+（I（i,2）-y0）*（I（i,2）-y0）/f）;
xs（2*i,6）=-（I（i,1）-x0）;
%L=[x0-I（i,1）,y0-I（i,2）];
L（2*i-1,1）=x0-I（i,1）;
L（2*i,1）=y0-I（i,2）;
%此时L为矩阵2n*1
end
%=[deltaXs;deltaYs;deltaZs;deltaq;deltaw;deltak];
N=xs'*xs;
X=pinv（N）*xs'*L;
deltaXs=X（1,1）;
deltaYs=X（2,1）;
deltaZs=X（3,1）;
deltaq=X（4,1）;
deltaw=X（5,1）;
deltak=X（6,1）;
Xs=deltaXs+Xs;
Ys=deltaYs+Ys;
Zs=deltaZs+Zs;
q=q+deltaq;
w=w+deltaw;
k=k+deltak;
zy=zy+1;
PD
（1）=deltaq/pi*60;
PD
（2）=deltaw/pi*60;
PD（3）=deltak/pi*60;
ifPD
（1）<0.1&&PD
（2）<0.1&&PD（3）<0.1
break;
end
end
Sigmav=0;
forzx=1:
n:
1
v（2*zx-1,1）=xs（1,:
）*X-L（2*zx-1,1）;
%v矩阵为2n行一列
Sigmav=Sigmav+v（2*zx-1,1）;
end
m0=sqrt（Sigmav/（n-6））;
fid111=fopen（'jieguo.txt','wt'）;
%%%%需要改文件名称的地方
fprintf（fid111,'%f%f%f%f%f%f\n',Xs,Ys,Zs,q,w,k）;%%%%ve：
需要导出的数据名称
4.7影像匹配原理
常用的影像匹配方法可分为两种：
基于影像灰度匹配的影像匹配算法和基于影像特征的匹配算法。
基于影像灰度的影像匹配算法以左、右像片含有相应影像的目标区和搜索区中的像元的灰度作为影像匹配的基础，利用某种相关度量，例如用协方差或相关系数为最大来判定左右影像中的相应像点。
影像匹配可以用二维的窗口，也可用一维窗口（按核线排列）的像元灰度参与计算。
一维相关时，搜索仅在一维方向进行，具有计算量小，精度好的特点。
基于灰度的影像匹配方法中，以最小二乘影像匹配算法精度最高。
基于影像特征的影像匹配算法利用了数字影像处理中的特征边缘提取技术，先检测出左右影像中的特征点，并利用匹配算法或其它方法匹配出特征点对。
基于特征的影像匹配有较高的可靠性，但匹配的精度低于基于灰度的最小二乘影像匹配算法。
一种比较好的方案是将基于特征的影像匹配与基于灰度的最小二乘影像匹配结合在一起，形成分层由粗到精的匹配方案。
从而使最终的匹配结果既有较高的可靠性，又有最小二乘影像匹配的高精度。
灰度匹配算法：
设两个随机变量分别代表左、右数字影像中的一个，两个像元灰度阵列为相应匹配影像阵列，其中点即为相应像点。
辐射畸变纠正
几何畸变纠正
重采样
4.8主程序
function[WpointR,BpointR]=LmsArray0824（Baseimg,Warpimg,Bpoint,Wpoint,Num_it,r）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%LMS为最小二乘影像匹配函数
%%%%%输入变量：
Baseimg为基础影像，Warpimg为待校正影像，Bpoint,Wpoint分别为基准影像和待匹配影像特征点对,[NX2]
%%%%%Num_it为迭代次数，
%%%%%r为匹配窗口半径。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%WpointR\BpointR为匹配后结果。
W=[];
W=imread（'C:
\Users\Administrator\Desktop\数字摄影测量编程实习资料-D方向\匹配\1BAK.bmp'）;
W=rgb2gray（W）;
Warpimg=W;
B=[];
B=imread（'C:
\Users\Administrator\Desktop\数字摄影测量编程实习资料-D方向\匹配\2BAK.bmp'）;
Bresample=[];
Bresample=imresize（B,0.5,'bicubic'）;
Baseimg=Bresample;
Bpoint=Bresamplepoint;
Num_it=10;
r=7;
if（size（Bpoint,1）==0）
error（'输入特征点个数不能为零'）;
end
if（size（Bpoint,1）~=size（Wpoint,1））
error（'输入错误'）;
end
BpointR=zeros（size（Bpoint,1）,2）;
WpointR=zeros（size（Wpoint,1）,2）;
BW=zeros（size（Bpoint,1）,1）;
if（size（size（Baseimg）,2）>2）
Baseimg=rgb2gray（Baseimg）;
end
if（size（size（Warpimg）,2）>2）
Warpimg=rgb2gray（Warpimg）;
end
[Bimagesizer,Bimagesizec]=size（Baseimg）;
[Wimagesizer,Wimagesizec]=size（Warpimg）;
Baseimg=double（Baseimg）;
Warpimg=double（Warpimg）;
%%%%%%%%%%%%%%%%%%高斯滤波
BA=fspecial（'gaussian',[11,11]）;
Baseimg=filter2（BA,Baseimg）;
Warpimg=filter2（BA,Warpimg）;
fornum_Bpoint=1:
size（Bpoint,1）
Lo=zeros（Num_it+1,1）;
%Lo
（1）=0;
%if（size（size（Baseimg）,2）>2）
%Baseimg=rgb2gray（Baseimg）;
%end
%if（size（size（Warpimg）,2）>2）
%Warpimg=rgb2gray（Warpim
升级会员 