水准网程序设计.docx

水准网程序设计.docx

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

水准网程序设计

水准网程序设计（总19页）

第五章水准网程序设计

概述

水准网是为了确定地面点的高程而布设的控制网，网中的观测值是高程控制点之间的高差。

为了统一全国的高程系统，我国采用黄海平均海水面作为全国高程系统的基准面，在该面上的任一点，其高程为零。

水准网中的任一点高程以及点与点之间的高差应属于正常高系统。

但是，水准测量中水准仪逐站测量并累加得到的原始高差并非正常高差，对于高精度的水准网，应在原始数据基础上添加尺长改正、正常水准网面不平行、球气差改正等系统改正，才能得到正常高差，受篇幅所限，本章不讨论观测值的归算问题，本章提到的观测值均假定是已经归算后的正常高差观测值。

由于观测值存在误差，实际工作中需要若干已知点作为平差的基准控制点。

水准网平差的目的就是求解各观测值的最佳估值及评定未知量的精度。

本章开始，我们将开始接触到测量数据处理程序的设计与开发。

测量数据因其数据量大、繁杂等特殊性，往往需要频繁对矩阵进行处理计算。

然而在MATLAB软件中，矩阵的计算变得格外简单，这为广大测绘工作者提供了实用性强、简单方便的数据处理平台。

平差程序就是将平差计算过程程序化，综合考虑，选择参数平差模型作为水准网平差的主要模型。

测量程序设计一般包含程序功能设计、平差模型选择、算法选择等内容。

在本章的结构的组织上，先对水准原理进行简要分析，在此基础上按不同功能模块设计函数，最后对程序进行分析验证。

水准路线处理程序设计

水准网平差函数类设计

1.类设计

Main%函数主体程序

calculateHO（）%计算近似高程

ca_ATPA（）%计算法方程系数A，权阵GP，常数项GL

ca_V（）%计算高差改正数

Printlevledata（）%输出结果

2.函数说明

具体各个成员变量的含义在后面程序中会有注释，在此对主函数main做简要强调。

主函数main中包括数据的读取、检查数据格式、提取相关观测值、各类成员函数和最小二平差。

是本程序运行的主要M文件。

原始数据文件格式设计

水准网平差所需的数据需要从txt文本中提取，本节主要介绍原始数据文件的内容与格式，需要说明的是，原始数据文件的设计并没有统一的、严格的标准，在方便程序设计的基础下，程序设计者可以自由设计数据文件的格式。

我们将网的数据分为三类：

网的概况信息、已知数据、观测数据。

网的概况信息包括总点数、已知点点数、观测者总数、验前单位权中误差。

已知数据包括已知点名、已知高程值。

观测数据包括高差起点点名、高差终点点名、观测值高差值、线路长度。

下面结合实例说明原始文件的具体数据格式

图5-2-1为典型的水准网

表5-2-1已知点数据

点名

高程/m

A

B

图5-1水准网示意图

表5-2-2观测高差

No.

起点

终点

h/m

S/km

No.

起点

终点

h/m

S/km

1

A

P1

5

P1

P2

2

A

P2

6

P1

P3

3

B

P1

7

P3

B

4

B

P2

利用以上数据进行水准网平差，数据文件的内容如下：

752

A1B1P1P2P3

A1

B1

A1P1

A1P2

B1P1

B1P2

P1P2

P1P3

P3B1

格式说明：

（1）第一行为网的概况信息：

观测值总数、总点数、已知点总数、验前单位权中误差。

（2）第二行是所有点点名依次书写。

（3）第三行是对应已知点点名、高程（单位为m），依次往下写、列出所有一直点名和高程。

（4）接下来是观测高差：

高差起始点名、高差终点点名、高差观测子和高差线路长度。

网中点名需要字母和数字结合，否则点名读不出来，点名之间空格；实际平差时，水准网规模和数据可能会不同，但是只要按照上面的格式和顺序输入相关数据即可用本章的程序。

数据的存储

1点数据的存储

平差程序用到的数组中，点名、高程值、高程改正数等数组与高程点名一一对应，这点数据称为点数据。

点数据数组的总长度为总点数。

高程值用数组Height保存，高差改正数用数据dX存放，近似值存放在数组caheight中。

点名用数组pname数组存放。

pname数组不仅存放数组名，同时也可看做点名地址数组。

点数据数组存放顺序必须保持严格的一致，例如某点点名地址存在pname数组中的第3单元，即pname（3）中，该点的高程值也就存放在Height（3）中，高程改正数也存放在dX（3）中。

当点数据存放一致时，可以更快捷的进行点数据存储，更利于程序设计。

2观测数据的存储

观测数据包括高差观测值、高差起始点的点名、高差的权、高差的改正数，它们与观测值一一对应，数组长度都等于总点数。

高差起点点名和终点点名在程序中的作用是缺点观测值的起始点号和终点点号，有了点号才能访问点数据数组，进行相应的计算。

利用函数strcmp将点名和点号一一对应，通过点号可以得到点名，当然也能通过点名得到点号，高差的起始点点号和终点点号分别用数组startp，endp来保存。

个观测数据数组的存放顺序也必须是一致，即L（i）、P（i）、V（i）、startp（i）、endp（i）均对应于同一个观测量。

近似高程计算

1近似高程值计算思路

第一步，设置未知点标志。

各点的高程值用数组Height保存，近似值高程计算之前，由于正常的高程值不可能小于，将Height数组中未知点高程赋值为，根据高程数组中的值是否小于判断某点是否需要计算近似高程，当该点近似高程计算出来之后就会被取代。

第二步，计算高程值。

近似高程计算的基本思路是，遍历观测值，找到每一个观测值起点号、终点号，根据点号从Height数组中提取起点和终点的高程值,当高差的一端是高程已知点，另一端是未知点时，就由已知点高程和观测高差算出未知点高程，并将计算结果赋值给高程数组，如果两端都是未知点，就跳到下一循环。

在遍历观测值时，排在前面的观测值总是最先用于高程近似值计算，当观测值的排序不同时，近似高程计算用到的观测值也不同，计算结果也不同。

有时候一次遍历观测值可能还有部分点的高程值仍未计算出来，还需要再次遍历观测值，直至算出所有高程值为止。

2函数文件calculateHO（）代码

functiony=calculateHO（a,b,c,d,e,f,g,h）

%计算近似高程

%caheight=calculateHO（m_Knumber,m_Pnumber,m_Lnumber,height,startP,endP,pname,L）;

jj=0;

fori=（a+1）:

b

d（1,i）=;%未知点的初始高程为

end

fori=1:

（b-a）%循环的次数不超过（m_Pnumber-m_Knumber）

forj=1:

c%按照观测组进行循环

K1=e（j）;%起点点号

K2=f（j）;%终点点号

ifd（1,K1）>&&d（1,K2）<%起点未知，终点已知

d（1,K2）=d（1,K1）+h（1,j）;

jj=jj+1;

elseifd（1,K1）<&&d（1,K2）>%起点已知，终点未知

d（1,K1）=d（1,K2）-h（1,j）;

jj=jj+1;

end

end

y=d;

ifjj==（b-a）

break

end

ifi>（b-a）%生成产生错误的txt

fid=fopen（'','a+'）;

fprintf（fid,'\n\n下列点无法计算出概略高程:

\\n'）;

fori=1:

b

ifd（i）<

fprint（fid,'\n%s',f{i}）;

end

end

fclose（fid）

errordlg（'computationFailed!

'）;

break

end

End

组成法方程

组成法方程式是参数平差的关键步骤，由参数平差模型可知，在误差方程系数矩阵

，误差方程自由项

及观测值权阵

确定时，法方程系数矩阵和自由项分别为

和

。

计算法方程的ca_ATPA（）代码如下：

function[xyz]=ca_ATPA（a,b,c,d,e,f,g,h）

%[AGPGL]=ca_ATPA（m_Pnumber,m_Lnumber,startP,endP,L,P,caheight,m_Knumber）

%A为系数矩阵，GP为权阵，GL为常系数矩阵

A=zeros（b,a）;

R=zeros（b,b）;

L=zeros（b,1）;

%获取观测值的常系数阵

fork=1:

b

i=c（1,k）;

j=d（1,k）;

L（k,1）=e（1,k）-（g（1,j）-g（1,i））;

end

%获取观测值的权阵

fori=1:

b

P（i,i）=f（1,i）;

end

%获取法方程的系数矩阵

fork=1:

b

i=c（1,k）;

j=d（1,k）;

ifi>h&&j>h

A（k,i）=-1;

A（k,j）=1;

end

ifi<=h&&j>h

A（k,j）=1;

end

ifi>h&&j<=h

A（k,i）=-1;

end

end

x=A（:

（h+1）:

a）;%除去已知点的系数，使其成为满秩矩阵

y=P;

z=L;

残差计算

残差也叫观测值的平差改正数，在参数平差中观测值的平差值可以直接用参数平差值计算出来，所以计算残差的目的不是用来改正观测值，而是主要用来进行精度估计。

参加计算由M文件ca_V完成，计算结果存在数据V中，并返回[pvv]值（用于计算单位权中误差）。

1算法

以观测值序号

为循环变量，观测值循环，由式（5-2-1）计算各观测值的残差

，第

次循环中所要进行的工作包括：

（1）获得高差的起点号

和终点号

；

（2）获得起点和终点的高程值

和

（3）计算残差

，存在

中；

（4）计算

。

2ca_V函数源代码

function[xy]=ca_V（a,b,c,d,e,f）

%计算高差改正数

%[Vpvv]=ca_V（m_Lnumber,startP,endP,Height,L,P）

pvv=;

fori=1:

a

k1=b（i）;

k2=c（i）;

V（i）=d（1,k2）-d（1,k1）-e（1,i）;

end

x=V;

y=V.*V*f';

精度估计与平差结果输出

精度估计与平差结果输出由函数printlevledata（）完成。

该函数的平差结果写入结果文件。

输出内容包括高程点名，高程平差值及其中误差、高差端点点名、高差平差值、高差改正数、高差平差值的中误差。

M文件printlevledata的源代码如下：

functionprintlevledata（pathName,Filename,m_Lnumber,m_Pnumber,m_Knumber,m_Sigma,pname,height,startP,endP,L,V,P,pvv,m_mu,caheight,dx,Height,qii,Q）

%输出结果

[FilenamepathName]=uiputfile（'*.txt','savingfileas'）;

ifisequal（Filename,0）||isequal（pathName,0）,

return;

end

str=[pathNameFilename];

%依次输入结果到txt文件

fid1=fopen（str,'wt'）;

fprintf（fid1,'TotalObserved:

%dTotalpoints:

%dTotalknownpoints:

%d\n',m_Lnumber,m_Pnumber,m_Knumber）;

fprintf（fid1,'\nPrioriunitweightweightmeanerror:

%',m_Sigma）;

fprintf（fid1,'\n\nKnownHeight:

\n'）;

fori=1:

m_Pnumber

fprintf（fid1,'\n%5d%8s',i,pname{i}）;

fprintf（fid1,'%',height（i））;

end

fprintf（fid1,'\n\nThedifferenceofObservations\n'）;

fori=1:

m_Lnumber

fprintf（fid1,'\n%5d%8s%8s',i,pname{startP（i）},pname{endP（i）}）;

fprintf（fid1,'%%',L（i）,P（i））;

end

fprintf（fid1,'\n\nReslutsofLeastSquares:

\n[pvv]=%',pvv）;

fprintf（fid1,'\nu=%1f',m_mu）;

fprintf（fid1,'\n\n\n====TheHeightAdjustmengtandPrecision====\n'）;

fprintf（fid1,'\nPnameApproHeightVAdjustHeightRMSE\n'）;

fori=1:

m_Pnumber

fprintf（fid1,'\n%6s',pname{i}）;

fprintf（fid1,'%%%%',caheight（1,i）,dx（1,i）,Height（1,i）,qii（1,i））;

end

fprintf（fid1,'\n\n\n====TheObservationsAdjustmengtandPrecision====\n'）;

fprintf（fid1,'\nNO.SPointEPointDiffObservationsV'）;

fprintf（fid1,'AdjustObservationsWObservationsRMSE\n'）;

fori=1:

m_Lnumber

k1=startP（i）;

k2=endP（i）;

ml=0;

ifk1<=m_Knumber&&k2>m_Knumber

ml=qii（1,k2）;

end

ifk1>m_Knumber&&k2<=m_Knumber

ml=qii（1,k1）;

end

ifk1>m_Knumber&&k2>m_Knumber

ml=sqrt（Q（k1-m_Knumber,k1-m_Knumber）+Q（k2-m_Knumber,k2-m_Knumber）-2*Q（k1-m_Knumber,k2-m_Knumber））*m_mu;

end

fprintf（fid1,'\n%5d%5s%5s',i,pname{k1},pname{k2}）;

fprintf（fid1,'%%%%%',L（1,i）,V（1,i）,L（1,i）+V（1,i）,P（1,i）,ml）;

end

fclose（fid1）;

ifexist（str）==2

msgbox（'savedatasuccessfully'）;

else

errordlg（'Havenotsaveddata','Failed'）;

return;

end

水准网平差算例及分析

要进行最小二乘平差计算，还需要编写一个主函数M文件，进行数据提取，连接上面提到的函数文件，才能完成特定的平差任务。

1计算步骤

（1）已知点提取处理。

（2）计算近似高程。

（3）计算法方程式。

（4）阶段法方程，最小二乘平差。

（5）计算残差。

（6）最后成果输出，包括高程平差值及其精度、高差平差值及其精度。

2函数源代码

%读取起始文件

[FilenamepathName]=uigetfile（{'*.txt','TxtFiles|*.txt'},'chooseaFile'）;%打开文件

%检查读取文件格式

[m_Lnumber,m_Pnumber,m_Knumber,m_Sigma,unpnumber,P1,P2,P3,P4,pname]=levelcheck（Filename,pathName）;

%提取已知点数据

fori=1:

m_Knumber

pname1{i}=P1{i};

point2（i）=str2num（cell2mat（P2（i）））;

end

%提取观测值及每个观测值的起终点号，路线长度

fori=1:

m_Lnumber

cpname1{i}=P1{i+m_Knumber};

cpname2{i}=P2{i+m_Knumber};

point3（i）=P3（i+m_Knumber）;

point4（i）=P4（i+m_Knumber）;

end

%按已知点点高程按着该点在点名数组中未知，将高程值赋给对应位置

height=zeros（1,m_Pnumber）;

fori=1:

m_Knumber

c=strmatch（pname1（i）,pname）;

height（c）=point2（i）;

end

%将观测值高差的起点点名转换成数组位置

fori=1:

m_Lnumber

startP（i）=strmatch（cpname1（i）,pname）;

end

%将观测值高差的终点点名转换成数组位置

fori=1:

m_Lnumber

endP（i）=strmatch（cpname2（i）,pname）;

end

%获取每条线路的高差值及其对应的权

fori=1:

m_Lnumber

L（i）=point3（i）;

P（i）=point4（i）;

P（i）=P（i）;

end

%计算近似高程

caheight=calculateHO（m_Knumber,m_Pnumber,m_Lnumber,height,startP,endP,pname,L）;

%计算法方程系数A，权阵GP，常数项GL

[AGPGL]=ca_ATPA（m_Pnumber,m_Lnumber,startP,endP,L,P,caheight,m_Knumber）;

%最小二乘平差

dX=inv（A'*GP*A）*A'*GP*GL;%未知点的改正数

dx=zeros（1,m_Pnumber）;%所以点的改正数为定义为0

%将改正数为零的未知点的改正数替换成计算值

fori=（m_Knumber+1）:

m_Pnumber

dx（1,i）=dX（i-m_Knumber,1）;

end

%计算未知点的最佳估值

fori=（m_Knumber+1）:

m_Pnumber

caheight（i）=caheight（i）+dX（i-m_Knumber,1）;

end

Height=caheight;

%计算高差改正数

[Vpvv]=ca_V（m_Lnumber,startP,endP,Height,L,P）;

%计算验后单位权中误差

m_mu=sqrt（pvv/（m_Lnumber-（m_Pnumber-m_Knumber）））;

%计算未知点的点位中误差

Q=inv（A'*GP*A）;

qii=zeros（1,m_Pnumber）;%所有点的点位中误差定义为零

fori=（1+m_Knumber）:

m_Pnumber

qii（1,i）=sqrt（Q（i-m_Knumber,i-m_Knumber））*m_mu;

end

%输出结果

printlevledata（pathName,Filename,m_Lnumber,m_Pnumber,m_Knumber,m_Sigma,pname,height,startP,endP,L,V,P,pvv,m_mu,caheight,dx,Height,qii,Q）;3最小二乘算例

水准网如图5-3-1所示，共有7个高程点，其中

、

、

、

为已知点，

、

、

为未知点，共观测了6段高差，已知高程和观测高差分别见表5-3-1和表5-3-2，观测值的每千米高差中误差为

。

试以此数据进行最小二乘平差。

表5-3-1已知点数据

点名

高程/m

A

B

C

D

图5-3-1

表5-3-2观测高差

No.

起点

终点

h/m

S/km

No.

起点

终点

h/m

S/km

1

A

E

4

F

P2

2

E

B

5

G

P2

3

E

F

6

G

P3

过程

（1）

准备数据文件

752

A1B1P1P2P3

A1

B1

A1P1

A1P2

B1P1

B1P2

P1P2

P1P3

P3B1

（2）算例函数

[FilenamepathName]=uigetfile（{'*.txt','TxtFiles|*.txt'},'chooseaFile'）;

[m_Lnumber,m_Pnumber,m_Knumber,m_Sigma,unpnumber,P1,P2,P3,P4,pname]=levelcheck（Filename,pathName）;

fori=1:

m_Knumber

pname1{i}=P1{i};

point2（i）=str2num（cell2mat（P2（i）））;

end

fori=1:

m_Lnumber

cpname1{i}=P1{i+m_Knumber};

cpname2{i}=P2{i+m_Knumber};

point3（i）=P3（i+m_Knumber）;

point4（i）=P4（i+m_Knumber）;

end

height=zeros（1,m_Pnumber）;

fori=1:

m_Knumber

c=strmatch（pname1（i）,pname）;

height（c）=point2（i）;

end

fori=1:

m_Lnumber

startP（i）=strmatch（cpname1（i）,pname）;

end

fori=1:

m_Lnumber

endP（i）=strmatch（cpname2（i）,pname）;

end

fori=1:

m_Lnumber

L（i）=point3（i）;

P（i）=point4（i）;

P（i）=P（i）;

end

caheight=calculateHO（m_Knumber,m_Pnumber,m_Lnumber,height,startP,endP,pname,L）;

[AGPGL]=ca_ATPA（m_Pnumber,m_Lnumber,startP,endP,L,P,caheight,m_Knumber）;

dX=inv（A'*GP*A）*A'*GP