测绘工程matlab论文.docx
- 文档编号:16729337
- 上传时间:2023-07-16
- 格式:DOCX
- 页数:16
- 大小:345.98KB
测绘工程matlab论文.docx
《测绘工程matlab论文.docx》由会员分享,可在线阅读,更多相关《测绘工程matlab论文.docx(16页珍藏版)》请在冰点文库上搜索。
测绘工程matlab论文
课程编号:
课程性质:
MATLAB及其应用
课程论文
学院:
测绘学院
专业:
测绘工程
:
学号:
注:
这是《误差理论与测量平差基础》这一课程在学习其基本原理之后为解决实际问题而编写的一个程序。
一、题目内容
要求:
对于给出的导线控制网图2,根据已知条件和观测数据,设计平差方案,编制平差程序,给出平差结果,评定精度。
题目:
A、B为已知点,已知点坐标已给出如下。
为待定点。
同精度观测了24个角度,测角中误差。
测量了17条边长,其观测结果及中误差列于表4中,试按间接平差法求:
〔1〕平差后单位权中误差;
〔2〕14个待定点的的坐标平差值及中误差;
〔3〕平差后边的相对中误差。
已知点坐标:
=730.024,
=126.040,
,
表一:
导线网角度观测值
编号
角度观测值L
编号
角度观测值L
编号
角度观测值L
1
9
17
2
10
18
3
11
19
4
12
20
5
13
21
6
14
22
7
15
23
8
16
24
表二:
导线网边长观测值及精度
编号
边观测值S〔m〕
中误差〔mm〕
编号
边长观测值S〔m〕
中误差〔mm〕
编号
边长观测值S〔m〕
中误差〔mm〕
1
13
7
12
13
2
14
8
14
3
9
15
4
10
16
5
11
17
6
12
图一:
二、程序编写思路
1、数学模型
◆函数模型观测方程
误差方程
◆随机模型
◆方程的解与精度评定
令
法方程即为
可得解为:
故可得平差结果:
从而解得单位权中误差为
2、程序解析
◆计算框图
三、程序运行结果
◆平差后单位权中误差:
=5.30(由于程序取的先验单位权中误差为10,两者相比较,故可以认为这个结果是可以接受的)
◆14个待定点的坐标平差值及中误差见下表:
表三:
点号
坐标平差值〔m〕
中误差〔mm〕
点号
坐标平差值〔m〕
中误差〔mm〕
P1
X
P8
X
Y
Y
P2
X
P9
X
Y
Y
P3
X
P10
X
Y
Y
P4
X
P11
X
Y
Y
P5
X
P12
X
Y
Y
P6
X
P13
X
Y
Y
P7
X
P14
X
Y
Y
表七:
点号
中误差〔mm〕
点号
中误差〔mm〕
点号
中误差〔mm〕
1
6
11
2
7
12
3
8
13
4
9
14
5
10
◆平差后
边的相对中误差为:
四、程序源代码
这是要输入的文件
◆%%从文件中输入数据
fid1=fopen('E:
\学习\平差\S3.txt','rt');
[S1,count]=fscanf(fid1,'%f',[724]);%#ok
fclose(fid1);
S1=S1';
fid2=fopen('E:
\\学习\平差\S4.txt','rt');
[S2,count]=fscanf(fid2,'%f',[517]);
fclose(fid2);
S2=S2';
[S3]=convert(S1);
[X0]=daoxian(S2,S3);
[B,l,P]=BBll(X0,S2,S1,S3);
x=inv(B'*P*B)*B'*P*l;
X=X0'+x/1000;
disp(X);
◆%%计算B,l,P矩阵
function[B,l,P]=BBll(X0,S2,S1,S3)%#ok
X10=[
855.111
];
fori=1:
1:
7
ifi<2
S=sqrt((X0
(1)-X10
(1))^2+(X0
(2)-X10
(2))^2);
B(1,1)=(X0
(1)-X10
(1))/S;
B(1,2)=(X0
(2)-X10
(2))/S;
l(1,1)=1000*(S2(1,2)-S);
elseifi>6
S=sqrt((X10
(1)-X0(11))^2+(X10
(2)-X0(12))^2);
B(i,11)=-(X10
(1)-X0(11))/S;%#ok
B(i,12)=-(X10
(2)-X0(12))/S;%#ok
l(i,1)=1000*(S2(7,2)-S);%#ok
elseifi>1&&i<7
S=sqrt((X0(2*i-1)-X0(2*(i-1)-1))^2+(X0(2*i)-X0(2*(i-1)))^2);
B(i,2*(i-1)-1)=-(X0(2*i-1)-X0(2*(i-1)-1))/S;%#ok
B(i,2*(i-1))=-(X0(2*i)-X0(2*(i-1)))/S;%#ok
B(i,2*i-1)=(X0(2*i-1)-X0(2*(i-1)-1))/S;%#ok
B(i,2*i)=(X0(2*i)-X0(2*(i-1)))/S;%#ok
l(i,1)=1000*(S2(i,2)-S);%#ok
end
end
fori=8:
1:
16
ifi<9
S=sqrt((X0(13)-X10(3))^2+(X0(14)-X10(4))^2);
B(8,13)=(X0(13)-X10(3))/S;
B(8,14)=(X0(14)-X10(4))/S;
l(8,1)=1000*(S2(8,2)-S);
elseifi>15
S=sqrt((X10(3)-X0(27))^2+(X10(4)-X0(28))^2);
B(16,27)=-(X10(3)-X0(27))/S;
B(16,28)=-(X10(4)-X0(28))/S;
l(16,1)=1000*(S2(16,2)-S);
elseifi<16&&i>8
S=sqrt((X0(2*(i-1)-1)-X0(2*(i-1)-3))^2+(X0(2*(i-1))-X0(2*(i-1)-2))^2);
B(i,2*(i-1)-3)=-(X0(2*(i-1)-1)-X0(2*(i-1)-3))/S;
B(i,2*(i-1)-1)=(X0(2*(i-1)-1)-X0(2*(i-1)-3))/S;
B(i,2*(i-1)-2)=-(X0(2*(i-1))-X0(2*(i-1)-2))/S;
B(i,2*(i-1))=(X0(2*(i-1))-X0(2*(i-1)-2))/S;
l(i,1)=1000*((S2(i,2)-S));
end
end
S=sqrt((X0(19)-X0(5))^2+(X0(20)-X0(6))^2);
B(17,5)=-(X0(19)-X0(5))/S;
B(17,19)=(X0(19)-X0(5))/S;
B(17,6)=-(X0(20)-X0(6))/S;
B(17,20)=(X0(20)-X0(6))/S;
l(17,1)=1000*(S2(17,2)-S);
row=206.265;
fori=1:
1:
24
[M1]=fangweijiao(S1(i,6),S1(i,5),X0);
[M2]=fangweijiao(S1(i,6),S1(i,7),X0);
ifS1(i,7)>0
B(17+i,2*S1(i,7)-1)=-row*M2
(1);
B(17+i,2*S1(i,7))=row*M2
(2);
end
ifS1(i,5)>0
B(17+i,2*S1(i,5)-1)=row*M1
(1);
B(17+i,2*S1(i,5))=-row*M1
(2);
end
ifS1(i,6)>0
B(17+i,2*S1(i,6)-1)=row*M2
(1)-row*M1
(1);
B(17+i,2*S1(i,6))=-(row*M2
(2)-row*M1
(2));
end
L0=M2(3)-M1(3);
ifL0<0
L0=360+L0;
end
l(17+i,1)=3600*(S3(i,2)-L0);
end
fori=1:
1:
17
P(i,i)=100/S2(i,3)^2;%#ok
end
fori=18:
41
P(i,i)=1;%#ok
end
end
◆%%单位转换
function[S3]=convert(S1)
fori=1:
24
S3(i,1)=S1(i,1);
S3(i,2)=(S1(i,2)+S1(i,3)/60+S1(i,4)/3600);
end
◆%%计算坐标近似值
function[X0]=daoxian(S2,S3)%#ok
X10=[
855.111
];
t=180+atand((X10
(2)-X10(4))/(X10
(1)-X10(3)));
t1=t;
fori=1:
1:
6
ifS2(i,5)>0
t1=t1+S3(S2(i,4),2)+S3(S2(i,5),2)-180;
else
t1=t1+S3(S2(i,4),2)-180;
end
ift1<0
t1=t1+360;
end
ifi<2
X0
(1)=X10
(1)+S2(1,2)*cosd(t1);
X0
(2)=X10
(2)+S2(1,2)*sind(t1);
end
ifi>=2
X0(2*i-1)=X0(2*(i-1)-1)+S2(i,2)*cosd(t1);
X0(2*i)=X0(2*(i-1))+S2(i,2)*sind(t1);
end
end
t2=t;
t2=t2-S3(12,2);
ift2<0
t2=t2+360;
end
fori=8:
1:
15
ifi<9
X0(2*(i-1)-1)=X10(3)+S2(i,2)*cosd(t2);
X0(2*(i-1))=X10(4)+S2(i,2)*sind(t2);
end
ifi>8
t2=t2+S3(S2(i,4),2)-180;%#ok
ift2<0
t2=t2+360;
end
X0(2*(i-1)-1)=X0(2*(i-2)-1)+S2(i,2)*cosd(t2);
X0(2*(i-1))=X0(2*(i-2))+S2(i,2)*sind(t2);
end
end
◆%%计算方位角
function[M]=fangweijiao(a,b,X0)%ÓÃ-1´ú±íAµã£¬ÓÃ0´ú±íBµã
X10=[
855.111
];
ifa>-1&&a<1
X1=X10
(1);%#ok
Y1=X10
(2);%#ok
end
ifb>-1&&b<1
X2=X10
(1);%#ok
Y2=X10
(2);%#ok
end
ifa>-2&&a<0
X1=X10(3);%#ok
Y1=X10(4);%#ok
end
ifb>-2&&b<0
X2=X10(3);%#ok
Y2=X10(4);%#ok
end
ifa>0
X1=X0(2*a-1);
Y1=X0(2*a);
end
ifb>0
X2=X0(2*b-1);
Y2=X0(2*b);
end
deltaX=X2-X1;
deltaY=Y2-Y1;
S=deltaX^2+deltaY^2;
xishux=deltaY/S;
xishuy=deltaX/S;
ifdeltaX>0&&deltaY>0
T=atand(deltaY/deltaX);
end
ifdeltaX<0&&deltaY>0
T=180-atand(deltaY/(-deltaX));
end
ifdeltaX<0&&deltaY<0
T=270-atand((-deltaX)/(-deltaY));
end
ifdeltaX>0&&deltaY<0
T=360-atand(deltaY/(-deltaX));
end
M
(1)=xishux;
M
(2)=xishuy;
M(3)=T;
end
◆%%输出结果
fid=fopen('C:
\DocumentsandSettings\Administrator\桌面\输出结果2.txt','wt');
fprintf(fid,'点号X(m)Y(m)中误差X(cm)中误差Y(cm)点位中误差\n');
fprintf(fid,'%d%11.4f%11.4f%5.4f%5.4f%5.4f\n',(jieguo2)');
fclose(fid);
◆关于matlab的画图的应用,由于本程序不需要用到,故,单独写了个程序如下:
t=-3*pi:
pi/10:
3*pi;%自变量数组中,存在0值
y=sin(t)./t;%在t=0处,将产生非数。
tt=t+(t==0)*eps;%逻辑数组参与运算,使0元素被一个“机器零”小数代替
yy=sin(tt)./tt;%用sin(eps)/eps近似替代sin(0)/0极限
subplot(1,2,1),plot(t,y),axis([-9,9,-0.5,1.2]),%非数在绘图中的剪裁作用
xlabel('t'),ylabel('y'),title('残缺图形')
subplot(1,2,2),plot(tt,yy),axis([-9,9,-0.5,1.2])
xlabel('tt'),ylabel('yy'),title('正确图形')
运行如下:
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 测绘 工程 matlab 论文