程序及运算结果.docx

资源描述

程序及运算结果.docx

《程序及运算结果.docx》由会员分享，可在线阅读，更多相关《程序及运算结果.docx（22页珍藏版）》请在冰豆网上搜索。

程序及运算结果.docx

程序及运算结果

MATLAB程序展示——水准网间接平差

function[v,ch,cx]=szw（s,h,B,x0,d,n,t,f）%此函数用于计算改正数，高差中误差，高程中误差

p=diag（1./s）;%定义权阵

disp（'P='）

disp（p）

l=h-B*x0-d;

W=B'*p*l;

Nbb=B'*p*B;

disp（'NBB='）

disp（Nbb）

x=（inv（Nbb））*W;

disp（'改正数'）

v=B*x-l;%改正数

disp（v）;

c0=sqrt（（v'*p*v）/（n-t））;%单位中误差

Nbb=B'*p*B;

Qh=f*（inv（Nbb））*f';%h5的协因数阵

disp（'高差平差值中误差'）

ch=c0*sqrt（Qh）;%高差平差值中误差

disp（ch）;

Nbbn=inv（Nbb）;%求逆矩阵

disp（'高程平差值中误差'）

cx=c0*sqrt（diag（Nbbn））;%高程平差值中误差

disp（cx）;

return

loadB.txt;

loadd.txt;

loads.txt;

loadx0.txt;

loadh.txt;

loadn.txt;

loadt.txt

loadf.txt

[v,ch,cx]=szw（s,h,B,x0,d,n,t,f）;

计算结果展示

MATLAB程序展示——导线网间接平差

loadn.txt%观测总数

loadt.txt%必要观测数

loaddws.txt%导入点位数

loadcsgs.txt%导入参数个数

loadbs.txt%导入待测边数

loaddcdh1.txt%%第一个待测点号

loadx1y1.txt%导入已知点位一的坐标

loadxydws.txt%导入末点位的坐标

loaddeg.txt%导入由观测角计算而得的方位角

loads.txt%导入边长观测值

loadgcj.txt%各点位观测角

digits（12）;

digits（15）;

x1=x1y1（1,1）;

y1=x1y1（1,2）;

x0=zeros（dws,1）;

y0=zeros（dws,1）;

a=zeros（1,3）;

（1）=x1;

（1）=y1;

fori=2:

（dws-1）

a=deg（i-1,:

）;

[x2,y2]=zbzs（x1,y1,s（i-1）,a）;

x1=x2;

y1=y2;

x0（i）=x1;

y0（i）=y1;

end

x1=xydws（1,1）;

y1=xydws（1,2）;

x0（dws）=x1;%各点位x的近似值

y0（dws）=y1;%各点位y的近似值

s0=zeros（dws-1,1）;%各各点位之间距离的近似值

fori=1:

（dws-1）

s0（i）=sqrt（（x0（i）-x0（i+1））.^2+（y0（i）-y0（i+1））.^2）;

end

[B1,L1]=jsB1L1（x0,y0,s0,s,bs,dcdh1,csgs）;

[B2,L2,ajs]=jsB2L2（x0,y0,s0,dws,dcdh1,csgs,gcj,bs）;

B=[B1;B2];

L=[L1;L2];

juw=5;%测角中误差

buw=0.5*sqrt（s'）;%测边中误差

P=zeros（1,length（gcj））;

fori=1:

length（gcj）

P（i）=1./（juw.^2）;

end

P=[1./（buw.^2）,P];

P=diag（P）*（juw.^2）;

NBB=B'*P*B;

W=B'*P*L;

x=pinv（NBB）*W;

disp（'[xEyExFyF]='）

disp（x）

X=zeros（1,length（x））;

X=[x

（1）+x0

（2）x

（2）+y0

（2）x（3）+x0（3）x（4）+y0（3）];

disp（'E点坐标平差值='）

disp（X（1,1:

2））

disp（'F点坐标平差值='）

disp（X（1,3:

4））

V=B*x-L;

uwef=zwc（V,P,n,t,B）;%各点位精度

disp（'点位精度uweuwf='）

disp（uwef）

disp（'各观测量改正数[Vb1Vb2Vb3Vb4Vj1Vj2Vj3]='）

disp（V）

spcz=s+V（1:

3,1）;%边长观测值平差值

disp（'边长观测平差值S1S2S3='）

disp（spcz）

gcjpcz=gcj+V（4:

7,1）;

disp（'角度观测平差值J1J2J3J4='）

disp（gcjpcz）

functionA=deg2rad（a）%角度转弧度

A=a*[pi/180;pi/10800;pi/648000];

Return

function[x2,y2]=zbzs（x1,y1,s,a）%坐标正算

A=deg2rad（a）;

x2=x1+s*cos（A）;

y2=y1+s*sin（A）;

Return

function[B1,L1]=jsB1L1（x0,y0,s0,s,bs,dcdh1,csgs）%计算边长系数矩阵B

B1=zeros（bs,2*（bs+1））;

L1=s-s0;

fori=1:

k=1;

forj=1+2*（i-1）:

（2*csgs+2*（i-1））

B1（i,j）=（x0（i+1）-x0（i））*（-1）^k./s0（i）;

B1（i,j+1）=（y0（i+1）-y0（i））*（-1）^k./s0（i）;

k=k+1;

end

B1=B1（1:

bs,2*dcdh1-1:

2*dcdh1-2+2*csgs）;

return

function[B2,L2,ajs]=jsB2L2（x0,y0,s0,dws,dcdh1,csgs,gcj,bs）%计算角度的系数矩阵

B2=zeros（dws,2*dws）;

axs=zeros（dws+1,4）;

fori=2:

dws

axs（i,1）=3600*（y0（i）-y0（i-1））./（s0（i-1）.^2）;

axs（i,2）=-3600*（x0（i）-x0（i-1））./（s0（i-1）.^2）;

axs（i,3）=3600*（y0（i）-y0（i-1））./（s0（i-1）.^2）;

axs（i,4）=-3600*（x0（i）-x0（i-1））./（s0（i-1）.^2）;

end

fori=1:

dws

ifi==1

B2（i,1:

4）=axs（i,:

）+axs（i+1,:

）;

else

ifi==dws

B2（i,2*dws-3:

2*dws）=axs（i,:

）+axs（i+1,:

）;

else

B2（i,1+2*（i-2）:

4+2*（i-2））=axs（i,:

）;

B2（i,1+2*（i-1）:

4+2*（i-1））=B2（i,1+2*（i-1）:

4+2*（i-1））+axs（i+1,:

）;

end

B2=B2（1:

dws,2*dcdh1-1:

2*dcdh1-2+2*csgs）;

ajs=zeros（bs+2,1）;

ajs

（1）=46.74972222;

fori=2:

bs+1

ajs（i）=atan（（y0（i-1）-y0（i））./（x0（i-1）-x0（i）））*180./pi;%近似方位角

end

ajs（bs+2）=144.7675;

L2=zeros（bs+1,1）;

fori=1:

bs+1

L2（i）=gcj（i）-180+ajs（i）-ajs（i+1）;

ifL2（i）>90

L2（i）=L2（i）-180;

elseifL2（i）<-90

L2（i）=L2（i）+180;

end

Return

function[vs,cx]=gzs（B,L,P）%求观测值的改正数及坐标参数的改正数

W=B'*P*L;

Nbb=B'*P*B;

cx=（pinv（Nbb））*W;

vs=B*cx-L;

return

functionuwef=zwc（v,P,n,t,B）%计算点位精度

uw0=sqrt（v'*P*v/（n-t））;

Qxx=pinv（B'*P*B）;

Qxx=diag（Qxx）;

uwef=sqrt（uw0.^2*[（Qxx

（1）+Qxx

（2））;（Qxx（3）+Qxx（4））]）;

return

计算结果展示

MATLAB程序展示——三角网间接平差

loadxqs.txt%各起始点x值

loadyqs.txt%各起始点y值

loadgcj.txt%观测角

loada0.txt%近似坐标方位角

loadL1.txt

loadL2.txt

loads0.txt

loadB1.txt

loadB2.txt

loadn.txt%观测总数

loadt.txt%必要观测数

[x0,y0]=jszb（xqs,yqs,gcj）;

P=DP（s0）;

disp（'p='）

disp（P）

B=[B1;B2];

L=[L1;L2];

NBB=B'*P*B;

disp（'NBB='）

disp（NBB）

W=B'*P*L;

x=（inv（NBB））*W;

disp（'坐标改正数'）

disp（x）

x=x./10;

X=[x0

（1）+x

（1）;y0

（1）+x

（2）;x0

（2）+x（3）;y0

（2）+x（4）];

disp（'坐标平差值[X1Y1X2Y2]T='）

disp（X）

x=x*10;

v=B*x-L;

S=s0+v（1:

7,1）./10;

disp（'边长平差值P1AP1BP1CP1P2P2AP2CP2D'）

disp（S）

GCJ=gcj;

GCJ（1:

18,3）=gcj（1:

18,3）+v（1:

18,1）;

disp（'观测角的平差值'）

disp（GCJ）

uwx=jszwc（v,n,t,NBB）;

uwp1=sqrt（uwx

（1）.^2+uwx

（2）.^2）;

disp（'P1点的点位中误差'）

disp（uwp1）

uwp2=sqrt（uwx（3）.^2+uwx（4）.^2）;

disp（'P2点的点位中误差'）

disp（uwp2）

functionA=deg2rad（a）%角度转弧度

A=a*[pi/180;pi/10800;pi/648000];

Return

function[x0,y0]=jszb（xqs,yqs,gcj）%近似坐标

gcj=deg2rad（gcj）;

x0=zeros（2,1）;

y0=zeros（2,1）;

（1）=（xqs

（2）./tan（gcj（16））+xqs（3）./tan（gcj（18））-yqs（3）+yqs

（2））./（1./tan（gcj（18））+1./tan（gcj（16）））;

（2）=（xqs（4）./tan（gcj（7））+xqs（3）./tan（gcj（9））-yqs（4）+yqs（3））./（1./tan（gcj（9））+1./tan（gcj（7）））;

（1）=（yqs

（2）./tan（gcj（16））+yqs（3）./tan（gcj（18））+xqs（3）-xqs

（2））./（1./tan（gcj（18））+1./tan（gcj（16）））;

（2）=（yqs（4）./tan（gcj（7））+yqs（3）./tan（gcj（9））+xqs（4）-xqs（3））./（1./tan（gcj（9））+1./tan（gcj（7）））;

Return

functionP=DP（s0）%定权

P=zeros（1,25）;

P（1,1:

7）=1./s0;

P（1,8:

25）=1;

P=diag（P）;

Return

functionuwx=jszwc（v,n,t,NBB）%计算中误差

v=v（8:

25,1）;

uw0=sqrt（v'*v./（n-t））;

QXX=diag（inv（NBB））;

uwx=uw0*sqrt（QXX）;

return

运算结果展示

[控制网概况]

计算软件：

南方平差易2005

网名：

三角网

计算日期：

2015年11月12日星期四

观测人：

记录人：

计算者：

检查者：

测量单位：

备注：

平面控制网等级：

国家四等，验前单位权中误差：

2.50（s）

已知坐标点个数：

未知坐标点个数：

未知边数：

最大点位误差[P2]=27.4580（m）

最小点位误差[P1]=27.2219（m）

平均点位误差=27.3400（m）

最大点间误差=49.3635（m）

最大边长比例误差=267

平面网验后单位权中误差=1053.86（s）

[闭合差统计报告]

序号:

<1>:

中点多边形

路径：

[]

极条件闭合差=1764.9,限差=±4.1

序号:

<2>:

中点多边形

路径：

[路径：

[]]

极条件闭合差=1196.1,限差=±4.3

序号:

<3>:

闭合导线

路径：

[B-P1-A]

角度闭合差=-2433.62（s）,限差=±8.66（s）

序号:

<4>:

闭合导线

路径：

[P1-P2-A]

角度闭合差=2439.16（s）,限差=±8.66（s）

序号:

<5>:

闭合导线

路径：

[P2-D-A]

角度闭合差=-36.66（s）,限差=±8.66（s）

序号:

<6>:

闭合导线

路径：

[P1-C-B]

角度闭合差=2405.56（s）,限差=±8.66（s）

序号:

<7>:

闭合导线

路径：

[C-P2-P1]

角度闭合差=-2404.72（s）,限差=±8.66（s）

序号:

<8>:

闭合导线

路径：

[C-D-P2]

角度闭合差=-30.28（s）,限差=±8.66（s）

[起算点数据表]

点名

X（m）

Y（m）

H（m）

备注

9684.2800

43836.8200

10649.5500

31996.5000

19063.6600

37818.8600

17814.6300

49923.1900

[方向观测成果表]

测站

照准

方向值（dms）

改正数（s）

平差后值（dms）

备注

0.000000

23.663028

-1611.89

23.393839

90.710084

-732.79

90.584805

122.157306

-404.95

122.092811

0.000000

29.882056

-2496.88

29.464368

59.978362

-2381.46

59.584216

0.000000

126.240028

-937.15

126.082313

246.382000

-475.96

246.302404

312.963861

-1847.72

313.055089

0.000000

22.045278

-534.32

21.555846

88.814222

-1530.49

88.561173

118.790750

-1900.37

118.472713

0.000000

130.053944

-51.73

130.044771

247.433444

324.44

247.485888

313.350333

36.79

313.354012

0.000000

31.172944

-740.67

31.050877

59.072750

-184.90

59.042260

[平面点位误差表]

点名

长轴（m）

短轴（m）

长轴方位（dms）

点位中误差（m）

备注

21.4602

16.7479

57.364786

27.2219

21.6957

16.8298

62.051262

27.4580

[平面点间误差表]

点名

长轴MT（m）

短轴MD（m）

D/MD

长轴方位T（dms）

平距D（m）

备注

27.2219

17.9886

411

57.364786

7399.3840

27.4580

18.4529

321

62.051262

5924.5510

27.2219

21.4005

276

57.364786

5903.4463

27.2219

17.9886

411

57.364786

5889.8917

27.2219

21.4005

276

57.364786

7486.2398

27.2219

18.6105

316

57.364786

7463.1500

34.9053

28.0773

267

55.590287

5940.3263

[控制点成果表]

点名

X（m）

Y（m）

H（m）

备注

9684.2800

43836.8200

已知点

10649.5500

31996.5000

已知点

13194.6667

37323.1368

19063.6600

37818.8600

已知点

15580.1434

44419.1418

17814.6300

49923.1900

已知点

展开阅读全文