MATLAB在测量控制网优化设计与平差中的应用

测绘出版社，１９９５ ［２］ 国家测绘局 ．全球定位系统（ＧＰＳ）测量规范．北 京 ：
测绘出版社，１９９２
第一作者简介 詹家民，男，测绘高级工程师，现任教于 徐州师范大学工学院，曾在《测绘通报》、《矿山测量》 《江苏测绘》等刊物上发表论文多篇。
参考文献 ［１］ 刘基余，李征航等．全球定位系统原理及应用．北京：
MATLAB 有两种常用的工作方式：一种是直接交互的
观测数据增删和观测权改变的依据。在这一过程中，其误 指令行操作方式；另一种是M 文件的编程方式。在前一种
差椭圆的绘制需要将测量坐标系化为数学坐标系，并以各 工作方式下，MATLAB 被当作一种高级的数学演算和计算
点为中心的坐标系绘制误差椭圆，其函数模型为：
幕供设计人员修改。该过程有两种：a 调整观测值；b 调
题之一。根据修改后的数据进一步进行优化，使其达到所 整观测值的权。反复调整到符合要求为止，打印输出观测
要求的技术指标。
表和设计图。
绘制控制网图与误差椭圆图是该阶段的主要任务，通
（5）MATLAB 部分源程序：
过绘制的网图进行控制网设计的调整，并将修改结果作为
第一作者简介 左廷英，女， 硕士，现在中南大学信息地理 工程学院从事教学工作，发表多篇论文。
（收稿日期：2003年1月6日）
63
况下，误差方程中的系数A 产生列亏。这种没有足够起始 数据的平差问题，称为秩亏自由网平差。
可根据最小二乘准则，组成法方程： ATPAX^ = ATPL
即 N = ATPA 则
（2）
协因数阵为： Q xx = N－ 1 = (ATPA )－ 1 （3）
由于系数阵 A 列亏 , 参数的最小二乘解就不唯一。为此
（收稿日期：２００２年１２月２日）
（上接第47页） for i1=1:1:4 j1=i1+1; s1=X(i1)/2; s2=X(j1)/2; s3=s1+s2; m1=Y(i1)/2+Y(i1+1)/2; tuoyuan(q2(j),a2(j),b2(j),m1,s3);%相对误差椭圆 end hold off
P=eye(n);
开始
L=[6275.17 6112.59 6723.85 7287.57 5485.04 6326.16];
输入初步设计或观测数据
X=[44 52.4 42.9 46.4 49.4]; Y=[13.7 18.5 25.3 19.5 12.9];
平差
计算平差值改正 数以及精度值
矩阵等的计算
第1 期 2003 年3 月
计算机应用
矿 山 测 量 MINE SURVEYING
№ .1 Mar.2003
MATLAB 在测量控制网优化设计
与平差中的应用
左廷英 刘庆元（中南大学信息物理工程学院 湖南 长沙 410083）
摘 要 文中主要阐述了MATLAB 在控制测量数据处理与控制网优化设计中的使用，主要目的是说明MATLAB 在矩阵计算与成图方面 所具有的方便、快捷的优势、以及推广使用价值。
下组成法方程
ATPA
VTPA 为奇异方阵，故该法方程是相容方程组，有
无穷多组解，故在
^
X
T
P
^
X = min
x
下求相容方程唯一的加权最小范数解，得
^
X
=
N
(NN)－ATPl
Q = N (NN )－N (NN )－N = N +
xx
两式合并则有：
其中 Px 为参数稳定程度的先验权阵，且令Px = I ，即采用
第 一 作 者 简 介 申东顺，女，工程师，1992年毕业于黑龙江矿业 学院，长期从事地质科研工作。
向，张性结构面走向NWW。一组扭性面走向NEE；另一组
（收稿日期：2002年7月30日）
（上接第39页） 数的98.63%，大部分控制在0~20 mm之间，最小 点位中误差为0。
从以上成果资料看，两个时段相比，其点位 中误差均小于10 cm的要求。说明利用RTK技术建 立图根控制，质量是可靠的，精度完全能够满足 图根控制点的精度要求，其方法是可行的。
参考文献 [1] 崔希璋，於宗俦等. 广义测量平差. 武汉：武汉测绘科技 大学出版社，2001 [2] 陶本藻. 广义逆矩阵与测量平差. 测绘工程，2000(2) [3] 周金萍，王冉，吴斌. MATLAB6实践与提高警惕. 北京：中 国电力出版社，2002. [4] 陈永奇等. 高等应用测量. 武汉测绘科技大学出版社，1996. [5] 周秋生. 测量控制网优化设计. 武汉：测绘出版社，1992
一期构造压性结构面（包括褶皱轴面和逆冲断层面） 走向 NWW；张性结构面（张性正断层）走向 NNE。两组 扭性面（平移断层），一组走向NNW（320°～330 °）； 另一组走向NEE（55°～75°）。上述平面应变图象显示地 应力的主压力近水平，方向为 NNE～SSW，主张应力为 NWW～SEE。应属于NW向构造体系。
数和权值。
行控制测量日常处理工作，特别是对一些工程控制网和变 形监测网有着重要意义。而以往编制这样的专业软件特别 是涉及矩阵求逆和矩阵运算，以及结果的可视化方面比较
2.2 平 差函数模型 平差函数模型同式（1）、式（2）相同，在处理时可分
为以下两种情况：
复杂和繁琐，而MATLAB 是一种集数值计算和可视化为一
（1）如果控制网中具有足够的起算数据，仅将其它的
体的数学软件，对数值分析，矩阵运算、信号和图形显示 待定点高程（或坐标）作为平差参数，此时误差方程中的
等数据处理具有较强功能。本文着重于用MATLAB 编制成 可供实际应用的控制网（高程网、平面网，以及GPS 网） 优化设计与平差程序。
系数阵为列满阵，由最小二乘准则可得：
附图 （1）输入初步设计数据，即将已知点和设计点坐标值 输入若是平差计算则需输入观测值数据。 （2）计算矩阵，即包括了不进行调整的控制点数据， 以及可调整点的设计数据。 若是进行控制网平差，则计算出法方程式系数阵与逆 阵以及相应点的误差椭圆。 （3）控制网优化设计中实时计算出的设计方案均需要
gtext(ff(l,1)); end juzhen(A,P,L);
矩阵的运算中又包含了矩阵的转置运算乘运算凯利逆以及在测?控制网优化设计matlab与平差中的应用左廷英刘庆元中南大学信息物?工程学院湖南长沙410083摘要文中主要阐述了在控制测?数据处?与控制网优化设计中的使用主要目的是说明在矩阵计算与成图方面matlabmatlab所具有的方?快捷的优势以及推广使用价值
点有滴水现象。可见该断层也经历了两次不同性质的构造 近水平，方向为NWW～SEE，主张应力方向为NNE～SSW,
活动。首先为水平错移，从其分支小断层可以断定其力学 应属于新华夏构造体系。
性质为顺扭，然后为上盘斜落，为张性反扭。 该方向的小断层巷道实见点较多，特征与上述的两断
层类似，多数都发育两组擦痕，断层较为平直光滑。
N
+AT1P1A1)－
1
=
(Q－xx
1+ATP1A1
)－
1
（4）
利用矩阵反演公式可得：
Q
' xx
=
Qxx－Qxx AT1
(P
－ 1
1+A1Q
xx
AT1
)－ 1A1Qxx
（5）
同理可得到删除观测值修改后的协因数阵
Q
' xx
=
Qxx+
QxxAT1
(P
－ 1
1－
A1Q xx AT1 )－ 1A1Qxx
关 键 词 MATLAB 控制网优化设计 控制网平差
1引 言 控制网优化设计与控制网平差是控制测量中较重要的两
项工作，随着计算机的日益广泛应用将两者统一为一体进
Q
' xx
=Qxx－QxxAT2
(
ΔP
－ 2
1+
A2Q
xx
AT2
)－ 1A2Qxx
（7）
其中A2 ，ΔP2 分别为改变了观测值权部分的误差方程系
中反复调用的这种通用程序段，多采用函数的形式编写， A=[0.382 0.924 0 0;-0.986 0.164 0 0;0.517 -0.856 0 0;-0.414
其总体框图如下附图。
0.910 0.414 -0.910;0 0 0.989 -0.147;0 0 -0.472 -0.881];
二期构造压性结构面（褶皱轴面和压性断层）呈NNE
NNE向断裂比较发育，力学性质为先张后压，表现为 正断层。
NNE向断层也比较发育，力学性质以扭性为主，位移 方式先水平顺扭，后张性反扭，表现为正断层。
NWW向断层较少，力学性质为先压后张，表现为正断 层。
NEE向断层较少发育，力学性质为先反扭后顺扭，位移 方式为水平移动，表现为正断层。
n=2; for j=1:1:2 tuoyuan(q1(j),a1(j),b1(j),X(j+3),Y(j+3)); end n=5; （下转第 63 页 )
47
第1 期 2003年3月
矿 山 测 量 MINE SURVEYING
№.1 Mar.2003
也为正断层，落差为1.5 m，水平擦痕极为明显，断层个别 张性兼扭性走向NNW。上述应变图象显示地应力主压应力
可视器来使用。对于简单问题，在MATLAB 的提示符下直
x = B sin T sin i + A cos T cos i
接输入命令是快速有效的。然而，当命令数量增加或希望
y = －B cos T sin i + A cos T cos i 式中
改变一个或几个变量的值或某个命令需要重复多次时，直 接输入就非常麻烦。这时可把许多可执行的MATLAB 命令
3 程序与框图 软件设计根据MATLAB 的特点，分为两大块，即优化
设计与平差。在每一块中都有些程序是通用的，故在程序
控制网和误差椭圆绘制程序 n=input('n=');% 输入观测数 t1=input('t1=');% 输入未知数 t2=input('t2=');% 输入未知数 ff=['1'; '2'; '3'; '4'; '5'];
3 构造预测 （1）褶皱轴向多数为NWW向和NNE向，多为宽缓的
由实际观测可以断定，该方向断裂也经历了两次不同 背向斜。
性质的构造活动，首先表现为压性顺扭，然后转化为张性
（2）主要大中型断层方位、力学性质、位移方式等
反扭，呈正平移断层，为一组复合式断裂。
特点如下：
2 构造体系隶属
从采区主要断裂构造的特征可知，各方位的断裂均经 历了两次不同性质的构造活动。
i —— 误差椭圆的参数变量，一般取0 °～360 °，i 放在 M 文件中，只要在 MATLAB 提示符下输入M 文件的
角增量大小取决于绘制椭圆的精度；
文件名，即可执行多个命令。以下是针对上述问题由
A、B —— 分别为椭圆的长、短半轴；
MATLAB 软件编辑器而建立的 M 文件。
T = 90 °－φ0，φ0 为长半轴方位角值。
^
X
=
(ATPA)－
1ATP
l
=N－
1ATP
l
Q = N－ 1
xx
2 控制网优化设计与平差
这就是经典的间接平差情形。 （2）如果控制网中没有足够的起算数据，而平差时需
2.1 增加或删除观测值以及修改观测值权的函数模型及算 将网中全部点的高程（或坐标）作为平差参数，在这种情
法：
^
由误差方程 V = AX －L ：P （1）
Q
' xx
=Qxxμ
QxxAT1
(P
－ 1
1±
A1Q xx AT1 )－ 1A1Qxx （6） 重心基准。
当增加观测值时取上面的符号，否则取下面的符号。
综上可知，其过程无非是求逆和矩阵运算过程。矩阵
当改变观测值权时有：
的运算中又包含了矩阵的转置运算、乘运算、凯利逆以及
46
第1 期
矿山测量
2003 年3 月
广义逆的运算等。
输出到屏幕供设计人员对网形及其精度指标检查、修改和
2.3 优化设计方案调整
调整。
不论是控制网优化还是控制网平差，都需将其网图。 绘制在屏幕上以便进行调整或以最终成果存储，而对于控 制网优化来说，在屏幕上直接进行修改，是至关重要的问
（4）实时矩阵即根据式（6）和式（7）分别计算出各 矩阵，求得 Q 'x x 的矩阵并获得各点的精度总值，输出到屏
优化 设计
计算实 时矩阵
global d; global q1; global Qx;
满足精度
绘图显示 设计结果
精度富余 删除观测值
或减少观测权
plot(Y,X,'OR'); axis('equal');
绘图输出、打印 观
精度 不够好
增加观测值 或增大观测权
hold on grid on
for l=1:1:t2
采用在最小二乘 VTPV = min 和加权最小范 数X^ T P
^
X = min
x
^
的准则下，来求未知参数的最佳估值 x 。即在
VTPV = min
若现增加一组观测值 L1，其权为 P1， 相 应 的 法 方 程 系
数为：
N 1 = (AT1P1A 1)
则修改后的协因数阵为：
Q
' xx
=
(