CPIII平面控制测量方法及程序

第6章 CPIII控制网数据处理
当前我国客运专线的建设多采用无碴轨道技术，由于设计速度高，为保证列车在高速运行时的安全性，以及乘客的舒适度，高速客运专线的轨道必须具有高平顺性和高稳定性。

除轨道结构的合理尺寸、良好的材质和制造工艺外，轨道的高精度铺设是实现轨道初始高平顺性的关键。

而高精度CPIII控制网是无碴轨道施工的保证，并为日后运营维护提供控制基准。

6.1 CPIII控制网基础知识
CPIII控制网是沿线路布设控制无碴轨道施工的三维施工控制网，起闭于上一级的基础平面控制网（CPI）或线路控制网（CPII）。

CPIII控制网点对称布设于线路两侧，每对间距约为15m左右，控制点间的纵向间距以50～60m为宜；CPIII平面网采用自由测站后方交会进行施测，其原始观测值为测站到测点的平距与方向，每两测站间有4对公共观测点，由此构成了一个控制网点间具有强相关性、精度分布较为均匀的边角交会网。

由于采用了全新的构网方式，需要发展相应的严密数据处理方法来对CPIII平面网观测数据进行处理。

6.1.1 CPIII相关概念
（1）工程独立坐标系：为满足铁路工程建设要求采用的以任意中央子午线和高程投影面进行投影而建立的平面直角坐标系。

（2）基础框架平面控制网CP0：为满足线路平面控制测量起闭联测的要求，沿线路每50km左右建立的卫星定位测量控制网，作为全线勘测设计、施工、运营维护的坐标基准。

（3）基础平面控制网CPⅠ：在基础框架平面控制网（CP0）或国家高等级平面控制网的基础上，沿线路走向布设，按GPS静态相对定位原理建立，为线路平面控制网起闭的基准。

在勘测阶段按静态GPS相对定位原理建立。

点间距为4km左右，测量精度为GPS B级网。

（4）线路平面控制网CPⅡ：在基础平面控制网（CPⅠ）上沿线路附近布设，为勘测、施工阶段的线路平面控制和轨道控制网起闭的基准。

可用GPS静态相对定位原理测量或常规导线网测量，在勘测阶段建立。

点间距为400~800m左右，测量精度为GPS C级网或三等导线。

（5）轨道控制网CPⅢ：沿线路布设的三维控制网，起闭于基础平面控制网（CPⅠ）或线路控制网（CPⅡ），一般在线下工程施工完成后进行施测，为轨道施工和运营维护的基准。

CPⅢ网按自由设站边角交会方法测量。

点间距为纵向60m左右、横向为线路结构物宽度，测量精度为相邻点位的相对点位中误差小于1mm。

（6）精密水准测量：高速铁路无碴轨道工程施工测量中，用于测量轨道控制网CPⅢ各
标志点高程的等级水准测量，其精度介于二等和三等水准测量之间，每公里高差测量的偶然中误差为2mm/km和全中误差为4mm/km。

（7）自由设站后方交会：在靠近线中线位置架设全站仪，测量线路两侧多对轨道控制网CPⅢ的方向和距离，以确定仪器中心点的平面和高程位置。

常用于无碴轨道板的粗调、轨道的精调和轨道线形的检测。

（8）三网合一：为保证控制网的测量成果质量满足高速铁路勘测、施工、运营维护三个阶段测量的要求，适应客运专线无碴轨道铁路工程建设和运营管理的需要，要求三阶段的平面、高程控制测量必须采用统一的基准。

即勘测控制网、施工控制网、运营维护控制网均采用CP0 控制网为基础平面控制网，二等水准基点网为基础高程控制网，简称为三网合一。

6.1.2 平面控制网测量网形
（1）测站间距为120m时，CPⅢ平面控制网测量网形示意图如图6-1所示。

图6-1 测站间距为120m时CPⅢ平面控制网测量网形示意图（2）测站间距为60m时，CPⅢ平面控制网测量网形示意图如图6-2所示。

图6-2 测站间距为60m时CPⅢ平面控制网测量网形示意图（3）采用测站间距120m的标准网形测量过程中如某CPⅢ点由于障碍物被挡，可以考虑采用由测站间距120m转测站间距60m的测量网形，如图6-3所示。

图6-3 测站间距由120m转60m时的平面控制网测量网形示意图（4）在实际测量过程中，如果CPⅠ或者CPⅡ点离线路较远，可以在线路外合适位置设
置辅助点，在辅助点上架设仪器，观测临近的CPⅢ点和CPⅠ或者CPⅡ点。

此时其测量网形示意图如图6-4所示。

图6-4 辅助点联测CPⅠ或者CPⅡ点平面控制网测量网形示意图
6.1.3 CPIII控制网的特点
(1)控制点数量众多。

沿线路方向通常每公里有16对即32个控制点。

(2)精度要求高。

每个控制点与相邻5个控制点的相对点位中误差均要求小于1mm。

(3)控制的范围长。

线路有多长，控制网的长度就有多长。

(4)是一个平面位置和高程位置共点的三维控制网。

目前CPⅢ三维网平面和高程是分开测量后合并形成共点的三维网，但其使用时却是平面和高程同时使用的。

(5)控制点的位置、CPⅢ测量标志较传统控制测量有很大不同。

控制点通常设置在接触网杆上（路基部分）、防撞墙上（桥梁部分）和围岩上（隧道部分）。

CPⅢ测量标志通常由永久性的预埋件、平面测量杆、高程测量杆和精密棱镜组成。

(6)CPⅢ平面网是一个边角控制网，但其测量方法较传统边角网测量有很大差异。

传统的边角网测量仪器都是架设在控制点上进行观测，距离必须进行往返观测，但CPⅢ平面网却采用自由设站进行边角交会测量，而其距离只能进行单程观测。

(7)CPⅢ控制网测量的仪器均采用高精度和自动化程度高的电子测量仪器。

其平面网测量要求全站仪具有电子驱动、目标自动搜索和操作系统功能的测量机器人（如 Leica TCA2003和 TCRA1201、Trimble S6和S8系列全站仪等）；高程测量一律采用电子水准仪（如Trimble DiNi12、Leica DNA03等）。

(8)测站和测点均强制对中，测点标志要求具有互换性和重复安装性。

(9)图形规则对称，多余观测数多，可靠性强。

(10)是一个标准的带状控制网，其纵向精度高、横向精度略差。

6.1.4 CPIII控制网测量的一般规定
(1) CPⅢ平面控制网测量前，应确保线路两侧50m范围内CPⅡ控制点的密度达到500m～700m，否则应同精度加密CPⅡ控制点；CPⅢ高程控制网测量前，应确保线路两侧50m范围内水准点的密度达到2000m左右，否则应同精度用水准测量的方法加密水准点。

(2) CP Ⅲ平面控制网的主要技术指标应满足表6-1的规定。

表6-1 CP Ⅲ平面网的主要技术指标
注：（1）可重复性测量精度指的是控制点两次测量，其X 、Y 方向坐标差的中误差。

（2）相对点位精度指的是相邻两点间相对点位误差椭圆长短轴平方和的开根号值。

6.2 基于非线性模型的CPIII 控制网平差
CPIII 控制网的观测值有方向观测值和距离观测值，待求参数为各点的平面坐标。

其平差计算过程主要包括待求参数近似值的计算，误差方程系数矩阵的计算，观测值权的估计三大部分。

6.2.1 待求参数近似值的计算
在计算误差方程系数矩阵前，首先要获取待求参数的近似值。

由于在各测站进行观测时，没有设定起始方向,因此需从联测有已知点的测站开始，利用平面直角坐标变换模型借助相邻测站公共点进行各待求参数近似值的计算。

（6-1）式即为平面直角坐标变换模型，其中t x 和t y 为平移参数，k 为缩放系数，为旋转分量。

⎥⎦
⎤⎢⎣⎡⨯⎥⎦⎤⎢⎣⎡-⨯++⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡i i y x i i t v y x cos sin sin cos )k 1t u θθθθ（（6-1） 通常情况下，距离联测已知点的测站越远，近似值计算精度越低。

近似值的偏差会影响到待求参数估计值的精度，为减小这项误差，通常需要将上次平差结果作为待求参数近似值反复迭代计算，以获取最佳估计值。

待求参数近似值偏差越大，迭代次数越多。

近似值计算是传统CPIII 控制网平差中关键的一个步骤。

6.2.2 误差方程系数矩阵的计算
CPIII 控制网的观测值有方向观测值和距离观测值两类。

设在某测站点j 对待求点k 进行观测，其中方向观测值对应函数模型为（6-2）式，由于是非线性模型，需要进行线性化。

（6-3）式为方向观测值对应的误差方程。

)(
~~~~~~~j k j k jk j jk X X Y Y arctg L Z --=+=α （6-2） 0
00
00
0000000
)(cos sin cos sin j
jk jk jk jk
Z X X Y Y arctg L x S x S y S x S z V j k j k
jk k jk k jk j jk j jk
j jk ---+-+--+-=αααα（6-3）
（6-4）、（6-5）分别为距离观测值所对应的函数模型和误差方程。

2
~
~2~~~)()(j k j k jk Y Y X X S -+-= （6-4）
jk
j k j k k k j j jk S Y Y X X y x y x V jk jk jk jk --+-+++--=20~~020~0~0000)()(sin cos sin cos αααα （6-5）
式中jk L 、 jk S 分别为方向观测值和距离观测值，jk L ~、jk S ~为各自平差值。

j 、k 两点最或然坐标值分别为）（j X X ~j ~,、）（k X X ~k ~,，j 测站的定向角为~
j Z 。

6.2.3 观测值权的估计
CPIII 控制网有两类观测值，需利用仪器标称精度对观测值定权并确定其权比关系。

通常, 第一次平差前按经验给出的两类观测值的权P 1和P 2通常是不合理的, 也可理解为与它们
对应的单位权方差不相等, 记为201δ和202δ。

赫尔默特方差分量估计可以利用各次平差后的
两类观测值改正数的加权平方和111V P V T 和222V P V T
来重新估计201δ和202δ。

从而确定更为合理的权阵。

根据文献[3],给出方差分量估计公式
θ~
W θS = （6-6）
（6-6）式中 ⎥⎦⎤⎢⎣
⎡+-+-=--------22121221112111211111)()(2)()()()(2S N N tr N N tr n N N N N tr N N N N tr N N tr N N tr n ， []T 202201~δδθ=，[]T
22T 211T 1θV P V V P V W = 为提高估计精度，需迭代计算，直至201δ和202δ比等于1为止。

6.2.4 代码详解
%函数功能：约束平差主程序
%基于非线性模型的CPIII数据处理程序在设计初采用线性模型平差结果作为近似值的办法function [X,tQ,tmo,tQf,tmof,XY,V]=nonelinarpc(gcJZ,XY,czs,ceZB,xzJD,Q)
format long e; %设置变量类型
pi=3.1415926535898; %π变量
rou=180*360/pi; %秒和弧度互化常数
[row,col]=size(XY); %获取XY矩阵行列数
for i=1:row %在XY矩阵中循环
XY(i,15)=XY(i,11);
XY(i,16)=XY(i,12); %将线性模型平差结果作为近似值end
[row,col]=size(gcJZ); %获取观测矩阵行列数
jsq=0; %计数器置零
for i=1:row %在观测矩阵中循环
if gcJZ(i,1)~=0 %如果不是分割行
jsq=jsq+1; %计数器加1
end
end
surveysum=jsq;
maxnum=max(XY);
wznum=maxnum(5); %未知点个数
canfinish=0; %迭代标志
while canfinish==0 %如果标志变量为0，继续迭代
jsq=0; %计数器置零
czindex=0; %测站序号置零
dist=0; %观测边长置零
angle=0; %观测方位置零
index=0;
px=0;
py=0;
cpx=0;
cpy=0;
jsfw=0;
B1=0;
B2=0;
B=0;
P1=0;
P2=0;
P=0;
l1=0;
l2=0;
l=0;
for j=1:row %在观测矩阵中循环
if gcJZ(j,1)==0 %如果是分割行
czindex=gcJZ(j,2); %读取测站序号
else %如是观测值所在行
jsq=jsq+1; %计数器加1
pname=gcJZ(j,1); %获取观测点名称
wzindex=getindex(XY,pname,5); %获取当前观测点对应未知数序号 angle=dfmtohd(gcJZ(j,2)); %获取观测方位
dist=gcJZ(j,3); %获取观测距离
index=getindex(XY,pname,6); %获取点索引
px=XY(index,15);
py=XY(index,16); %获取目标点近似坐标
cpx=ceZB(czindex,1);
cpy=ceZB(czindex,2); %获取测站点坐标
dist0=sqrt((px-cpx)^2+(py-cpy)^2); %计算近似距离
jsfw=getfw((px-cpx),(py-cpy)); %计算近似方位
jj=xzJD(czindex,1); %获取测站旋转角度
while jj>=2*pi %如果角度大于2π
jj=jj-2*pi; %反复减去2π，使之小于2π
end
while jj<0 %如果旋转角度大于2π
jj=jj+2*pi; %反复加2π
end
det=0;
if jsfw>jj %如果近似方位大于夹角
k=-1;
det=0; %补偿值为0
end
if jsfw<=jj %如果近似方位小于等于夹角
k=-1;
det=2*pi; %补偿值为2π
end
if wzindex~=0 %如果未知数序号不为0，即未知点 B1(jsq,2*wzindex-1)=-rou*sin(jsfw)/dist0/1000;
B1(jsq,2*wzindex)=rou*cos(jsfw)/dist0/1000;
B1(jsq,2*wznum+2*czindex-1)=rou*sin(jsfw)/dist0/1000;
B1(jsq,2*wznum+2*czindex)=-rou*cos(jsfw)/dist0/1000;
%方位观测值误差方程中测站坐标之系数
B1(jsq,2*wznum+2*czs+czindex)=k; %起始方位旋转夹角 l1(jsq,1)=-rou*(det+jsfw-angle+jj*k); %常数项
P1(jsq,jsq)=1; %方位权
B2(jsq,2*wzindex-1)=cos(jsfw);
B2(jsq,2*wzindex)=sin(jsfw);
B2(jsq,2*wznum+2*czindex-1)=-cos(jsfw);
B2(jsq,2*wznum+2*czindex)=-sin(jsfw);
%距离观测值误差方程中测站坐标之系数
B2(jsq,2*wznum+2*czs+czindex)=0; %起始方位夹角 l2(jsq,1)=-1000*(dist0-dist); %常数项
P2(jsq,jsq)=1/(1+1*dist0/1000)^2; %距离观测权 else %如果目标点是已知点
B1(jsq,2*wznum+2*czindex-1)=rou*sin(jsfw)/dist0/1000;
B1(jsq,2*wznum+2*czindex)=-rou*cos(jsfw)/dist0/1000;
%方位观测值误差方程中测站坐标之系数
B1(jsq,2*wznum+2*czs+czindex)=k;
l1(jsq,1)=-rou*(det+jsfw-angle+jj*k);
P1(jsq,jsq)=1;
B2(jsq,2*wznum+2*czindex-1)=-cos(jsfw);
B2(jsq,2*wznum+2*czindex)=-sin(jsfw);
%距离观测值误差方程中测站坐标之系数
B2(jsq,2*wznum+2*czs+czindex)=0;
l2(jsq,1)=-1000*(dist0-dist);
P2(jsq,jsq)=1/(1+1*dist0/1000)^2;
end
end
end
%以上分别获取方位观测值和距离观测值对应误差方程系数矩阵以及常数项
[temr,temc]=size(B1); %获取B1矩阵行数
for i=1:surveysum %根据B1矩阵行数循环
for j=1:temc
B(i,j)=B1(i,j);
B(i+surveysum,j)=B2(i,j);
end
end %由B1和B2矩阵合并成B矩阵
for i=1:surveysum
l(i,1)=l1(i,1);
l(i+surveysum,1)=l2(i,1);
end %由l1和l2矩阵合并成l矩阵
for i=1:surveysum
for j=1:surveysum
if i==j
P(i,j)=P1(i,j);
P(i+surveysum,j+surveysum)=P2(i,j);
end
end
end %由由P1和P2矩阵合并成P矩阵 X=inv(B'*P*B)*B'*P*l; %解算未知参数矩阵
[row1,col1]=size(XY); %获取XY矩阵行列数
for k=1:row1 %在XY矩阵中循环
wzindex=XY(k,5); %获取第k个点未知数序号
if wzindex~=0 %如果为未知点
XY(k,15)=X(2*wzindex-1,1)/1000+XY(k,15);
XY(k,16)=X(2*wzindex,1)/1000+XY(k,16); %利用平差值对近似值加改正 end
end
[row2,col2]=size(ceZB); %获取测站坐标矩阵行列数
for k=1:row2
ceZB(k,1)=X(2*wznum+2*k-1,1)/1000+ceZB(k,1);
ceZB(k,2)=X(2*wznum+2*k,1)/1000+ceZB(k,2);
end %测站坐标进行改正
for k=1:row2
temfw=xzJD(k,1);
temfw=temfw+X(2*wznum+2*czs+k,1)/rou;
while temfw>=2*pi
temfw=temfw-2*pi;
end
while temfw<0
temfw=temfw+2*pi;
end
xzJD(k,1)=temfw;
end %对起始方位夹角做修正
for y=1:2*wznum
temPT1(y,1)=X(y,1);
end
maxx1=max(abs(temPT1)); %获取x未知数改正量最大值
for y=1:2*czs
temPT2(y,1)=X(2*wznum+y,1);
end
maxx2=max(abs(temPT2)); %获取y未知数改正量最大值
if maxx1>maxx2
maxx=maxx1;
else
maxx=maxx2;
end %获取所有未知数改正量最大值
if maxx<0.0001 %如果未知数改正量最大值小于该值 canfinish=1; %修改循环标志变量为1，即可以退出 end
end
tXf=X;
V=B*X-l; %获取残差矩阵
[n,t]=size(B);
tmof=sqrt((V'*P*V)/(n-t)); %解算单位权中误差
tQf=inv(B'*P*B); %未知数协因数矩阵
tMf=tmof^2*tQf; %未知数方差阵
%以上通过迭代计算获取无定权情况下未知参数的解
XI(1,1)=1;
XI(2,1)=1; %两类观测值单位权方差初值都取1 notfirst=1; %一循环次数变量初值取1
while notfirst==1 || abs(XI(1,1)-XI(2,1))>0.000001
P=0;
P1=XI(1,1)*inv(XI(1,1)*inv(P1));
P2=XI(1,1)*inv(XI(2,1)*inv(P2)); %获取两类观测值权矩阵
for i=1:surveysum
for j=1:surveysum
if i==j
P(i,j)=P1(i,j);
P(i+surveysum,j+surveysum)=P2(i,j);
end
end
end %获取统一的观测值权矩阵
X=inv(B'*P*B)*B'*P*l; %解算未知参数矩阵
N=B'*P*B;
V=B*X-l; %获取残差值矩阵
notfirst=2; %循环次数变量设为1
V1=B1*X-l1;
V2=B2*X-l2;
[row1,col1]=size(B1);
[row2,col2]=size(B2);
N1=B1'*P1*B1;
N2=B2'*P2*B2;
W1=V1'*P1*V1;
W2=V2'*P2*V2;
W(1,1)=W1;
W(2,1)=W2;
S(1,1)=row1-2*trace(inv(N)*N1)+trace((inv(N)*N1)^2);
S(1,2)=trace(inv(N)*N1*inv(N)*N2);
S(2,1)=trace(inv(N)*N1*inv(N)*N2);
S(2,2)=row2-2*trace(inv(N)*N2)+trace((inv(N)*N2)^2);
XI=inv(S)*W;
%以上代码利用验后方差估计定权，再进行平差，不再详解代码
end
tX=X;
[n,t]=size(B);
tmo=sqrt((V'*P*V)/(n-t)); %获取单位权中误差
tQ=inv(B'*P*B); %未知数协因数矩阵
tM=tmo^2*tQ; %未知数方差阵
[row,col]=size(XY);
for j=1:row %在XY 矩阵中进行循环
wzindex=XY(j,5); %获取第j 个要素对应未知数序号 if wzindex~=0 %如果是未知点
XY(j,15)=tXf(2*wzindex-1,1)/1000+XY(j,15);
XY(j,16)=tXf(2*wzindex,1)/1000+XY(j,16);
XY(j,17)=tMf(2*wzindex-1,2*wzindex-1);
XY(j,18)=tMf(2*wzindex,2*wzindex); %无定权平差坐标及精度 XY(j,19)=tX(2*wzindex-1,1)/1000+XY(j,11);
XY(j,20)=tX(2*wzindex,1)/1000+XY(j,12);
XY(j,21)=tM(2*wzindex-1,2*wzindex-1);
XY(j,22)=tM(2*wzindex,2*wzindex); %有估权平差坐标及精度
end
end
end
6.3 基于线性模型的CPIII 控制网平差
6.3.1 CPIII 控制网平差线性模型推理
图6-5:基线变换示意图
如图6-5所示，设在测站点j 对待测点k 进行观测，利用其观测方位和观测距离可获取基线向量k A j A ，设平差后该基线向量为k B j A 。

则
k k k B B A k j j A A A -= (6-7)
在k B j A 上取一点0A ，使得1A 0j =A ，过0A 做k B k A 平行线交k B j A 于0B 点，则
θ B o
A o
B k A k A j
k k
B B A A k 00j 0
j A A A A = (6-8) 由(6-8)式得 00j k A A A B A B k k ⨯= (6-9)
则
00j k A A A B A B k k ⨯= (6-10) 又由于0j 00A R A A B ⨯=，其中
⎥⎦
⎤⎢⎣⎡-=⎥⎦⎤⎢⎣⎡-⨯=a b b a θθ
θθcos sin sin cos k R (6-11) 将(6-11)式代入(6-10)得 k k k A R A R A B j 0j j k A A A A ⨯=⨯⨯= (6-12)
将(6-12)代入(6-7)式得
k k k A R B A j j j A A A ⨯-= (6-13)
则
⎥⎦⎤⎢⎣⎡∆∆⨯⎥⎦⎤⎢⎣
⎡--⎥⎦⎤⎢⎣⎡-⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡∆∆jk jk jk jk Y X Y X a b b a V U V U j j k k (6-14) ⎥⎦
⎤⎢⎣⎡⨯⎥⎦⎤⎢⎣⎡∆∆∆-∆-⎥⎦⎤⎢⎣⎡-⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡∆∆b a V U V U j j k k jk jk jk jk jk jk X Y Y X Y X (6-15) (6-15)式即为CPIII 控制网平差的线性模型，其中左侧⎥⎦
⎤⎢⎣⎡∆∆jk jk Y X 为随机量，而右侧⎥⎦
⎤⎢⎣⎡∆∆jk jk Y X 为非随机量。

(6-15)对应误差方程为
⎥⎦
⎤⎢⎣⎡∆∆-⎥⎦⎤⎢⎣⎡⨯⎥⎦⎤⎢⎣⎡∆∆∆-∆-⎥⎦⎤⎢⎣⎡-⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡jk jk jk jk jk jk jk jk Y X X Y Y X b a V U V U V V j j k k (6-16) 设总观测点数（包括重复点）为n ，则有2n 个误差方程，写成矩阵形式为
l BX -=V
其中
[]T
Yn Xn Y1X1V V V V V = ⎥⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎣⎡------= n n n n
ΔX ΔY ΔY ΔX ΔX ΔY ΔY ΔX B 10100101101001011111 []T
11j1j1k1k1b a V U V U X =
[]T n n 11ΔY ΔX ΔY ΔX l = ⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=∆∆∆∆∆∆∆∆∆∆∆∆∆∆∆∆ Yn Yn Yn Xn Yn Xn Xn Xn Y1Y1Y1
X1Y1X1X1X1P P 000
P P 000
000P P 000P P B 设待定点个数为m ，则未知参数个数为2m 。

在观测值总数大于或等于未知参数的情况下，
方程有解，即2n> 2m ，得n>m 。

利用间接平差原理可得 Pl B PB B T
T 1)(X -=。

协因数阵1)(Q -=PB B T 。

单位权方差因子为f
PV V n T )(0=。

转换后各待求参数精度矩阵为Q n M 0=。

从以上推理可知，基于线性模型的CPIII 控制网平差模型不需要计算近似值，不需要迭代计算，因此计算量将大为减小，同时线性模型在形式上也更加简单。

6.3.2 基线向量观测值权的确定
线性模型采用了基线向量作为虚拟观测值，因此本节首先推导出基线向量和原始观测值的微分关系，再利用最小范数二次无偏估计确定基线向量观测值的权比关系。

设测站点为j ，待测点为k 。

则基线向量 ⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡∆∆jk jk
jk jk jk jk sin cos Y X ααD D (6-17)
对(6-17)式取微分得
jk jk jk jk jk jk jk jk jk jk sin cos 1000cos sin Y X dD d D D d d ⎥⎦⎤⎢⎣
⎡+⎥⎦⎤⎢⎣⎡-=⎥⎦⎤⎢⎣⎡∆∆ααρααα (6-18) 设所测点(包括重复点)个数为n ，则可写出n 个(6-18)式，合写为矩阵形式为 dD H dA H dL 21+= (6-19) (6-19)式中
[]T X X Y X L n n 11d d d d d ∆∆∆∆= []T n d d αα 1dA = []T n dD dD D 1d = ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--=ραραραραn n n n D D D D H cos 100000sin 1000000000cos 100000sin 100011111 ⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=n n H ααααsin 00cos 000000sin 00cos 112 则由文献可得
T H H T 111⨯=，T H H T 222⨯=，21T T T += 11111)(------=T B B T B B T T C T T ⎥⎦⎤⎢⎣⎡=)()()()(22122111CT CT tr CT CT tr CT CT tr CT CT tr S ⎥⎦⎤⎢⎣⎡=CV CT V CV CT V W T T 21 W S 1-=θ
其中
⎥⎦
⎤⎢⎣⎡=22D δδθα
2α
δ和2D δ分别为方位观测值和距离观测值的单位权方差因子，则其中误差矩阵分别为1M 和2M 。

E 21αδ=M ，E 2d 2δ=M
则两类观测值权分别为
1121M -=αδP ，12
22M -=d P δ 由此可得基线向量的协因数阵为
T T H P H H P H Q 21221111--+=
可得基线观测向量的权阵-1Q =P 。

由于最小范数二次无偏估计不需要进行迭代计算，因此计算工作量比较小。

6.3.3 模拟数据验证
为验证本节所提出线性模型的可行性，利用程序产生模拟数据，并对模拟数据进行计算处理。

通过分析各点的点位中误差以及每对相邻点相对中误差的情况，来判断数据处理的可靠性。

首先利用计算机产生CPIII 点，按规范每隔50～60 m 左右布设一对，每对间距约为15m 。

且每500 m 左右设置一对CPII 点，作为CPIII 控制网的约束点[1]。

为模拟真实的情况，CPIII 控制网的位置和方向都是随机确定的。

然后利用CPIII 点坐标值来产生模拟观测数据。

第一站在开始的4个CPIII 点间设站，测站位置由4个CPIII 点中心点坐标加上一定随机数产生，然后利用各CPIII 点坐标和测站点坐标反算方位和距离观测值。

对该测站所有反算方位角加上同一定向角，定向角随机产生。

假定观测仪器标称精度已知，测角精度为1”，测距为1mm+1ppm 。

根据仪器精度对产生的模拟观测数据附加随机误差值。

第二测站观测在开始的8个CPIII 点间进行，第三测站在开始的12个点间进行，模拟观测数据产生同前。

以后每测站要求观测12个CPIII 点，且保证对所有CPIII 点观测3次。

另外测段末尾3站要按初3站的逆顺序进行，即倒数第三站观测12个CPIII 点，倒数第二站观测CPIII8个点，倒数第一站观测4个CPIII 点。

本次模拟实验产生CPIII 点60个，CPII 点8个，每测站观测3个测回。

图6-6为各CPIII 点点位中误差曲线图，图6-7为各相邻点间中误差曲线图。

图6-6 平差后各点点位中误差示意图
图6-7 平差后各相邻点相对中误差示意图
从图6-7可知，平差后各点点位中误差大部分都在1mm 之内，利用最小范数二次无偏估计对基线向量观测值重新定权后，各点点位精度略有提高。

而从图6-7可知，经过二次无偏估计重新定权后，相邻点相对中误差更加集中，且都小于1mm ，满足规范要求。

通过大量模拟实验可知，本节提出的线性模型对于CPIII 控制网的平差是可行的。

本节在分析传统CPIII 控制网平差方法的基础上，提出一个线性模型。

由于不需要计算近似值，不需要迭代计算，因此计算工作量相对于传统处理方法有较大优势。

同时该模型在形式上比传统模型更简单。

为了验证该线性模型的正确性，笔者编写了程序代码，利用模拟数据进行了大量验证，由实验结果证明该方法的可行性。

由于受篇幅的限制，本节只列出了其中一次模拟实验的部分数据。

点号点位中误差
相邻点对序号相邻点相对中误
差
6.3.4 代码详解
程序主界面如图6-8所示，包括4个消息激发按钮，分别为产生模拟CPIII点、产生模拟观测数据、无约束平差、约束平差。

图6-8 CPIII平差程序主界面
%函数功能：主程序Main，以下包括四个按钮的消息函数
function pushbutton3_Callback(hObject, eventdata, handles)
ptnum=60; %确定CPIII控制点个数，可以修改jjsq=0; %奇数点计数器置零
ojsq=0; %偶数点计数器置零
for i=1:ptnum %进行循环
if fix(i/2)==i/2 %当i为偶数时
XY(i,1)=i; %点号保存到第一列
XY(i,2)=ojsq*60+normrnd(0,1) ; %获取纵坐标，并加随机数
XY(i,3)=30+ normrnd(0,1); %获取横坐标，并加随机数
ojsq=ojsq+1; %偶数计数器加1
else %当i为奇数
XY(i,1)=i; %保存点号
XY(i,2)=jjsq*60+ normrnd(0,1); %第二列保存纵坐标
XY(i,3)=15+ normrnd(0,1); %第三列保存横坐标
jjsq=jjsq+1; %奇数计数器加1
end
end
%以上代码产生CPIII点阵列坐标数据
angle=rand(1)*2*pi(); %随机产生一角度angle
xzJZ(1,1)=cos(angle);
xzJZ(1,2)=-sin(angle);
xzJZ(2,1)=sin(angle);
xzJZ(2,2)=cos(angle); %生成二维旋转矩阵
for j=1:ptnum %在点要素表中循环
temXY(j,1)=XY(j,2);
temXY(j,2)=XY(j,3); %将CPIII点坐标存入临时变量中end
getXY=(xzJZ*temXY')'; %将CPIII点坐标进行坐标旋转XY=0; %将XY矩阵置零
for i=1:ptnum %对旋转后点坐标进行循环处理 XY(i,1)=i; %获取点号
XY(i,2)=getXY(i,1); %第二列保存纵坐标
XY(i,3)=getXY(i,2); %第三列保存横坐标
end
%以上代码将XY矩阵中坐标值按随机角度进行旋转
[XY]=addCP(XY,ptnum,1,1);
[XY]=addCP(XY,ptnum,2,2);
[XY]=addCP(XY,ptnum,21,3);
[XY]=addCP(XY,ptnum,22,4);
[XY]=addCP(XY,ptnum,38,5);
[XY]=addCP(XY,ptnum,39,6);
[XY]=addCP(XY,ptnum,58,7);
[XY]=addCP(XY,ptnum,59,8);
%以上代码在CPIII点阵列中增加若干控制点
jsq=8; %控制点个数为8
jsq1=0; %计数器置零
for i=1:ptnum+jsq %在CPIII点表中进行循环
if XY(i,4)==0 %如果是未知点
jsq1=jsq1+1; %计数器加1
XY(i,5)=jsq1; %这个是未知点序号，放第5列 end
XY(i,6)=i; %第6列存放所有点序号
XY(i,16)=0; %第16列备用，此句可选
end
save ('datafile','XY'); %保存模拟点坐标
msgbox '成功产生模拟观测数据！'; %提示程序执行成功
%以上代码是产生CPIII模拟点，为后续产生模拟观测数据做好基础工作
%--------------------------------------------------
function pushbutton1_Callback(hObject, eventdata, handles)
clear; %清空已有变量
format long e; %定义变量类型
rou=3600*180/pi(); %弧度化秒参数
load('datafile','XY'); %加载模拟点坐标文件
[gcJZ]=creatsurvey(XY); %根据模拟点产生模拟观测数据
save ('gcfile','gcJZ'); %将模拟观测数据保存
clear; %清空变量
msgbox '成功产生模拟点！'; %提示执行成功
%以上代码根据模拟CPIII点产生模拟观测数据
%----------------------------------------------------
function pushbutton2_Callback(hObject, eventdata, handles)
format long e; %设置变量类型
load ('gcfile','gcJZ'); %加载观测数据
load('datafile','XY'); %加载CPIII点表
[base,czs]=getBase(gcJZ); %构建基线向量
[B,l]=getV(czs,base,XY,0); %解算误差方程系数矩阵
l=l*1000; %改正数和未知数参数单位取毫米
X=inv(B'*B)*B'*l; %解算未知数
[X,M,mo,XY]=lastpc(gcJZ,XY,B,l,base,czs);%进行无约束平差
draw(X,czs); %绘制误差分布曲线
clear; %清除变量
msgbox '无约束平差完成！'; %提示解算完毕
%以上代码实现CPIII控制网的无约束平差
%--------------------------------------------------
function pushbutton4_Callback(hObject, eventdata, handles)
clear; %清空现有变量
format long e; %设置变量类型
load ('gcfile','gcJZ'); %加载观测数据
load('datafile','XY'); %加载CPIII点表
[base,czs]=getBase(gcJZ); %获取基线向量
[B,l]=getV(czs,base,XY,1); %解算误差方程系数矩阵
l=l*1000; %改正数和未知参数坐标单位取mm
ceZB=0; %测站坐标矩阵
[X,Q,mo,Qf,mof,XY,V,ceZB,xzJD]=lastpc3(gcJZ,XY,B,l,base,czs); %进行约束平差[X,tQ,tmo,tQf,tmof,XY,V]=nonelinarpc1(gcJZ,XY,czs,ceZB,xzJD,Q); %进行约束平差[Mbp]=getbetweenM(XY,mo,Q); %相邻点方差矩阵（定权后）
[Mbpf]=getbetweenM(XY,mof,Qf); %相邻点方差矩阵（无定权）
[r,c]=size(Mbp); %获取方差矩阵行数
jsq=0; %计数器置零
x=0; %绘图横坐标矩阵
y1=0; %绘图纵坐标矩阵（估权）
y2=0; %绘图纵坐标矩阵（无估权）
for i=1:r/2 %在方差矩阵中循环
jsq=jsq+1; %计数器加1
mx=Mbp(2*i-1,3); %无估权纵坐标方差
my=Mbp(2*i,3); %无估权横坐标方差
mxf=Mbpf(2*i-1,3); %估权纵坐标方差
myf=Mbpf(2*i,3); %估权横坐标方差
x(jsq,1)=jsq; %绘图横坐标矩阵
y1(jsq,1)=sqrt(mx^2+my^2); %绘图纵坐标估权点位中误差
y2(jsq,1)=sqrt(mxf^2+myf^2); %绘图纵坐标无估权点位中误差
end
xaxis='相邻点对序号'; %绘图横坐标标签
yaxis='相邻点相对中误差'; %绘图纵坐标标签
figh1=0;
[figh1]=drawcurve(x,y1,y2,xaxis,yaxis,figh1); %调用绘图函数
%以上代码绘制相邻点位中误差分布图
[r,c]=size(XY); %获取XY矩阵行数
jsq=0; %计数器置零
x=0; %绘图横坐标变量置零
y1=0; %绘图纵坐标变量置零（估权）
y2=0; %绘图纵坐标变量置零（无估权）
for i=1:r-8 %循环处理非已知点
jsq=jsq+1; %计数器加1
mx=XY(i,13); %纵坐标中误差（估权）
my=XY(i,14); %横坐标中误差（估权）
mxf=XY(i,9); %纵坐标中误差（无估权）
myf=XY(i,10); %横坐标中误差（无估权）
x(jsq,1)=jsq; %绘图横坐标矩阵变量
y1(jsq,1)=sqrt(mx^2+my^2); %绘图纵坐标矩阵变量（估权）
y2(jsq,1)=sqrt(mxf^2+myf^2); %绘图纵坐标矩阵变量（无估权）
end
xaxis='点号'; %绘图横坐标标签
yaxis='点位中误差'; %绘图纵坐标标签
figh2=0;
[figh2]=drawcurve(x,y1,y2,xaxis,yaxis,figh2); %调用绘图函数
msgbox '约束平差完成！';
%-------------------------------------------------
%函数功能：添加控制点子程序
%输入变量：jsq为控制点序号，pn为插入控制点之参照点，即在哪个CPIII点附近插入一%控制点，ptnum为CPIII点个数，XY为CPIII点表矩阵。

输出变量亦为XY
function [XY]=addCP(XY,ptnum,pn,jsq)
x0=XY(pn,2);
y0=XY(pn,3); %获取参照点点名
x1=XY(pn+1,2);
y1=XY(pn+1,3); %获取参照点下一个点，作为定向点
detx=x1-x0;
dety=y1-y0; %获取定向点和参照点坐标增量
fw=getfw(detx,dety); %获取定向点相对于参照点的方位角
dist=sqrt(detx^2+dety^2)+50*randn(1); %获取定向点相对于参照点的距离并加随机数XY(ptnum+jsq,1)=1000+jsq; %确定控制点点名
XY(ptnum+jsq,2)=x1+dist*cos(fw);
XY(ptnum+jsq,3)=y1+dist*sin(fw); %获取控制点纵横坐标值
XY(ptnum+jsq,4)=1; %未知点标志列值设为1
end
%-------------------------------------------------
%函数功能：联测控制点，决定在哪些测站需要引入控制点
%输入变量：CPIII点表，控制点名，当前测站点观测之CPIII点矩阵，输出变量相同
function [mytem,name,XY]=connectContr(mytem,name,XY)
[row0,col0]=size(mytem); %获取当前测站观测点矩阵行列数
[r,c]=size(XY); %获取CPIII点表行列数
sumx=0;
sumy=0;
jsq=0; %计数器置零
for v=1:row0 %在mytem矩阵中循环
sumx=sumx+mytem(v,1);
sumy=sumy+mytem(v,2); %将该测站所以CPIII点纵横坐标求和
end
stnx=sumx/row0;
stny=sumy/row0; %此即为该测站CPIII点阵之中心坐标
for i=1:r %在XY矩阵中循环
if XY(i,4)==1 && XY(i,16)~=-1 %如果第i个点是已知点
pc=XY(i,1);
xc=XY(i,2);
yc=XY(i,3); %获取第i个点的点名已经坐标值
dist=sqrt((xc-stnx)^2+(yc-stny)^2);
%获取第i个点相对于测站距离
if dist<150 %如果距离小于150米
jsq=jsq+1; %计数器加1
XY(i,16)=-1; %第16列设为-1
mytem(row0+jsq,1)=xc;
mytem(row0+jsq,2)=yc;
name(row0+jsq,1)=pc; %将第i个点保存到mytem矩阵中
end
end
end
end
%以上代码实质上是将控制点进行划分，根据具体距离划分到最近测站，以便后面统一计算%观测数据，其中XY矩阵第16列存放一个标志变量，默认值为0，如果是-1，表明这个点已%经被联测过了，在别的测站不需要考虑这个点。

%-------------------------------------------------
%函数功能：产生模拟观测数据子程序
%输入变量为CPIII点坐标矩阵，输出变量为观测值矩阵gcJZ
%本函数循环从模拟的CPIII点阵中提取若干点作为本测站观测目标，然后根据模拟坐标值%产生一测站点，再根据CPIII点和测站点模拟坐标反算观测数据，直至所有CPIII点观测%次数达3次
function [gcJZ]=creatsurvey(XY)
format long e; %设置变量类型
rou=3600*180/pi(); %弧度化秒常数
[row,col]=size(XY); %获取XY矩阵行列数
xydata=0; %xydata矩阵置零
for i=1:row %在XY矩阵中循环。