基于MATLAB的水准网平差程序设计
ｒｏｗ，ｃｏｉｌ＝ｓｉｚｅ（１ｉｎｅＴａｂ）Ｉ
ｆｏｒ ｉ＝１：ｒｏｗ
ｌｂｎａｍｅ＝ｌｉｎｅＴａｂ（ｉ，２）；
ｌｅｎａｍｅ＝ｌｉｎｅＴａｂ（ｉ，３）；
ｌｈｈ＝ｌｉｎｅＴａｂ（ｉ，４）；
ｌｂｉｎｄｅｘ＝ｎａｍｅＴｏｉｎｄｅｘ（１ｂｎａｍｅ，ｐｔＴａｂ）； ｌｅｉｎｄｅｘ＝ｍｍｅＴｏｉｎｄｅｘ（１ｅｎａｍｅ，ｐｔＴａｂ）； ｌｂＨ＝ｐｔＴａｂ（１ｂｉｎｄｅｘ。２）；
程序流程图如图３所示．
图３ 概略高程计算流程图
２．２
Ｆｉｇ．３ Ｈｏｗ ｃｈａｒｔ ｏｆ ｃａｌｃｕｌａｔｉｎｇ ｅｌｅｖａｔｉｏｎ
误差方程系数矩阵构建程序
程序由已知点出发，通过搜索计算获得每一个
未知点的概略高程．但由于没有考虑网的整体情况，
点的误差累积比较严重，因此点的概略高程不能满 足测量应用的要求，需要利用最ｄｘ－－乘法对水准网 进行整体平差，以获取点的最或然高程值．
３１
ｅｌｓｅ
ｉｆ ｐｔＴａｂ（１ｂｉｎｄｅｘ，４）＝一Ｏ＆ｐｔＴａｂ（１ｅｉｎｄｅｘ，４）～一Ｏ
Ｂ（ｉ，ｐｔＴａｂ（１ｅｉｎｄｅｘ，４））＝１；
ｌ（ｉ，１）一（一１）＊（１ｈｈ＋ｌｂＨ）｝
ｅｌｓｅ
ｉｆ ｐｔＴａｂ（１ｂｉｎｄｅｘ，４）～＝Ｏ＆ｐｔＴａｂ（１ｅｉｎｄｅｘ，４）－－－－＝Ｏ
Ｂ（ｉ，ｐｔＴａｂ（１ｂｉｎｄｅｘ，４））＝－－１； ｌ（ｉ，１）＝ｌｂＨ－－ｌｈｈ；
以上两表数据存储在ＭＡＴＩＡＢ矩阵对象中， 保存为ｍａｔ文件，程序运行时加载处理．
４
ｏ已知点 图１水准网示意图 Ｆ皓１ Ｓｋｅｔｃｈ ｍａｐ ｏｆｌｅｖｅｌ ｍｅｔ
２程序设计 水准网的处理过程可以分为３个步骤，即概略
坐标计算、误差方程系数阵的建立以及未知数解算 与精度评定．其中未知数解算与精度评定，在ＭＡＴ＿ Ｉ．ＡＢ只需几行代码就可以完成，因此本文重点阐述 前两个问题． ２．１概略高程计算程序
摘要：设计了水准网数据结构，并根据该数据结构将水准网边和点的相对位置关系以及观测数据、已知数据存储
在ＭＡＴＬＡＢ矩阵对象中．在ＭＡＴＬＡＢ中编制程序，利用其强大的矩阵计算功能，获取点的最或然高程值．
关键词：水准网平差；ＭＡＴＬＡＢ）最小二乘法
中图分类号：Ｐ２０９
文献标识码：Ａ
ＭＡＴＩ。ＡＢ简单易学，其强大的矩阵计算功能 尤其适合于处理测量中各种复杂的数据计算问题． 本文通过在ＭＡＴＬＡＢ中编制程序解决水准网平差 问题，来探索利用ＭＡＴＩ。ＡＢ解决测量问题的一种 新思路．
计算每点的概略高程，需要从已知点开始计算， 因此程序首先要查找到一条起算边．所谓起算边就 是边的一个端点高程已知而另一端点高程未知的观 测边．
图２观测边圈
Ｆｉｇ．２ Ｓｋｅｔｃｈ ｍａｐ ｏｆ ｌｉｎｅｓ
如图２所示，设ｉ，Ｊ两点问高差为ｈ口．如歹点为
收稿日期：２００８—０５—１２ 作者简介：李建章（１９７４一），男，甘肃会宁人，讲师．
参考文献：
［１］阮沈勇，王永利，桑群芳．ＭＡＴＬＡＢ程序设计［Ｍ］．北 京：电子工业出版社，２００３．
Ｅｚ］潘雄，付宗堂．ＭＡＴＬＡＢ软件在测量平差教学中的 应用［Ｊ］．测绘工程，２００７（１）：７６—７８．
［３］聂桂根．ＭＡＴＬＡＢ在测量数据处理中的应用［ＪＪ．测绘 通报，２００１（２）：３９－４０．
ｆｕｎｃｔｉｏｎ［ｏｋ，ｐｔＴａｂ，ｌｉｎｅＴａｂ］＝ｇｅｔＨ０（ｐｔＴａｂ，ｌｉｎｅＴａｂ） ％在边表中查找起算边． ［１ｉｎｅｌｎｄｅｘ］一ｇｅｔＢｅｇｉｎ（ｐｔＴａｂ，ｌｉｎｅＴａｂ）； ｉｆ ｌｉｎｅｌｎｄｅｘ＝＝０
ｏｋ＝ｆａｌｓｅ； ｒｅｔｕｒｎ； ｅｎｄ ％由起算边开始计算． ［ｉｓｏｋ，ｐｔＴａｂ，ｌｉｎｅＴａｂ］一ｃａｃｕｌａｔｅＨ０（ｐｔＴａｂ，ｌｉｎｅＴａｂ， ｌｉｎｅｌｎｄｅｘ）； ％在点表中查询是否有概略高程未知点． ｉｆ ｎｏｔ（ｈＯｋ（ｐｔＴａｂ）） ［ｏｋ，ｐｔＴａｂ，ｌｉｎｅＴａｂ］＝ｇｅｔＨ０（ｐｔＴａｂ，ｌｉｎｅＴａｂ）； ｅｎｄ ｏｋ＝ｔｒｕｅ．
ｅｎｄ
ｅｎｄ
’
ｅｎｄ
Ｐ（ｉ，ｉ）＝１０／ｌｉｎｅＴａｂ（ｉ，５）； ｅｎｄ
３程序算例
为检验程序的正确性，特利用某工程水准测量 实例［５３进行验证．该水准网网形如图１所示，已知数 据、观测数据和网的连接关系数据保存在表１，２中．
表１点表数据
Ｔａｂ．１ Ｄａｔａ ｏｆ ｐｏｉｎｔｓ
边号起点点号终点点号
高差／ｍ
笔者参阅了利用ＭＡＴＩＡＢ处理水准网平差问 题的一些文献，发现大多数是手工计算误差方程系 数矩阵，再利用ＭＡＴＬＡＢ解算未知数，这种方式工 作量大，可靠性和工作效率比较低．有些将数据存储 在二进制文件中进行处理，由于对文件的操作比较 复杂，编程工作量也会随之大幅增加．本文首先定义 了水准网数据结构，依该数据结构将数据保存在 ＭＡＴＩ。ＡＢ矩阵对象中，然后在程序中对各矩阵对 象直接进行操作处理，降低了编程复杂度，减少了代 码数量．
第２８卷第３期 ２００９年６月
兰州交通大学学报 Ｊｏｕｒｎａｌ ｏｆ Ｌａｎｚｈｏｕ Ｊｉａｏｔｏｎｇ Ｕｎｉｖｅｒｓｉｔｙ
文章编号：１００１—４３７３（２００９）０３—００２９－０３
Ｖ０１．２８ Ｎｏ．３ Ｊｕｎｅ ２００９
基于ＭＡＴＬＡＢ的水准网平差程序设计
李建章
（兰州交通大学土木工程学院，甘肃兰州 ７３００７０）
（Ｓｃｈｏｏｌ ｏｆ Ｃｉｖｉｌ Ｅｎｇｉｎｅｅｒｉｎｇ，ｌ丑ｎｍｈｏｕ Ｊｉａｏｔｏｎｇ Ｕｎｉｖｅｒｓｉｔｙ．１矗ｎｚｈｏｕ ７３００７０，Ｃｈｉｎａ）
Ａｂｓｔｒａｃｔ：Ｏｎ ｔｈｅ ｂａｓｉｓ ｏｆ ｄａｔａ ｓｔｒｕｃｔｕｒｅ ｄｅｓｉｇｎｅｄ，ｔｈｅ ｒｅｌａｔｉｏｎ ｏｆ ｐｏｉｎｔｓ ａｎｄ ｌｉｎｅｓ ｏｆ ｌｅｖｅｌ ｎｅｔ，ｔｈｅ ｓｕｒｖｅｙｉｎｇ ｄａｔａ ａｎｄ ｔｈｅ ｋｎｏｗｎ ｄａｔａ ａｒｅ ｓｔｏｒｅｄ ｉｎ ｍａｔｒｉｘ ｏｂｊｅｃｔ ｉｎ ＭＡＴＩ，ＡＢ．Ａｎｄ ａ ｐｒｏｇｒａｍ ｉｓ ｄｅｓｉｇｎｅｄ ｉｎ ＭＡＴＩ。ＡＢ ｔｏ ｇｅｔ ｔｈｅ ｖａｌｕｅ ｏｆ ｍｏｓｔ ｐｒｏｂａｂｌｅ ｂｙ ｉｔｓ ｓｔｒｏｎｇ ａｂｉｌｉｔｙ ｏｆ ｃａｌｃｕｌａｔｉｎｇ ｍａｔｒｉｘ． Ｋｅｙ ｗｏｒｄｓ：ａｄｊｕｓｔｍｅｎｔ ｏｆ ｌｅｖｅｌ ｎｅｔ；ＭＡＴＩ。ＡＢ；ｌｅａｓｔ ｓｑｕａｒｅ ｍｅｔｈｏｄ
ｌｅＨ＝ｐｔＴａｂ（１ｅｉｎｄｅｘ，２）； ｉｆ ｐｔＴａｂ（１ｂｉｎｄｅｘ，４）＇－－＝０ ８Ｌ ｐｔＴａｂ（１ｅｉｎｄｅｘ，４）～一Ｏ
Ｂ（ｉ，ｐｔＴａｂ（１ｂｉｎｄｅｘ，４））＝一１；
Ｂ（ｉ，ｐｔＴａｂ（１ｅｉｎｄｅｘ，４））＝１； Ｉ（ｉ，１）＝（一１）＊ｌｈｈ；
万方数据
第３期
李建章：基于ＭＡ’ＦＬＡＢ的水准网平差程序设计
万方数据
３０
兰州 交通大学学报
第２８卷
未知点，ｉ点为已知点，其高程为Ｈｉ，则歹高程为Ｈ』 ＝Ｈｉ＋ｈｄ．如ｉ点为未知点，歹点为已知点，其高程 为％，则ｉ点高程为Ｈｔ＝Ｈｊ—ｈ口．
程序获取起算边后，由该边已知点出发向前传 递高程值，直至另一概略高程已知的点（不一定是已 知点）停止．然后再判断网中是否有概略高程未知 的点，如无，程序结束，反之调用函数自身．如下为计 算未知点概略高程主程序［１￣３］．
［４］武汉大学测绘学院测量平差学科组．误差理论与测量 平差基础［Ｍ］．武汉：武汉大学出版社，２００３．
Ｅ５］姚德新．土木工程测量学教程［Ｍ］．北京：中国铁道出 版社，２００３．
Ａｄｊｕｓｔｍｅｎｔ Ｐｒｏｇｒａｎｕｎｉｎｇ ｏｆ Ｌｅｖｅｌ Ｎｅｔ ｏｎ ｔｈｅ Ｂａｓｉｓ ｏｆ ＭＡＴＬＡＢ
ＬＩ Ｊｉａｎ－ｚｈａｎｇ
如图２所示，设ｉ，歹两点间高差观测值为ｈｄ，平
差值为五ｉ，改正数为可ｄ，两点最或然高程值分别为
．７‘１７ｉ，；ｆ，近似值分别为ｚ∞，ｘｊｏ，改正数分别为如，如，
观测线路里程为晶，则可得：
ｈｏ＝ｈｏ＋％一岛一疋＋（巧。一ｚ幻） （１）
％＝岛一如＋ｚ
Hale Waihona Puke （２）其中，Ｚ—Ｘｊｏ一．２７而一ｈ＃．
如果ｉ点是已知点，Ｊ是未知点，则误差方程为
本文将水准网各边点的相对位置关系、观测数 据和已知数据通过规定的数据结构保存在ＭＡＴ— ＬＡＢ矩阵对象中，程序存取数据直观简单，进一步 缩小了程序编写的工作量．在解算水准网各点概略 高程时采用了函数自身迭代法，从而使得程序能够对 各种形状的水准网进行有效处理．由于导线网、三角 网与水准网在图形构成上有着相似性，因此本文所用 方法同样也可以应用到这些数据处理问题中去．
％一岛＋ｚ．如果歹点是已知点，ｉ是未知点，则误差
方程为铂ｆ一一站＋ｚ．
如果有九个观测值，则有ｎ个这样的误差方程．
写成矩阵形式为
Ｖ＝ＢＸ＋ｌ
（３）
仇
６１ｌ
…６１ｆ
其中
砚 Ｖ—
：
●
巩
ｂ２１ Ｂ＝＝
●
：
…６２ｆ
●
●
：
：
…ｂ。
疋。
屯
Ｘ：＝
：
●
ａ：ｔ
Ｐ１１
０一
ｏ
０ Ｐｚ２
，Ｐ＝
●
●
：
：
０ Ｏ ～．．．一 ｏ；加
线路长／ｋｉｎ
利用本程序对该水准网数据进行平差处理．程 序计算结果和实际计算结果如表３所示．
１１ａｂ．３
表３计算结果对照
Ｄｉｆｆｅｒｅｎｃｅ ｂｅｔｗｅｅｎ ｔｗｏ ｒｅｓｕｌｔｓ
程序计算结果和原结果非常吻合，说明程序设 计思路是正确的．
４结束语
利用ＭＡＴＩ。ＡＢ编程解算水准网问题，其最大 的优势在于程序编制者无须再为复杂的矩阵计算而 编写大量的代码．且ＭＡＴＩ。ＡＢ语言简单易懂，稍有 编程基础的人都可以轻松的学习它、利用它来解决 实际当中的一些计算问题．
基于MATLAB的控制网平差程序设计--第六章源代码
近似坐标计算的函数-calcux0y0函数（126页）function [x0,y0]=calcux0y0(x0,y0,e,d,sid,g,f,dir,s,t,az,pn,xyknow,xyunknow,point,aa,bb,cc)%本函数的作用是计算待定点的近似坐标format short; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% time=0;prelength=length(xyknow);non_orient=0;point_angle=0;while length(xyunknow)>0%考虑的计算方法有：1.极坐标;2.前方交会；3.测边交会；4.后方交会；%5.无定向导线的两种情况：（1）已知两个点；（2）分离的已知点与方位角；基本思路：%采用循环的方法逐一对每一个未知点进行以上各种方法条件的搜索，满足后即解算。
aa0=[];bb0=[];cc0=[];%记录搜索到两条观测边但需用户给顺序的点，注意要放在while 里面。
time=time+1; % 用于统计循环次数。
way=0;for i=xyunknow %依次循环向量中的各元素%============================================================= =====%方法1.极坐标条件搜索与计算-->way=1,基本思路：找到或求出一个方位角，找出一条边。
temp1=[]; temp2=[]; temp3=[]; temp4=[]; temp5=[]; temp6=[]; temp7=[]; temp8=[];temp9=[]; temp10=[];A=[];B=[];P=[];%第一步：寻找观测条件：两种情况：一是有已知方位角；二是由两个已知点及方向观测值推出方位角。
实验一matlab完成水准网平差
实验一matlab完成水准网平差实验一 matlab完成水准网平差实验数据：水准网有2个已知点，3个未知点，7个测段。
已知点高程H1=5.016M H2=6.016h1=1.359; h2=2.009; h3=0.363; h4=1.012; h5=0.657; h6=0.238; h7=-0.595;S1=1.1 S2=1.7 S3=2.3 S4=2.7 S5=2.4 S6=1.4 S7=2.6求解（1）求个待定点高程,H5的高差中误差；3、4号点的高程中误差。
课程设计内容1、平差程序设计思路：使用间接平差法求解（1）由题意知必要观测数t=3，选取3、4、5号点高程X1、X2、X3为参数。
（2）误差方程：V1=x1v2=x2v3=x1v4=x2v5=x2-x1+h2-h1-h5v6=x3-x1v7=-x3(3)取1M 的观测高程为单位权观测，即 p=1/s；(4)求法方程：Nbbx-W=0 Nbb=b’pbW=b’pl(5)求的平差值x=Nbb^-1*W L=l+V V=bx-l （6）高差权函数式：k=-x1+x2(6)求中误差：单位权中误差δ0，协因数阵Nbb^-1.求得中误差δ2、平差程序流程代码说明：h1=1.359;h2=2.009;h3=0.363;h4=1.012;h5=0.657;h6=0.238;h7=-0.595;H1=5.016H2=6.016h=[h1 h2 h3 h4 h5 h6 h7]'s=[1.1 1.7 2.3 2.7 2.4 1.4 2.6]'B=[1 0 0 ;0 1 0;1 0 0;0 1 0 ; -1 1 0 ; -1 0 1 ;0 0 -1 ] p=diag(1./s)l=[0;0;4;3;7;2;0]W=B'*p*lNbb=B'*p*Bx=inv(Nbb)*WV=(B*x-l)H=h+V/1000Q=inv(Nbb)n=7;t=3;j=V'*p*Vd= sqrt(j/4)f=[-1 1 0]'q=f'*Q*fD=d*sqrt(q)D1=d*sqrt(Q)(3) 平差程序流程代码说明: clc cleardisp(‘观测高差，单位m’)h1=1.359;h2=2.009;h3=0.363;h4=1.012;h5=0.657;h6=0.238;h7=-0.595;H1=5.016 % 已知点高程，单位mH2=6.016 % 已知点高程，单位mh=[h1 h2 h3 h4 h5 h6 h7]'s=[1.1 1.7 2.3 2.7 2.4 1.4 2.6]' %S是线路长度disp(‘系数矩阵B、l’)B=[1 0 0 ;0 1 0;1 0 0;0 1 0 ; -1 1 0 ;-1 0 1 ;0 0 -1 ]p=diag(1./s) %定义权阵l=[0;0;4;3;7;2;0]W=B'*p*lNbb=B'*p*Bdisp(‘参数的解’)x=inv(Nbb)*WV=(B*x-l) % 误差方程(mm)。
基于MATLAB的水准网和测边网平差程序设计
基于MATLAB的水准网和测边网平差程序设计摘要MATLAB是目前在研究机构广泛应用的一种数值计算及图形工具软件，它的特点是语法结构简明、数值计算高效、图形功能完备，特别适合非专业编程员完成数值计算、科学试验处理等任务。
以往的测量数据处理方法需要编制特定的处理矩阵运算程序，而且程度复杂，难度大。
本文介绍一种基于MATLAB的水准网和测边网的程序设计方法，与其它算法语言相比，具有编程简单，运算速度快的特点。
文中分别阐述了水准网和测边网程序的理论基础、实现步骤和运行结果。
通过实例的分析，总结出利用MATLAB对测量数据处理有很大的应用价值，它缩短了编程的时间，提高工作效率。
关键词：MATLAB；水准网；测边网；程序设计ABSTRAC TMATLAB is one species of numerical-values calculation and graphic tools software which is widely used to apply at research institutions at present. The particularities are: concise grammar-structure、highly efficient in numerical values calculating、complete function of graphs、especially it is adapted to evildoing professional programmer to accomplish the tasks that are numerical-values calculating and scientific experiments treating. The ancient methods of measured data-processing need establishing special proceedings of treating matrices operation, moreover, it is complex and greatly difficult.This article introduces one programming method dealing with leveling and measuring edge network based on MATLAB. Compared with other algorithm language, it has particularities which are simply programming and quickly operating. The article separately expatiate the theories basics、realizing steps and running results at leveling and measuring edge network. With the analysis of examples, it has prodigious application value in measured data-processing by use of MATLAB. Moreover, it shortens programming time and improves working effectiveness.Key words：MATLAB；leveling network；measuring edge network；programming目录绪论 (4)1. MATLAB软件简介 (5)2．MATLAB 在测量平差中的应用 (6)2.1测量平差原理的概述 (6)2.2平差程序总体方案 (7)3.1程序的功能 (8)3.2水准模型网的间接平差 (8)3.2.1 “权”值的确定 (8)3.2.2 水准路线的平差计算 (9)3.2.3 精度评定 (11)3.3水准网间接平差程序信息设计 (11)3.4 水准网程序与使用说明 (12)3.4.1 水准网程序流程图 (12)3.4.2 水准网程序的使用 (12)3.5案例 (13)4. 测边网平差程序设计 (15)4.1数学模型 (15)4.1.1 误差方程和法方程的组成 (15)4.1.2 边长观测的权 (15)4.1.3 解算法方程 (16)4.1.4 精度评定 (19)4.2 测边网平差信息设计 (20)4.2.1 主要的技术要求 (21)4．3利用MATLAB的绘图语句绘制网图 (21)4.4测边网程序和使用说明 (22)4.5 程序代码说明： (23)4.6程序的使用算例 (25)结论 (29)致谢 (30)参考文献 (31)附录一 (32)附录二 (36)附录三 (46)绪论作为一名测量技术人员，如果不掌握一门PC机编程语言与便携计算工具，要想提高测量工作的效率几乎寸步难行。
基于MATLAB的控制网平差程序设计--第四章源代码
chkdat函数（72页）function [n1,k]=chkdat(sd,pn,n1)n=length(n1);k=0;for i=1:ni1=0;for j=1:sdif(n1(i)==pn(j))i1=1;n1(i)=j;break;endendif(i1==0)% fprintf(fit2,'%5d %5d\n',i,n1(i)k=1;endendreturnreadlevelnetdata函数（73页）function [ed,dd,sd,gd,pn,h0,k1,k2,h1,s]=readlevelnetdata global filename filepath;global ed dd sd pn gd h0 k1 k2 h1 s k11 k12;k1=[];k2=[];h=[];s=[];[filename,filepath]=uigetfile('*.txt','选择高程数据文件');fid1=fopen(strcat(filepath,filename),'rt');if(fid1==-1)msgbox('Input File or Path is not correct','Warning','warn');return;ended=fscanf(fid1,'%f',1);dd=fscanf(fid1,'%f',1);sd=ed+dd;gd=fscanf(fid1,'%f',1);pn=fscanf(fid1,'%f',sd);h0=fscanf(fid1,'%f',ed);h0(dd+1:ed+dd)=h0(1:ed);heightdiff=fscanf(fid1,'%f',[4,gd]);heightdiff=heightdiff';k1=heightdiff(:,1);%起点k2=heightdiff(:,2);%终点k11=heightdiff(:,1);%起点k12=heightdiff(:,2);%终点h1=heightdiff(:,3);%高差s=heightdiff(:,4);%距离fclose('all');%点号转换[k1,k01]=chkdat(sd,pn,k1);[k2,k02]=chkdat(sd,pn,k2);h0(1:dd)=20000;ie=0;while(1)%计算近似高程for k=1:gdi=k1(k);j=k2(k);if(h0(i)<1e4&h0(j)>1e4)h0(j)=h0(i)+h1(k);ie=ie+1;endif(h0(i)>1e4&h0(j)<1e4)h0(i)=h0(j)-h1(k);ie=ie+1;endendif(ie==dd)break;endendh0=reshape(h0,length(h0),1); returnbm1函数（75页）function id=bm1(gd,dd,k1,k2)%计算一维压缩存放的数组idid=[];for i=1:ddk=i;for j=1:gdi1=k1(j);i2=k2(j);if(i1==i&i2<k)k=i2;endif(i2==i&i1<k)k=i1;endendid(i)=k;endfor i=2:ddid(i)=id(i-1)+i-id(i)+1;endreturn一维压缩存储法方程平差（76页）global pathname filenameglobal ed dd sd pn gd h0 k1 k2 h1 s dh;p=1./s;id=bm1(gd,dd,k1,k2);mm=id(dd);a(1:mm)=0;b(1:dd)=0;for k=1:gd %形成法方程i=k1(k);j=k2(k);h1=h1(k)+h0(i)-h0(j);if(i<=dd)ii=id(i)-i;a(ii+i)=a(ii+i)+p(k);b(i)=b(i)-h1*p(k);endif(j<=dd)jj=id(j)-j;a(jj+j)=a(jj+j)+p(k);b(j)=b(j)+h1*p(k);if(i<=dd)if(i>=j)a(ii+j)=a(ii+j)-p(k);elsea(jj+i)=a(jj+i)-p(k);endendendenda=gs5(dd,a,id);%变带宽下三角紧缩存储高斯消元法dh=cy6(a,b,id,dd,1);%常数项约化与回代子程序dh(dd+1:ed+dd)=0;hm(dd+1:ed+dd)=0;for i=1:sdh(i)=h0(i)+dh(i);endvv=0;for i=1:gdL(i)=h(k2(i))-h(k1(i));v(i)=h(k2(i))-h(k1(i))-h1(i);vv=vv+v(i)*v(i)/s(i);endu=sqrt(vv/(gd-dd));for i=1:ddb(1:dd)=0;b(i)=1.0;b=cy6(a,b,id,dd,i);hm(i)=sqrt(b(i))*u;endwritelevelnetdata(pn,k1,k2,h1,v',L',h0,dh',h',hm',u);gs5函数（77页）function a=gs5(dd,a,id)%变带宽高斯消去法for i=1:ddii=id(i)-i;if(i-1==0)li=1-ii;elseli=id(ii-1)-ii+1;endfor j=li:ijj=id(j)-j;if(j-1)==0lj=1-jj;elselj=id(j-1)-jj+1;endlk=li;if(li<lj)lk=lj;endfor k=lk:j-1kk=id(k);a(ii+j)=a(ii+j)-a(ii+k)/a(kk)*a(jj+k);endendendreturncy6函数（78页）function b=cy6(a,b,id,dd,k1)%常数项约化与回代子程序for i=k1:ddii=id(i)-i;if(i==1)nd=id(i);elsend=id(i)-id(i-1);ende=0;for k=1:i-1if((i-k)<nd)e=e+a(ii+k)*b(k);endendb(i)=(b(i)-e)/(a(ii+i);endfor i=dd-1:-1:k1ii=id(i);for k=i+1:ddkk=id(k)-k;nk=id(k)-id(k-1);if(k-i<nk)b(i)=b(i)-a(kk+i)/a(ii)*b(k);endendendreturn上三角存储法方程平差程序（79页）mm=(dd+1)*dd/2;a(1:mm)=0;b(1:dd)=0;for k=1:gdi=k1(k);j=k2(k);h1=h1(k)+h0(i)-h0(j);if(i<=dd)ii=(i-1)*(dd-i/2);a(ii+i)=a(ii+i)+1/s(k);b(i)=b(i)+1./s(k)*h1;endif(j<=dd)jj=(j-1)*(dd-j/2);a(jj+j)=a(jj+j)+1/s(k);b(j)=b(j)-1./s(k)*h1;if(i<=dd)if(i<j)a(ii+j)=a(ii+j)-1/s(k);elsea(jj+i)=a(jj+i)-1/s(k);endendendenda=invsqr(a,dd);for i=1:dddh(i)=0;di=(i-1)*(dd-i/2);for j=1:dddj=(j-1)*(dd-j/2);if(j<i)dh(i)=dh(i)-a(dj+i)*b(j);elsedh(i)=dh(i)-a(di+j)*b(i);endendenddh(dd+1:ed+dd)=0;hm(dd+1:ed+dd)=0;for i=1:sdh(i)=h0(i)+dh(i);endvv=0;for i=1:gdL(i)=h(k2(i))-h(k1(i));v(i)=h(k2(i))-h(k1(i))-h1(i);vv=vv+v(i)*v(i)/s(i);enduw0=sqrt(vv/(gd-dd));for i=1:ddii=(i-1)*(dd-i/2);hm(i)=sqrt(a(ii+i))*uw0;endreturn输出数据函数（79页）function writelevelnetdata(pn,k1,k2,h1,v,L,h0,dh,h,hm,uw0)disp('待定点高程平差值及中误差:')disp('---点号----近似高程(m)-高程改正(m)-高程平差值(m)-中误差')[pn,h0,dh,h,hm]disp('高差观测值平差值:')disp('---点号------点号----观测高差(m)---高差改正(m)-平差高差(m)')[pn(k1),pn(k2),h1,v,L][filename1,pathname1]=uigetfile('*.txt','请选择输出文件');fid2=fopen(strcat(pathname1,filename1),'wt');if(fid2==-1)msgbox('Error by Opening Output File','Warning','warn');return;endfprintf(fid2,'待定点高程平差值及中误差：\n 点号--近似高程(m)--高程改正(m)-高程平差值(m)-中误差\n');fprintf(fid2,'%5d %10.4f %10.4f %10.4f %10.4f\n',[pn,h0,dh,h,hm]');fprintf(fid2,'高差观测值平差值:\n -点号---点号--观测高差(m)--高差改正(m)-平差高差(m)\n'); fprintf(fid2,'%5d %5d %10.4f %10.4f %10.4f\n',[pn(k1),pn(k2),h1,v,L]');fprintf(fid2,'单位权中误差：%10.4fm\n',uw0);% open(strcat(pathname1,filename1));fclose(fid2);return利用Matlab矩阵运算的平差程序（81页）function level3ticdisp('平差已经开始---->>>>')global ed dd sd pn gd h0 k1 k2 h1 s dh;[ed,dd,sd,gd,pn,h0,k1,k2,h1,s]=readlevelnetdata;[dh,h,V,L,uw0,uwh,uwl]=calculatelevelnet(ed,dd,sd,gd,pn,h0,k1,k2,h1,s);writelevelnetdata(pn,k1,k2,h1,V,L,h0,dh,h,uwh,uw0); %输出水准网解算结果yunxing=toc;disp(['平差过程的运行时间为',num2str(yunxing),'秒。
基于MATLAB的水准网平差程序设计与实现
ＡＰ Ａ—Ａ ｐＡｏ ｆ １ … ＋ ｌｎ１ 】 ｏ ＋Ａ ｐＡｌ－ ４ Ａ 一Ｐ－Ａ 一
Ａ Ｐ ＝Ａ ｐｌ＋Ａｆ １ ＋… ＋Ａ ｌ ｌ ｌ Ｉ ｏ ｏ ｆ Ｐ１ ｐ 一ｆ
加入法方程系数矩 阵可逆 ， 可得 ：
Ｘ＝一Ａ Ａ） Ａ ｌ （ Ｐ 一 ｐ
—
２ 最小二乘平差的计 算步骤 ． ５ （） 文件读取 已知 高程 和观测数据 ；２未知点近似高程计算 ； ） １ 从 （ ） （ 组 ３ 成法方 程式 ； ） 程系数 阵求 逆 ； ） 平差值计算 ； ） Ｖ及单 （ 法方 ４ （ 高程 ５ （ 残差 ６ 位 中误差 计算 ； ） 后成果（ 平差值 、 （最 ７ 高程 高差平 差值及它们 的中误 差） 计算及输 出。
基 孑 ＭＡ Ｌ Ｔ ＡＢ昀水准 网 蓥程序设计与实项
郑 州市规 划勘 测设 计研 究 院 陈永 星
［ 摘
王
蕾
ቤተ መጻሕፍቲ ባይዱ ；
要］ 本文首先讨论 了ＭＡ Ｌ Ｂ在测量平差 中的应 用现状 与在 国内外的研 究动 态， ＴＡ 对基 于间接平差的水准 网算 法进行 了分析 ， 在
水准网 程序设计
Ｘ＝ Ｘ。＋ ８ ｘ
２４精度估计 ． 单位权 中误差为 ：
＝
其 中 ， 为观 测值编 号 ； ｋ ｈ 是观测 高差 ； ｕ 是观测 值 的平 差改正 数 ， 叫残差 ．，表示 高差两端点 的编号（ 也 ｉ １ｊ 即点号） ｚ ３ 分别表示观 ； 、 ｊ ２ 测 高差起 点和终点 的高程平差值 ， 即平 差中的未知数 。实际平差 时还
±
Ｊ
要 引人 参数近似值 ， 设 、 为 五 、 ｝ ＿ 的近似值 ， 、 ， ｚ 妇 为平差值 与 近似值 的差 ， 也叫改 正数 ， Ｉ＝ ， ｚ ＋ 即ｚ ｚ十 而＝ ？ ， 入误差 代
基于MATLAB的控制网平差程序设计--第五章源代码 - 副本
观测数据读入程序，rddat1函数（85页）global net ed dd sd dd1 pn x0 y0 m1 m2 m3 ms pp e d sid md g f dir ni si ma s t az aa bb cc rt rr tt global pathname filenamex0=[];y0=[];e=[];d=[];sid=[];g=[];f=[];dir=[];si=[];ni=[];s=[];t=[];az=[];pn=[];[filename,pathname]=uigetfile('*.txt','请选择原始数据');fit1=fopen(strcat(pathname,filename),'rt');if(fit1==-1)msgbox('Input File or Path is not correct','Warning','warn');return;endnet=fscanf(fit1,'%d',1);[a]=fscanf(fit1,'%d',3);ed=a(1);dd=a(2);dd1=a(3);sd=ed+dd;[pn]=fscanf(fit1,'%d',sd);[a]=fscanf(fit1,'%f',2*ed);for i=1:edx0(i)=a(2*i-1);y0(i)=a(2*i);end[a]=fscanf(fit1,'%d',3);m1=a(1);m2=a(2);m3=a(3);isid=0;[a]=fscanf(fit1,'%f',2);ms=a(1);pp=a(2);[a]=fscanf(fit1,'%d %d %f',3*m1);for i=1:m1e(i)=a(3*i-2);d(i)=a(3*i-1);sid(i)=a(3*i);end[e,i1]=chkdat(sd,pn,e);[d,i2]=chkdat(sd,pn,d);i3=0;isid=i1+i2+i3;idir=0;md=fscanf(fit1,'%f',1);[a]=fscanf(fit1,'%d %d %f',3*m2);for i=1:m2n1(i)=a(3*i-2);n2(i)=a(3*i-1);unk(i)=a(3*i);end[n1,i1]=chkdat(sd,pn,n1);[n2,i2]=chkdat(sd,pn,n2);i3=0;ik=1;si(1)=1;for i=1:sdii=0;for j=1:m2if(n1(j)==j)ii=ii+1;g(ik)=n1(j);f(ik)=n2(j);dir(ik)=unk(j);ik=ik+1;endendni(i)=ii;si(i+1)=si(i)+ni(i);endidir=i1+i2+i3;iaz=0;if(m3>0)ma=fscanf(fit1,'%f',1);[a]=fscanf(fit1,'%d %d %f',3*m3);for i=1:m3s(i)=a(3*i-2);t(i)=a(3*i-1);az(i)=a(3*i);end[s,i1]=chkdat(sd,pn,s);[t,i2]=chkdat(sd,pn,t);i3=0;iaz=i1+i2+i3;endkk=isid+idir+iaz;if(kk>0)msgbox('Error by function rddat1','Warning','warn');return;endfclose('all');open(strcat(pathname,net_name,b_datafile));return误差方程与法方程的组成函数-obnorm函数（90页）function obnormglobal ed dd dd1 ni si e d g f s tglobal m1 m2 m3 ms pp md ma x0 y0 sid dir az c fit1 fit2 global a q1 pa3 qls wlo=2062.648062470964;m=m1+m2+m3;n=2*dd;sum=n*(n+1)/2.0;sd=ed+dd;a(1:m,1:9)=0.0;for i=1:sdii=4*(ni(i)+1);pa3(i,1:ii)=0.0;endc(1:sum)=0.0;w(1:n)=0.0;for i=1:m1 %边长观测误差方程dx=x0(d(i))-x0(e(i));dy=y0(d(i))-y0(e(i));ss=sqrt(dx*dx+dy*dy);cosa=dx/ss;sina=dy/ss;a(i,1)=2*e(i)-1-2*ed+1.0e-9;a(i,2)=-cosa;a(i,3)=a(i,1)+1;a(i,4)=-sina;a(i,5)=2*d(i)-1-2*ed+1.0e-9;a(i,6)=cosa;a(i,7)=a(i,5)+1;a(i,8)=sina;a(i,9)=100.0*(ss-sid(i));q1(i)=(ms^2+(ss*pp*0.0001)^2);endq1(m1+1:m2+m1)=md*md;for i=1:sdif(ni(i)==0)continue;endjj=5;z0=0.0;zal=0;for j=si(i):s(i)+ni(i)-1dx=x0(f(j))-x0(g(j));dy=y0(f(j))-y0(g(j));a0=alfa(dx,dy);z1=a0-dir(j);if(z1<=0.0)z1=z1+2.0*pi;endzal=zal+1./ql(m1+j);z0=z0+z1/q1(m1+j);endz0=z0/zal;for j=si(i):si(i)+ni(i)-1dx=x0(f(j))-x0(g(j));dy=y0(f(j))-y0(g(j));ss=dx*dx+dy*dy;a0=alfa(dx,dy);ai=-dy/ss*lo;bi=dx/ss*lo;ii=m1+j;a(ii,1)=2*g(j)-1-2*ed+1.0e-9;a(ii,2)=-ai;a(ii,3)=a(ii,1)+1;a(ii,4)=-bi;a(ii,5)=2*f(j)-1-2*ed+1.0e-9;a(ii,6)=ai;a(ii,7)=a(ii,5)+1;a(ii,8)=bi;ss=dir(j)+z0;if(ss>=2.0*pi)ss=ss-2.0*pi;enda(ii,9)=(a0-ss)*lo*100.0;pa3(i,jj)=a(ii,5);pa3(i,jj+1)=a(ii,6)/q1(ii);pa3(i,jj+2)=a(ii,7);pa3(i,jj+3)=a(ii,8)/q1(ii);pa3(i,2)=pa3(i,2)+a(ii,2)/q1(ii);pa3(i,4)=pa3(i,4)+a(ii,4)/q1(ii);jj=jj+4;endpa3(i,1)=a(ii,1);pa3(i,3)=a(ii,3);qls(i)=-zal;endfor i=1:m3dx=x0(t(i))-x0(s(i));dy=y0(t(i))-y0(s(i));a0=alfa(dx,dy,a0);ss=dx*dx+dy*dy;ai=-dy/ss*lo;bi=dx/ss*lo;ii=m1+m2+i;a(ii,1)=2*s(i)-1-2*ed+1.0e-9;a(ii,2)=-ai;a(ii,3)=a(ii,1)+1;a(ii,4)=-bi;a(ii,5)=2*t(i)-1-2*ed+1.0e-9;a(ii,6)=ai;a(ii,7)=a(ii,5)+1;a(ii,8)=bi;if((a0-az(i))>pi)a(ii,9)=(a0-az(i)-2.0*p)*lo*100;elsea(ii,9)=(a0-az(i))*lo*100.0;endq1(ii)=ma*ma;endfor i=1:m %形成法方程for j=1:4jj=fix(a(i,2*j-1));if(jj<=0)continue;endw(jj)=w(jj)+a(i,2*j)*a(i,9)/q1(i);di=(jj-1)*(n-jj/2.0);for k=1:4kk=fix(a(i,2*k-1));if(kk<=0|jj>kk)continue;endc(di+kk)=c(di+kk)+a(i,2*k)*a(i,2*j)/q1(i);endendendif(m2>0) %和误差方程形成法方程for i=1:sdif(ni(i)==0)continue;endfor j=1:2*(ni(i)+1)jj=fix(pa3(i,2*j-1));if(jj<=0)continue;enddi=(jj-1)*(n-jj/2);for k=1:2*(ni(i)+1)kk=fix(pa3(i,2*k-1));if(kk<=0|jj>kk)continue;endc(di+kk)=c(di+kk)+pa3(i,2*k)*pa3(i,2*j)/qls(i);endendendendreturn平差值与精度评定（94页）global net ed dd sd dd1 pn x0 y0 m1 m2 m3 ms pp e d sid md g f dir ni si ma s t az global aa bb cc rt rr ttglobal a q1 pa3 qls w c x y uw0global pathname net_name s_datafile a_datafile;fit2=fopen(strcat(pathname,net_name,a_datafile,'wt');if(fit2==-1)msgbox('Input File or Path is not correct','Warning','warn');return;endk=1;while(k)m=m1+m2+m3;obnorm;c=invsqr(c,2*dd);[uw0,k]=adjxy(fit2);endellipse(uw0,fit2);n=2*dd;sum=n*(n+1)/2.0;n1=2*(ed+dd);sum1=n1*(n1+1)/2.0;for i=sum1:-1:sum1-sum+1c(i)=c(i-(sum1-sum));endfor i=1:2*eddi=(i-1)*(n1-i/2.0);for j=i:n1if(j==i)c(di+j)=0.00000001;elsec(di+j)=0.0;endendendif(m1>0)adjs(uw0,fit2);endif(m2>0)adjd(uw0,fit2);endif(m3>0)adja(uw0,fit2);endfclose(fit2);open(strcat(pathname,net_name,a_datafile));坐标改正数计算及单位权中误差计算函数-adjxy函数（95页）function [uw0,k]=adjxy(fit2)global ed dd dd1 ni si e d g f s t pn x yglobal m1 m2 m3 ms md ma x0 y0 sid dir az cglobal a q1 pa3 qls wsd=ed+dd;n=2*dd;k=0;for i=1:ndxy(i)=0.0;di=(i-1)*(n-i/2.0);for j=1:ndj=(j-1)*(n-j/2.0);if(j<i)dxy(i)=dxy(i)-c(dj+i)*w(j);elsedxy(i)=dxy(i)-c(di+j)*w(j);endendif(abs(dxy(i))>1.0)k=1;enddxy(i)=dxy(i)/100.0;endx(1:ed)=x0(1:ed);y(1:ed)=y0(1:ed);for i=1:ddx(ed+i)=x0(ed+i)+dxy(2*i-1);y(ed+i)=y0(ed+i)+dxy(2*i);endx0(1:sd)=x(1:sd);y0(1:sd)=y(1:sd);for i=1:sdif(i<=ed)vx(i)=0.0;vy(i)=0.0;elsevx(i)=dxy(2*(i-ed)-1);vy(i)=dxy(2*(i-ed));endendfprintf(fit2,' adjusted coordinates\n');fprintf(fit2,' pn vx x vy y\n');for i=1:sdfprintf(fit2,' %6d %8.4f %14.4f %8.4f %14.4f\n',pn(i),vx(i),x(i),vy(i),y(i));endpvv=0.0;for i=1:npvv=pvv+w(i)*dxy(i)*100.0;endm=m1+m2+m3;for i=1:mpvv=pvv+a(i,9)*a(i,9)/ql(i);endif(m2>0)ii=0;for i=1:sdif(ni(i)~=0;ii=ii+1;endendenduw0=sqrt(pvv/((m-n-ii)*1.0e0));return计算各点误差椭圆-ellipse函数（98页）function ellipse(uw0,fit2)glosbal ed dd pn c x y ai bi fifprintf(fit2,' parameter of error ellipse\n');fprintf(fit2,' pn(i) mx my mm a b fi\n'); n=2.0*dd;maxmm=0.0;smm=0.0;for i=1:ddii=ed+i;di=(2*i-2)*(n-(2*i-1)/2.0);dj=(2*i-1)*(n-i);q1=c(di+2*i-1);q2=c(dj+2*i);q3=c(di+2*i);d1=sqrt(abs(q1+q2-q3));d=uw0*dl;xx=q1-q2;yy=2*q3;zz=q1+q2;mx1=sqrt(q1);my1=sqrt(q2);mm1=sqrt(zz);if(abs(xx)<1d-10)fi(i)=sign(90.0,q3);elsefi(i)=atan(yy/xx)*57.2958;endif(xx>=0&yy>=0)fi(i)=fi(i)/2.0;elseif(xx>=0&yy<=0)fi(i)=(fi(i)+360)/2.0;elseif(xx<0)fi(i)=(fi(i)+180)/2.0;endww=sqrt(xx*xx+yy*yy);a1=sqrt((zz+ww)/2.0);b1=sqrt((zz-ww)/2.0);ab1=a1-b1;mx=uw0*mx1;my=uw0*my1;mm=uw0*mm1;if(mm>maxmm)maxmm=mm;i1=ii;endsmm=smm+mm;ai(i)=uw0*a1;bi(i)=uw0*b1;ab=ai(i)-bi(i);fprintf(fit2,' %10d %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f\n',pn(ii),mx,my,mm,ai(i),bi(i),fi(i));endsmm=smm/dd;fprintf(fit2,' mse of unit weight= %9.6f\n',uw0);fprintf(fit2,' the maximum station error mm= %8.3f(cm) pn= %4d\n',maxmm,pn(i1));fprintf(fit2,' the average station error mm= %8.3f\n',smm);return边长观测值平差值改正数及精度评定-adjs函数（101页）function adjs(uw,fit2)global ed dd sd pn m1 e d sid x yfprintf(fit2,'adjusted sides\n');fprintf(fit2,'i pn(e) pn(d) side vs(cm) side+vs ms\n');for i=1:m1dx=x(d(i))-x(e(i));dy=y(d(i))-y(e(i));if(e(i)<d(i))[maa,mss]=trel(uw,e(i),d(i));else[maa,mss]=trel(uw,d(i),e(i));endss=dx*dx+dy*dy;ss=sqrt(ss);vs=(ss-sid(i))*100;sid1=sid(i)+vs/100;fprintf(fit2,' %3d %8d %8d %15.4f %10.4f %15.4f %8.2f\n',i,pn(e(i)),pn(d(i)),sid(i),vs,sid1,mss);endreturn方向观测值的平差与精度评定-adjd函数（102页）function adjd(uw,fit2)global ed dd dd1 sd pn ni si e d g f s t netglobal m1 m2 m3 ms pp md ma x0 y0 x y sid dir az cglobal a q1 pa3 qls wfprintf(fit2,'adusted directions and their accuracy\n');fprintf(fit2,' i pn(g) pn(f) dir vd(") dir+vd\n');lo=206264.8062470964;for j=1:m2q1(j+m1)=md*md;endfor i=1:sdzi=0.0;zal=0.0;for j=si(i):si(i)+ni(i)-1dx=x(f(j))-x(g(j));dy=y(f(j))-y(g(j));a0=alfa(dx,dy);z1=a0-dir(j);if(z1<0.0)z1=z1+2.0*pi;endzal=zal+1./ql(m1+j);zi=zi+z1/ql(m1+j);endif(ni(i)~=0)zi=zi/zal;endfor j=si(i):si(i)+ni(i)-1dx=x(f(j))-x(g(j));dy=y(f(j))-y(g(j));a0=alfa(dx,dy);ss=dir(j)+zi;if(ss>=2.0*pi)ss=ss-2*pi;endvd=(a0-ss)*lo;dir1=dir(j)+vd/lo;dir1=rad_dms(dir1);dir(j)=rad_dms(dir(j));fprintf(fit2,' %3d %8d %8d %14.5f %9.2f %16.5f\n', j,pn(g(j)),pn(f(j)),dir(j),vd,dir1);endendfprintf(fit2,'adusted directions\n');fprintf(fit2,' i pn(g) pn(f) dir az ma\n');for i=1:sdfor j=si(i):si(i)+ni(i)-1dx=x(f(j))-x(g(j));dy=y(f(j))-y(g(j));if(f(j)<g(j))[maa,mss]=trel(uw,f(j),g(j));else[maa,mss]=trel(uw,g(j),f(j));enda0=alfa(dx,dy);if(j==si(i))a00=a0;endss=a0-a00;if(ss<0)ss=ss+2.0*pi;endss=rad_dms(ss);a0=rad_dms(a0);fprintf(fit2,' %3d %10d %10d %14.5d %14.5d %10.4f\n',j,pn(g(j)),pn(f(j)),ss,a0,maa);endendreturn方位角观测值改正数与精度评定-adja函数（104页）function adja(uw,fit2)global m3 x y s t azl0=206264.8062470964;fprintf(fit2,'adjusted azimath \n i pn(s) pn(t) az va(") az+va ma');for i=1:m3dx=x(t(i))-x(s(i));dy=y(t(i))-y(s(i));a0=alfa(dx,dy);if(t(i)<s(i))[maa,mss]=trel(uw,t(i),s(i));else[maa,mss]=trel(uw,s(i),t(i));endif((a0-az(i))>pi)va=(a0-az(i)-2*pi)*lo;elseva=(a0-az(i))*lo;endaz1=az(i)+va/lo;az1=wg(az1);az(i)=wg(az(i));fprintf(fit2,' %3d %8d %8d %13.5f %8.2f %12.5f %8.2f',i,pn(s(i)),pn(t(i)),az(i),va,az1,maa);endreturn调用的的trel函数function [maa,mss]=trel(uw,rr,tt)global ed dd x0 y0 pn cn=2*(dd+ed);dx=x0(tt)-x0(rr);dy=y0(tt)-y0(rr);ss=sqrt(dx*dx+dy*dy);a0=alfa(dx,dy);aij=2062.648*sin(a0)/ss;bij=-2062.648*cos(a0)/ss;b=a+1;f=2*tt-1;g=f+1;da=(a-1)*(n-a/2.0);db=(b-1)*(n-b/2.0);df=(f-1)*(n-f/2.0);dg=(g-1)*(n-g/2.0);q1=c(da+a)+c(df+f)-2.0*c(da+f);q2=c(db+b)+c(dg+g)-2.0*c(db+g);q3=c(da+b)+c(df+g)-c(da+g)-c(db+f);x=q1-q2;y=2.0*q3;z=q1+q2;qs=q1*cos(a0)^2+q2*sin(a0)^2+q3*sin(2.0*a0); qa=q1*aij*aij+q2*bij*bij+2.0*q3*aij*bij;ms=sqrt(qs)*uw;maa=sqrt(qa)*uw;return。
matlab在测量控制网优化设计与平差中的应用
Ｎ －ＳＥ Ｗ Ｅ ．应属于 Ｗ向 Ｗ Ｎ 构造休系． 二期构造压性结构而 （ 褶皱轴面和压性断层） Ｎ 呈Ｎ Ｅ 向， 张性结构Ｉ走向Ｎ ｂ ｌ ＷＷ。 组扭性面走向Ｎ Ｅ Ｅ ；另一组
ｘ二Ｂｓｎ ｎｉ ｉ ＋Ａｃ ｓ ｏ ｉ ｉ ７ｓ ｏ Ｔｃ ｓ
Ｙ二 一 ｃｓ ｎ ＋ ｃｓ ｃｓ ｏ ７ ｉ ｉ Ａ ｏ ７ ｏ ｉ Ｂ ｓ 式， １ ，
ｉ － 误差椭圆的参数变觉， 一 般取０．＂３０。，ｉ －６ 角增量大小取决于绘制椭圆的粕度： Ａ Ｂ— 分别为椭圆的长、短半轴； ．
了 ９ 小 ， 中 为 长半轴方位角值 。 ＝ ０。一 ． ３ 粗俘与框图
Ｍ ＴＡ 有两种常用的＿作方式：一种是直接交互的 ＡＬＢ ［ 指令行操作方式；另 一 种是Ｍ文件的编程方式。在前一种 工作方式下，Ｍ Ｔ Ａ Ａ Ｌ Ｂ被当作一种高级的数学演算和计算 可视器来使用。对于简单问题，在ＭＡ Ｌ Ｂ的提示符下直 ＴＡ 接输入命令是快速有效的。然而，当命令数量增加或希望 改变一个或几个变量的值或某个命令需要重复多次时，直 接输入就非常麻烦。这时可把许多可执行的ＭＡ Ｌ Ｂ命令 ＴＡ 放在Ｍ文件中，只要在ＭＡ Ｌ Ｂ提示符下愉入Ｍ文件的 ＴＡ 文件名，即可执行多个命令。以下是针对＿述问题由 Ｌ ＭＡ Ｌ Ｂ软件编辑器而建立的Ｍ文件． ＴＡ
由干系敬阵 Ａ列亏 ．位的．小二乘解就不唯一 为此
ＡＰＸ ＡＰ Ｔ ＝Ｔ Ａ Ｌ 即Ｎ Ａ Ａ 二Ｔ 则 Ｐ
协因数阵为：
（） ２
采 在 小乘 －ｉ 加最 范 气 一ｉ 用最 二 内Ｖ 和 权 小 批ＴＸｍ ｍ ｎ ｎ
的准则下，来求未知参致的最佳估值ｘ． 即在 Ｑ ＝ Ｉ（Ｐ）１ ， Ｎ Ａ Ａ一 ｕ ＝ Ｔ （ ３ ） Ｖ Ｖ ｍｎ Ｐ＝ ｉ 若现增加一组观测值Ｌ 二其权为Ｐ．相应的法方程系数 ｉ 下组成法方程
通用水准网形秩亏自由网和拟稳平差程序设计
MATLAB设计任何网形的秩亏自由网平差和拟稳平差程序设计长安大学王省超2015.5程序介绍：程序适合于任何网形的水准网平差，原始数据输入到连个excel表格程序界面：原始数据录入表格：（1）DH表（2）GXLB表程序代码：function xsz=XSZ(num1,num2)%函数功能提取误差系数A[m1,n1]=size(num1);[m2,n2]=size(num2);n=0;for i=1:m1 %用来判断参数个数if num1(i,2)==1n=n+0;elsen=n+1;endendxsz=zeros(m2,n); %建立系数阵，全为零for i=1:m2 %提取系数阵q=num2(i,1);z=num2(i,2);xsz(i,z)=1;xsz(i,q)=-1;endend%---------------常数项L-------------------------------------------------------------- function l=L(num1,num2)[m1,n1]=size(num1);[m2,n2]=size(num2);l=zeros(m2,1);for i=1:m2 %计算lq=num2(i,1);z=num2(i,2);l(i,1)=num2(i,4)-num1(z,4)+num1(q,4);endl=l*1000; %把l从米换算为毫米end%-------------求平差权阵P------------------------------------------------------------ function p=P1(num1,num2)[m1,n1]=size(num1);[m2,n2]=size(num2);p=zeros(m2,m2);for i=1:m2p(i,i)=1/num2(i,5);endend%------------秩亏自由网平差--------------------------------------------------------------------- function[v,Qxx]=PTPCclear allclcglobal H1;global K1;[num1]=xlsread('DH');K=num1(:,4);K1=K';[m1,n1]=size(num1); %用来判断num1得行列[num2]=xlsread('GXLB');[m2,n2]=size(num2);A=XSZ(num1,num2) %误差方程系数r=rank(A); %矩阵的秩d=1; %秩亏数l=L(num1,num2) %误差方程常数项单位mmP=P1(num1,num2); %权N=A'*P*A %法方程系数W=A'*P*lNN=N*NN0=NN(1:r,1:r);NN_=blkdiag(inv(N0),zeros(d,d));%广义逆求逆Nm=N*NN_x=Nm*W%-----------高程平差值------H=x/1000+num1(:,4)%单位统一为MH1=H';%----精度评定-------v=A*x-l %单位mmQxx=N*NN_*N*NN_*Nmsgbox('秩亏自由网普通平差完成')%---------------拟稳平差--------------------------------------------------------------------------------------- function[v,Qxx]= NWPCclear allglobal H2;global K1;[num1]=xlsread('DH');K=num1(:,4);K1=K';;[m1,n1]=size(num1); %用来判断num1得行列[num2]=xlsread('GXLB');[m2,n2]=size(num2);A=XSZ(num1,num2) %误差方程系数r=rank(A); %矩阵的秩d=m1-r; %秩亏数l=L(num1,num2) %误差方程常数项单位mmP=P1(num1,num2); %权for i=1:m1 %找出稳定点if num1(i,7)==1f(i)=1;elsef(i)=0;endendc1=find(f==1);%找出稳定点的下标c0=find(f==0);%非稳定点下标A1=A(:,c0)A2=A(:,c1)N11=A1'*P*A1;N12=A1'*P*A2;N21=N12';N22=A2'*P*A2;M=N22-N21*inv(N11)*N12;r=rank(M);d=1;MM=M*M;M0=MM(1:r,1:r);MM_=blkdiag(inv(M0),zeros(d,d));%广义逆求逆Mm_=M*MM_;a=(A2'-N21*inv(N11)*A1')';a_=Mm_*a';b_=inv(N11)*(A1'-N12*a_);%-------结算未知数-------------------------------------------------------------x2=a_*P*l;x1=b_*P*l;x=[x1;x2]H=x/1000+num1(:,4)%单位统一为MH2=H';%------计算改正数---------------------------------------------------------------------------v=A*x-l%------计算位置参数协因数矩阵-------Qxx=[b_*P*b_' b_*P*a_';a_*P*b_' a_*P*a_']msgbox('秩亏自由网拟稳平差完成')-%-------程序界面代码----------------------------------------------------------------------------------------- function varargout = GUI(varargin)gui_Singleton = 1;gui_State = struct('gui_Name', mfilename, ...'gui_Singleton', gui_Singleton, ...'gui_OpeningFcn', @GUI_OpeningFcn, ...'gui_OutputFcn', @GUI_OutputFcn, ...'gui_LayoutFcn', [] , ...'gui_Callback', []);if nargin && ischar(varargin{1})gui_State.gui_Callback = str2func(varargin{1});endif nargout[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});elsegui_mainfcn(gui_State, varargin{:});endfunction GUI_OpeningFcn(hObject, eventdata, handles, varargin)handles.output = hObject;% Update handles structureguidata(hObject, handles);% --- Outputs from this function are returned to the command line. function varargout = GUI_OutputFcn(hObject, eventdata, handles)varargout{1} = handles.output;function radiobutton1_Callback(hObject, eventdata, handles)[v,Qxx]=PTPC;set(handles.uitable1,'Data',v);set(handles.uitable2,'Data',Qxx);global H1;set(handles.edit4,'string',num2str(H1));global K1;set(handles.edit3,'string',num2str(K1));function radiobutton2_Callback(hObject, eventdata, handles)[v,Qxx]=NWPC;set(handles.uitable1,'Data',v);set(handles.uitable2,'Data',Qxx);global H2;set(handles.edit4,'string',num2str(H2));global K1;set(handles.edit3,'string',num2str(K1));function edit1_Callback(hObject, eventdata, handles)function edit1_CreateFcn(hObject, eventdata, handles)if ispc && isequal(get(hObject,'BackgroundColor'),get(0,'defaultUicontrolBackgroundColor'))set(hObject,'BackgroundColor','white');endfunction edit2_Callback(hObject, eventdata, handles)% --- Executes during object creation, after setting all properties. function edit2_CreateFcn(hObject, eventdata, handles)if ispc && isequal(get(hObject,'BackgroundColor'),get(0,'defaultUicontrolBackgroundColor'))set(hObject,'BackgroundColor','white');end% --- Executes when entered data in editable cell(s) in uitable1. function uitable1_CellEditCallback(hObject, eventdata, handles)--------------------------------------------------------------------function Untitled_1_Callback(hObject, eventdata, handles)--------------------------------------------------------------------function Untitled_2_Callback(hObject, eventdata, handles)h=msgbox({'´Ë³ÌÐòÔ-ʼÊý¾ÝµÄ¼ÈëÔÚDH(µãºÅ±í)ºÍGXLB£¨¹ØϵÁбí)ÖÐ,ÒÀ¾Ý±í¸ñ˵Ã÷½øÐÐÊý¾Ý¼È룬';'ÆäÖÐ"1"±íʾÊÇ"0"±íʾ·ñ'})% ÐÞ¸Ä×ÖÌåah = get( h, 'CurrentAxes' );ch = get( ah, 'Children' );set( ch, 'FontSize', 8 );function edit3_Callback(hObject, eventdata, handles)% --- Executes during object creation, after setting all properties. function edit3_CreateFcn(hObject, eventdata, handles)if ispc && isequal(get(hObject,'BackgroundColor'),get(0,'defaultUicontrolBackgroundColor'))set(hObject,'BackgroundColor','white');endfunction edit4_Callback(hObject, eventdata, handles)% --- Executes during object creation, after setting all properties. function edit4_CreateFcn(hObject, eventdata, handles)if ispc && isequal(get(hObject,'BackgroundColor'),get(0,'defaultUicontrolBackgroundColor'))set(hObject,'BackgroundColor','white');end% --- Executes on button press in togglebutton1.function togglebutton1_Callback(hObject, eventdata, handles)set(handles.edit3,'string','')set(handles.edit4,'string','')set(handles.uitable1,'data','')set(handles.uitable2,'data','')% --- Executes when entered data in editable cell(s) in uitable2. function uitable2_CellEditCallback(hObject, eventdata, handles)%--------------------------------------------------------------------function Untitled_3_Callback(hObject, eventdata, handles)function Untitled_4_Callback(hObject, eventdata, handles)。
实验三-利用matlab程序设计语言完成某工程导线网平差计算
实验三利用mat lab程序设计语言完成某工程导线网平差计算实验数据；某工程项目按城市测量规范(CJJ8-99)不设一个二级导线网作为首级平面控制网，主要技术要求为：平均边长200cm，测角中误差±8,导线全长相对闭合差＜1/10000,最弱点的点位中误差不得大于5cm，经过测量得到观测数据，设角度为等精度观测值、测角中误差为山=±8秒，鞭长光电测距、测距中误差为m二± Vsmm，根据所学的‘误差理论与测量平差基础'提出一个最佳的平差方案，利用matlab完成该网的严密平差级精度评定计算；平差程序设计思路：1采用间接平差方法，12个点的坐标的平差值作为参数.利用matlab进行坐标反算，求出已知坐标方位角；根据已知图形各观测方向方位角；2计算各待定点的近似坐标，然后反算出近似方位角，近似边. 计算各边坐标方位角改正数系数；3确定角和边的权,角度权Pj=1 ；边长权Ps=100/S；4计算角度和边长的误差方程系数和常数项，列出误差方程系数矩阵 B，算出Nbb=B’ PB,W=B’ Pl,参数改正数 x=inv(Nbb)*W;角度和边长改正数V=Bx-l； 6建立法方程和解算x,计算坐标平差值，精度计算;程序代码以及说明:s10=;s20=;s30=;s40=;s50=;s60=;s70=;s80=;s90=;s100=;s110=;s120=;s130=;s140=; %已知点间距离Xa=；Ya二；Xb=;Yb=;Xc=;Yc=;Xd=;Yd=;Xe=;Ye=;Xf=;Yf=; %已知点坐标值a0=atand((Yb-Ya)/(Xb-Xa))+180;d0=atand((Yd-Yc)/(Xd-Xc));f0=atand((Yf-Ye)/(Xf-Xe))+360; %坐标反算方位角a1=a0+(163+45/60+4/3600)-180a2=a1+(64+58/60+37/3600)-180;a3=a2+(250+18/60+11/3600)-180;a4=a3+(103+57/60+34/3600)-180;a5=d0+(83+8/60+5/3600)+180;a6=a5+(258+54/60+18/3600)-180-360;a7=a6+(249+13/60+17/3600)-180;a8=a7+(207+32/60+34/3600)-180;a9=a8+(169+10/60+30/3600)-180;a10=a9+(98+22/60+4/3600)-180;a12=f0+(111+14/60+23/3600)-180;a13=a12+(79+20/60+18/3600)-180;a14=a13+(268+6/60+4/3600)-180;a15=a14+(180+41/60+18/3600)-180; %推算个点方位角 aa=[a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a12 a13 a14 a15]'X20=Xb+s10*cosd(a1);X30=X20+s20*cosd(a2);X40=X30+s30*cosd(a3);X50a=X40+s40*cosd(a4);X60=Xd+s50*cosd(a5);X70=X60+s60*cosd(a6);X80=X70+s70*cosd(a7);X90=X80+s80*cosd(a8);X100=X90+s90*cosd(a9);X50c=X100+s100*cosd(a10);X130二Xf+s110*cosd(a12);X140=X130+s120*cosd(a13);X150=X140+s130*cosd(a14);X50e=X150+s140*cosd(a15); %各点横坐标近似值X0=[X20 X30 X40 X60 X70 X80 X90 X100 X130 X140 X150 X50a X50c X50e]'Y20=Yb+s10*sind(a1);Y30=Y20+s20*sind(a2);Y40=Y30+s30*sind(a3);Y50a=Y40+s40*sind(a4);Y60=Yd+s50*sind(a5);Y70=Y60+s60*sind(a6);Y80=Y70+s70*sind(a7);Y90=Y80+s80*sind(a8);Y100=Y90+s90*sind(a9);Y50c=Y100+s100*sind(a10);Y130=Yf+s110*sind(a12);Y140=Y130+s120*sind(a13);Y150=Y140+s130*sind(a14);Y50e=Y150+s140*sind(a15); %个点从坐标近似值Y0=[Y20 Y30 Y40 Y60 Y70 Y80 Y90 Y100 Y130 Y140 Y150 Y50a Y50c Y50e]'P=[X0 Y0];X50=(X50a+X50c+X50e)/3Y50=(Y50a+Y50c+Y50e)/3s4二sqrt((Y40-Y50)"2+(X40-X50厂2);si二sqrt((Y100-Y50厂2+(X100-X50厂2);s14二sqrt((Y150-Y50)"2+(X150-X50厂2);A1=[cosd(a1) cosd(a2) cosd(a3) cosd(a4) cos(a5) cosd(a6) cosd(a7) cosd(a8) cosd(a9) cosd(a10) cosd(a12) cosd(a13) cosd(a14) cosd(a15)]';B11=[sind(a1) sind(a2) sind(a3) sind(a4) sin(a5) sind(a6) sind(a7) sind(a8) sind(a9) sind(a10) sind(a12) sind(a13) sind(a14) sind(a15)]';s=blkdiag(s10,s20,s30,s4,s50,s60,s70,s80,s90,s10',s110,s120,s130,s14);a=*inv(s)*B11b=*inv(s)*A1ab4=atand((Y50-Y40)/(X50-X40))+180;ab10=atand((Y50-Y100)/(X50-X100));ab14=atand((Y50-Y150)/(X50-X150))+360;m4=ab4-a3+180;m10=ab10-a9+180;m11=ab4-ab10;m15=ab14-a14+180;m16=ab10-ab14+360;m04=103+57/60+34/3600;m010=98+22/60+4/3600;m011=94+53/60+50/3600;m015=180+41/60+18/3600;m016=ab10-ab14+360;l=[0 0 0 m4-103-57/60-34/3600 0 0 0 0 0 m10-98-22/60-4/3600 m11-94-53/60-50/3600 0 0 0 m15T80-41/60T8/3600m16-103-23/60-8/3600 0 0 0 s40-s4 0 0 0 0 0 s100-s1 0 0 0 s140-s14]';e1=(abs(X20-Xb))/s10;e2=(abs(X30-X20))/s20;e3=(abs(X40-X30))/s30;e4=(abs(X50-X40))/s4;e5=(abs(X60-Xd))/s50;e6= (abs(X70-X60))/s60;e7=(abs(X80-X70))/s70;e8=(abs(X90-X80))/s80;e9=(abs(X100-X90))/s90;e10=(abs(X50-X100))/s1;e11=(abs(X130-Xf))/s110;e12=(abs(X140-X130 ))/s120;e13=(abs(X150-X140))/s130;e14=(abs(X50-X150))/s 14;e=[e1 e2 e3 e4 e5 e6 e7 e8 e9 e10 e11 e12 e13 e14]' m1=(abs(Y20-Yb))/s10;m2=(abs(Y30-Y20))/s20;m3=(abs(Y40-Y30))/s30;m4=(abs(Y50-Y40))/s4;m5=(abs(Y60-Yd))/s50;m6= (abs(Y70-Y60))/s60;m7=(abs(Y80-Y70))/s70;m8=(abs(Y90-Y80))/s80;m9=(abs(Y100-Y90))/s90;m10=(abs(Y50-Y100))/s1;m11=(abs(Y130-Yf))/s110;m12=(abs(Y140-Y130 ))/s120;m13=(abs(Y150-Y140))/s130;m14=(abs(Y50-Y150))/s 14;m=[m1 m2 m3 m4 m5 m6 m7 m8 m9 m10 m11 m12 m13 m14]' % 以上为求得误差方程系数B=[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0%系数矩阵B0 0 ]P=blkdiag(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,100/s10,100/s 20,100/s30,100/s40,100/s50,100/s60,100/s70,100/s80,100/ s90,100/s100,100/s110,100/s120,100/s130,100/s140); %定义权矩阵Nbb二B'*P*BW=B'*P*l;x=inv(Nbb)*WV=B*x-l;inv(Nbb);Y=V'*P*V;O二sqrt(Y/6)*3600 %精度评定计算结果：平差值坐标X：+003 *Qx1= Qy1= Qx2= Qy2= ……Qx15= Qy15=。
用MATLAB解决-条件平差和间接平差
L1 A
L3 C
9
clc Disp(‘条件平差示例２’) Disp(‘三角形内角观测值’) L1 = [42 12 20] L2 = [78 9 9] L3 = [59 38 40] L = [L1; L2; L3] Disp(‘将角度单位由度分秒转换为弧度’) LL = dms2rad(mat2dms(L))
设误差Δ和参数X的估计值分别为V 和 Xˆ
15
则有
V AXˆ l
为了便于计算，通常给参数估计一个充分接近的近似值 X 0
Xˆ X 0 xˆ 则误差方程表示为
V Axˆ l
其中常数项为
l L (AX 0 d)
16
由最小二乘准则，所求参数的改正数应该满足
V T PV min
目标函数对x求一阶导数，并令其为零
if(sum(LL) == pi) disp(‘检核正确’)
else
V = A'*Ka
disp(‘检核错误’) end
11 例 《误差理论与测量平差基础》P75
在下图中，Ａ、Ｂ为已知水准点，其高程为 HA=12.013m, HB = 10.013m, 可视为无误差。为 了确定点Ｃ及Ｄ点的高程，共观测了四个高差，高差 观测值及相应的水准路线的距离为：
disp(‘系数矩阵B’)
B = [1 0; 0 1; 1 0; 0 1; -1 1; -1 0]
l = [0; 0; 4; 3; 7; 2]
disp(‘C是单位权观测高差的线路公里数，S是线路长度’)
C = l*ones(1,6)
23
h6
E h3
h7 B
S = [1.1, 1.7, 2.3, 2.7, 2.4, 4.0]