利用便携式计算机编程快速计算道路中线坐标

利用便携式计算机编程快速计算道路中线坐标
摘要：本文介绍了道路中线坐标计算公式，实现整条线路统一连续计算中桩坐标的编程思路。

采用本文的思路编写程序，能快速求解中线上任一点坐标，给道路放样工作带来了便捷。

关键词：中桩坐标；计算公式；编程；统一计算；放样
Abstract: This paper introduces the calculation formula of coordinate center of road, realize the unity of the whole line continuous programming calculation of pile coordinates. By using this program, can quickly solve the midline of any point coordinates, to the road layout work brings convenient.
Key words: pile coordinate; calculation formula; programming; unified computing; lofting
1前言
近几年来，随着我国经济社会高速发展，道路交通建设步伐不断加快，各级别特别是高等级公路建成了很多，更多的道路工程在陆续上马。

在道路建设初期和施工期，都需要我们测绘人员对其进行中线放样。

中线一般由直线段、缓和曲线段、圆曲线段组成，设计文件中提供线路起点、交点、终点坐标以及曲线半径、缓和曲线长度、平曲线要素、曲线主点桩号等。

利用传统方法进行放样，特别是曲线放样，步骤烦，速度慢，放样精度不高，还要受通视条件影响，野外作业时，如果需要用坐标来精确放样中桩，得现场分步计算，输入数据多，费时且容易算错。

在生产实践中，笔者通过编程实现了整条线路统一连续计算中桩坐标，在道路中线放样时，只需以线路中线的桩号为引数就可求得中线上相应点的坐标，采用全站仪、GPS RTK等先进的测绘仪器进行线路中线放样极为便利，大大提高了中线放样的效率和准确度。

2计算公式
2.1 直线中线坐标计算
已知直线起点坐标（X0，Y0）和直线前进方位角F，则直线前进方向上距起点S的待求点坐标（Xi，Yi）为：
Xi=X0＋S·CosF
Yi=Y0＋S·SinF (1)
2.2 曲线中线坐标计算
在推导曲线中线坐标计算公式之前，先将曲线要素计算公式叙述如下：
无缓和曲线的圆曲线的曲线要素计算公式为：
T=R·Tanα2
L= R·α (2)
E=R·Secα2－R
q=2·T－L
（2）式中：R为曲线半径，α为线路转角（以弧度为单位，角度单位下同），T为切线长，L为曲线长，q为切曲差。

有缓和曲线的圆曲线的曲线要素计算公式为：
T=m＋(R＋p)·Tanα2β0= L0 /(2·R)
L= R·（α－2·β0）＋2·L0 (3) m= L0/2－L03/(240·R2)(4)
E=（R＋p）·S ecα2－Rp= L02/(24·R)
q=2·T－L
(3)式中L0为缓和曲线长（两边相等）；m为加设缓和曲线后使切线增长的距离；p为加设缓和曲线后，圆曲线相对于切线的内移量；β0为缓和曲线角度。

m、p、β0为缓和曲线参数，按(4)式计算。

2.2.1 无缓和曲线的圆曲线中线坐标计算
如图1所示，圆曲线上任一点P的坐标（Xi，Yi），可利用切点ZY的坐标、切线方位角（即ZY至JD方位角F）、弦线与切线夹角k以及弦长S来计算，而切点ZY的坐标可由已知的交点坐标、T长和切线方位角F来推算。

设Li为点P 至切点的弧长，则：
k=Li/(2·R)
S=2·R·Sink
Xi=XZY＋S·Cos(F＋w·k) (5)
Yi=YZY＋S·Sin(F＋w·k)
式中：w为线路偏转方向系数，线路左转w=－1，线路右转w=＋1，下同。

图1图2
2.2.2 缓和曲线中线坐标计算
缓和曲线可分为两种情况。

一是用来连接直线和圆曲线的缓和曲线，其曲率半径从无穷大到圆曲线半径；二是用来连接不同半径圆曲线之间的缓和曲线。

在一般的公路线形中很少碰到后者（笔者暂时还未遇到这种复合曲线），在此，我们只讨论前者中线坐标计算公式。

如图2所示，以直缓点ZH为原点，过ZH的缓和曲线切线为x轴，ZH点上缓和曲线半径为y轴，建立相对直角坐标系。

设Li为缓和曲线上任一点P至ZH点的曲线长，则点P在此坐标系中的坐标为：
xi=Li－Li5/(40·R2·L02)＋…
yi=－w·[Li3/(6·R·L0)－Li7/(336·R3·L03)＋…](6)
对其进行坐标转换，得到该点在国家统一的或地方独立的平面坐标系中的坐标：
Xi=XZH＋yi·cosθ－xi·Sinθ
Yi= YZH＋yi·Sinθ＋xi·Cosθ (7)
(7)式中：θ为相对直角坐标系逆时针旋转角度，等于ZH至JD方位角F－π/2。

（XZH，YZH）为ZH点的坐标，可由已知的交点坐标、T长和F来推算。

2.2.3有缓和曲线的圆曲线中线坐标计算
对于有缓和曲线的圆曲线中线坐标计算，关键是先求出圆曲线起点坐标及该点处的切线方位角，然后就可按（5）式来计算。

此时，Li为圆曲线上任一点至
圆曲线起点的弧长。

如图2所示，圆曲线起点HY坐标可由ZH点坐标和缓和曲线长L0按（6）式（7）式计算，切线方位角等于ZH点至JD的方位角F＋w·β0。

3编程思路
为了使整条线路能统一起来连续计算中桩坐标，首先得定义此条线路。

在程序中把交点的坐标、曲线半径、缓和曲线长度定义为一维数组，按设计文件输入这些已知数据，计算出交点间距、线路前进方位角、转角、曲线长、切线长以及交点里程等线路其它要素，也保存在一维数组中。

线路转角有正负性，为正角时线路右转，为负角时线路左转。

这样，整条线路的位置、线形、走向就在计算机中唯一地存在了。

接下来，输入线路上任一中桩的桩号，判断出该中桩位于整条线路中哪段中线上，调用相应计算模块（计算本文2中几种线形中桩坐标的过程或函数），求出其坐标。

如何判断输入的中桩桩号位于整条线路中哪段中线上呢？这是整个程序的关键。

根据里程，利用循环语句很容易判断出待求中桩位于哪两个交点间，然后判断是位于小号端的曲线上，还是位于大号端的曲线上，还是位于中间的直线上。

如在曲线上，判断其位于缓和曲线上，还是位于圆曲线上。

如图3所示，L(i)、L(i+1)分别是JD(i)、JD(i+1)的里程，LYH 、LHZ、LZH、LHY分别为YH点、HZ点、ZH点、HY点的里程，设待求中桩里程为LC，则上述判断过程可用图4表示。

在程序中加入放样计算模块，求出中桩坐标后，就可解算出放样数据——方位角和距离。

图3
图4：流程图
4结束语
1995年7月，笔者利用上述的方法和思路，在SHARP PC—E500袖珍计算
机固化BASIC环境中编写出程序，实现了整条线路中桩坐标统一连续计算，可计算线路上任意里程的中桩坐标和放样数据。

首次应用于苏州观山机动车驾驶员培训基地道路测量项目中，在生产实践中检验了程序。

此后，该程序广泛应用于各个道路测量项目中，为线路中线放样提供了便利，提高了外业效率。

当然，可进一步扩充程序的功能，如复合曲线中线坐标计算、线路边桩坐标计算等，使程序更为完善。

也可用其它高级语言编程，在其它类型掌上电脑或笔记本电脑上运行计算。

参考文献
[1]李青岳主编.工程测量学.北京：测绘出版社，1984年6月
说明：笔者1993年8月至2001年12月，在江苏苏州地质工程勘察院任职，从事工程测量工作。

附PC—E500 BASIC语言环境下程序代码：
8 CLEAR
10 A$=“#######.###”:B$=“########.###”:C$=“####.###”:D$=“###.####”
15 INPUT “Jiao Dian Ge Shu: “;N
20 DIM X(N+1),Y(N+1),L(N+1),D(N),Z(N),F(N),HC(N),R(N+1),L0(N),T(N)
25 WAIT 0
30 PRINT “QI DIAN ZHUO BIAO”:INPUT “X=“;X(0):INPUT “Y=“;Y(0)
40 FOR I=1TO N
45 PRINT “JIAO DIAN”;I;”ZHUO BIAO”
50 PRINT “X”;I;”=“;:INPUT X(I)
60 PRINT “Y”;I;”=“;:INPUT Y(I)
70 NEXT I
80 PRINT “ZHONG DIAN ZHUO BIAO”:INPUT “X=“;X(N+1):INPUT “Y=“;Y(N+1)
90 FOR I=0TO N
100 X=X(I+1)-X(I):Y=Y(I+1)-Y(I)
110 D(I)=POL (X,Y)
120 IF Z180 LET F(I)=ZI-360:GOTO 150
146 IF ZI0THEN PRINT “ L0=“;L0(I)
270 X=T(I):GOSUB 1000:PRINT “T=“;X
271 X=HC(I):GOSUB 1000:PRINT “L=“;X
272 X=E:GOSUB 1000:PRINT “E=“;X
273 X=Q:GOSUB 1000:PRINT “q=“;X
275
ZH=L(I)-T(I):HY=ZH+L0(I):QZ=ZH+HC(I)/2:HZ=ZH+HC(I):YH=HZ-L0(I) 278 IF ZH=HY THEN 290
280 PRINT “ ZH”:A=ZH:GOSUB 900
282 PRINT “ HY”:A=HY:GOSUB 900
283 PRINT “ QZ”:A=QZ:GOSUB 900
284 PRINT “ YH”:A=YH:GOSUB 900
285 PRINT “ HZ”:A=HZ:GOSUB 900:GOTO 300
290 PRINT “ ZY”:A=ZH:GOSUB 900
291 PRINT “ QZ”:A=QZ:GOSUB 900
295 PRINT “ YZ”:A=HZ:GOSUB 900
300 NEXT I
310 X=D(N):Y=DMS Z(N):GOSUB 1000
320 CLS :WAIT 0:PRINT “S(“;N;”-”;N+1;”)=“;USING C$;X;:PRINT “a(“;N;”-”;N+1;”)=“;USING D$;Y:WAIT :PRINT
330 WAIT 0:PRINT “ZHONG DIAN”
340 A=L(N+1):GOSUB 900
350 PRINT “X=“;USING A$;X(N+1);:PRINT “ Y=“;USING B$;Y(N+1):WAIT :PRINT
360 CLS :WAIT 0:PRINT “FANG YANG JI SUAN”
362 PRINT “JIA ZHAN DIAN”
363 PRINT “ZHONG XIAN SHANG?(Y OR N)”;:INPUT AA$
364 IF AA$=“Y”THEN GOSUB 385: XA=XB:YA=YB ELSE INPUT “XA=“;XA,”YA=“;YA
366 PRINT “DING XIANG DIAN”
367 PRINT “ZHONG XIAN SHANG ?(Y OR N)”;:INPUT AA$
368 IF AA$=“Y”THEN GOSUB 385ELSE INPUT “XB=“;XB,”YB=“;YB
369 GOSUB 1300:WAIT :PRINT
370 CLS :WAIT 0
372 GOSUB 385:GOSUB 1300
375 PRINT “HUAN ZHAN?(Y OR N)”;:INPUT LL$
380 IF LL$=“Y” THEN 360 ELSE 370
385 D=1000
390 PRINT “ZHONG XIAN ZHUANG HAO:”;:INPUT LC
400 IF L(0)-0.001HZ LET S=LC-HZ+T(I):XB=REC (S,Z(I))+X(I):YB=Z+Y(I):GOTO 710
600 IF LC=INT (HZ*D+0.5)/D LET XB=XN:YB=YN:GOTO 710
610 B=XN:C=YN:B0=O+90:O=O+180
620 R=R(I):L0=L0(I)
630 IF F(I)<0 LET W=1 ELSE W=-1
650 IF L0=0 LET L=HZ-LC:GOSUB 1500:GOTO 710
660 YH=HZ-L0:S=L0:GOSUB 1700:XYH=XN+TX:YYH=YN+TY
670 IF LC=INT (YH*D+0.5)/D LET XB=XYH:YB=YYH:GOTO 710
680 IF YH<LC AND LC<HZ LET S=HZ-LC:GOSUB 1700:XB=XN+TX:YB=YN+TY:GOTO 710
690 L=YH-LC:O=O+90*L0*W/PI /R:B=XYH:C=YYH
700 GOSUB 1500
710 PRINT “X=“;INT (XB*D+0.5)/D;”Y=“;INT (YB*D+0.5)/D
720 RETURN
840 M=L0(I)/2-L0(I) /(240*R(I) ):P=L0(I) /(24*R(I)):B0=90*L0(I)/(PI *R(I))
850 T(I)=(R(I)+P)*TAN (F/2)+M
860 HC(I)=PI *R(I)*(F-2*B0)/180+2*L0(I)
870 E=(R(I)+P)/COS (F/2)-R(I)
880 Q=2*T(I)-HC(I)
890 RETURN
900 A1=INT (A/1000):A1$=STR$ (A1):A2=A-A1*1000
910 PRINT “K”;A1$;”+”;USING C$;INT (A2*1000+0.5)/1000
920 RETURN
1000 X=INT (X*1000+0.5)/1000
1010 Y=INT (Y*10000+0.5)/10000
1020 RETURN
1300 G=XB-XA:H=YB-YA:S=POL (G,H):X=S
1310 IF Z<0 LET Z=Z+360
1320 Z=DMS Z:Y=Z
1330 GOSUB 1000
1340 PRINT “S(a-b)=“;X;”a(a-b)=“;Y 1700 DX=S-S /40/R/R/L0/L0 1350 RETURN 1710 DY=S /6/R/L0-S /336/R /L0
1500 K=L*90/PI /R1720 DY=DY*(-W)
1510 J=2*R*SIN K1730 TX=DY*COS B0-DX*SIN B0
1520 XB=REC (J,O+K*W)+B 1740 TY=DY*SIN B0+DX*COS B0 1530 YB=Z+C:RETURN1750 RETURN。