移动坐标

2024-10-24

移动坐标（精选4篇）

移动坐标篇1

0 引言

一般曲线坐标系下的二维水沙模型,当涉及河床变形及岸坡崩退时,崩岸的宽度与计算网格的尺寸存在一定的差距。如要精确考虑崩岸尺寸,崩岸处网格必须尽量小,但对于整个河道而言,会导致计算速度过慢。为解决实际存在的问题,国内外一些学者进行了移动坐标下可变网格的研究,并取得了较为理想的成果。日本Nagata[1]等学者首先提出了移动坐标下模拟河岸侵蚀过程的数值方法。Jang和Shimizu[2]采用移动坐标下的平面二维水沙数学模型对宽浅型的水流形态,河床变形及河岸的侧向冲蚀变形进行了数值模拟研究。杨建明[3]采用可变网格技术模拟挑流冲刷时因河床变形引起的流场变化。

河流弯道水流呈现明显的三维特性,平面二维水沙模型无法很好地模拟由于弯道横向输沙引起的凹岸冲刷、凸岸淤积的现象。三维数学模型的研究虽然取得了一些成果,但许多技术性问题尚待解决[4];另外一种途径就是对现有的平面二维数学模型进行修正[5],现有的研究成果在文献[5]中作了比较详细的阐述。本文在文献[7]建立的曲线平面二维数学模型基础之上,引入文献[8]中移动坐标,并采取文献[5]的研究方法,考虑悬移质的横向输沙,对河岸岸坡崩塌的模拟,引入夏军强[6]改进后的黏性土河岸变形计算方法。通过上述这些方法对已有平面二维数学模型进行修正。

1 数学模型

1.1 基本方程

移动坐标下,方程中的沿时间偏导项为:

$\frac{\partial}{\partial t} = \frac{\partial}{\partial τ} + \frac{\partial}{\partial ξ} \frac{\partial ξ}{\partial τ} + \frac{\partial}{\partial η} \frac{\partial η}{\partial τ} (1)$

因此在移动一般曲线坐标下,平面二维水沙数学模型的基本方程如下[8]:

连续方程:

$J \frac{\partial h}{\partial τ} + \frac{\partial (U + ξ_{t}) h}{\partial ξ} + \frac{\partial (V + η_{t}) h}{\partial η} = J q (2)$

ξ方向运动方程:

$\begin{array}{l} J \frac{\partial h U}{\partial τ} + \frac{\partial (U + ξ_{t}) Μ}{\partial ξ} + \frac{\partial (V + η_{1}) Μ}{\partial η} = \\ - g h (\frac{\partial Ζ}{\partial ξ} ξ_{x} + \frac{\partial Ζ}{\partial η} η_{x}) J + \frac{\partial}{\partial ξ} [D J (q_{11} Μ_{ξ} + q_{12} Μ_{η})] + \\ \frac{\partial}{\partial η} [D J (q_{12} Μ_{ξ} + q_{22} Μ_{η})] - \frac{g n^{2} Μ J \sqrt{u^{2} + v^{2}}}{h^{4 / 3}} + J q u_{0} (3) \end{array}$

η方向运动方程:

$\begin{array}{l} J \frac{\partial h V}{\partial τ} + \frac{\partial (U + ξ_{t}) Ν}{\partial ξ} + \frac{\partial (V + η_{1}) Μ}{\partial η} = \\ - g h (\frac{\partial Ζ}{\partial ξ} ξ_{x} + \frac{\partial Ζ}{\partial η} η_{y}) J + \frac{\partial}{\partial ξ} [D J (q_{11} Ν_{ξ} + q_{12} Ν_{η})] + \\ \frac{\partial}{\partial η} [D J (q_{12} Ν_{ξ} + q_{22} Ν_{η})] - \frac{g n^{2} Ν J \sqrt{u^{2} + v^{2}}}{h^{4 / 3}} + J q v_{0} (4) \end{array}$

泥沙连续性方程:

$\begin{array}{l} J \frac{\partial h S_{k}}{\partial τ} + \frac{\partial (U + ξ_{t}) h S_{Κ}}{\partial ξ} + \frac{\partial (V + η_{t}) h S_{Κ}}{\partial η} = \\ \frac{\partial}{\partial ξ} [D J (q_{11} h S_{Κ ξ} + q_{12} h S_{Κ η})] + \\ \frac{\partial}{\partial η} [D J (q_{12} h S_{Κ ξ} + q_{22} h S_{Κ η})] - J α W_{Κ} (S_{Κ} - S_{* Κ}) + J q s_{0} (5) \end{array}$

河床变形方程:

$γ^{'} \frac{\partial Ζ_{b}}{\partial τ} = \sum α ω_{Κ} (S_{Κ} - S_{* Κ}) (6)$

式中:x、y为物理平面坐标;ξ、η为贴体坐标下的坐标;D为紊动黏性系数;n为Manning糙率系数;u、v分别为x 和y方向的流速,m/s;J为雅可比变换系数,即J=xξyη-xηyξ;U和V分别为曲线坐标中流速在ξ、η方向的分量,即U=ξxu+ξyv, V=ηxu+ηyv;ξt=-xtξx-ytξy,ηt=-xtηx-ytηy;h为水深,m;z为水位,m;q为单位面积上水流的源汇强度,m/s;M=uh, N=vh;q11=ξ $_{x}^{2}$ +ξ $_{y}^{2}$ ,q12=ξxηx+ξyηy;q22=η $_{x}^{2}$ +η $_{y}^{2}$ ;Mξ、Mη、Nξ、Nη、Skξ、Skη表示偏导数;ρ为水体密度,g/cm3;γ′为泥沙干密度,kg/m3;α为恢复饱和系数;ε为泥沙扩散系数;ω为泥沙颗粒沉速,m/s;ωK、SK、S*K分别对应分组泥沙颗粒沉速、含沙量和挟沙力。

本文引用文献[5]的研究成果,即在泥沙连续方程(5)中加入横向输沙项,根据图1得到修正后的泥沙连续性方程如下:

$\begin{array}{l} J \frac{\partial h S_{Κ}}{\partial τ} + \frac{\partial (U + ξ_{t}) h S_{Κ} + Q_{n Κ ξ}}{\partial ξ} + \frac{\partial (V + η_{t}) h S_{Κ} - Q_{n Κ η}}{\partial η} = \\ \frac{\partial}{\partial ξ} [D J (q_{11} h S_{Κ ξ} + q_{12} h S_{Κ η})] + \\ \frac{\partial}{\partial η} [D J (q_{12} h S_{Κ ξ} + q_{22} h S_{Κ η})] - J α W_{Κ} (S_{Κ} - S_{* Κ}) + J q s_{0} (7) \end{array}$

式中:QnKξ、QnKη曲线坐标下沿ξ、η方向的横向输沙率。

1.2 模型辅助方程及相关问题的处理

(1)水流挟沙力公式。

采用与实际吻合较好的张瑞瑾公式:

$S_{*} = k (\frac{\bar{u}^{3}}{g h \bar{ω}})^{m} (8)$

式中: $\bar{u}$ 为垂线平均流速,m/s;g为重力加速度,m/s2;h为水深,m; $\bar{w}$ 为泥沙平均沉速,m/s;k、m为系数。

其中, $S = \sum_{Κ = 1}^{n} S_{Κ}, S_{*} = \sum_{Κ = 1}^{n} S_{* Κ}$

(2)横向输沙率公式。

采用张红武公式[9],具体形式如下:

$Q_{n Κ} = \frac{75.86 \bar{u} h S_{Κ}}{Ν_{0} r} \exp (0.093 1 \frac{ω}{κ u_{*}} \arctan \sqrt{\frac{1}{η} - 1}) f (η) (9)$

其中:

$\begin{array}{l} f (η) = [(1 + 5.75 \frac{g}{C^{2}}) η^{1.875} - 0.88 η^{2.14} + (0.034 - \\ 12.5 \frac{g}{C^{2}}) η^{0.857} + 4.72 \frac{g}{C^{2}} - 0.088] \end{array}$

$Ν_{0} = \int_{0}^{1} η^{\frac{1}{7}} \exp (0.093 1 \frac{ω}{κ u_{*}} \arctan \sqrt{\frac{1}{η} - 1}) d η$

式中:r为流线曲率半径;C为Chezy系数;u*为摩阻流速。

(3)分组挟沙力及挟沙力级配。

先用挟沙力公式求出非均匀悬移质泥沙的总挟沙力,然后用挟沙力级配乘以总挟沙力,得到分组挟沙力公式:

S*K=PK*S*

挟沙力级配与泥沙来源的级配有关。水中泥沙有两个来源,一是上游来沙,而是床面补给,挟沙力级配采用如下公式:

$Ρ_{Κ} = \frac{S_{* Κ}^{´} + S_{Κ}}{\sum_{Κ \geq Κ_{d}} (S_{* Κ}^{´} + S_{Κ})}$

(4)床沙级配采用韦直林模式。

将河床组成概化为表层的交换层、中层的过渡层及底层的泥沙冲刷极限层。规定在每一个计算时段内,各层间的界面都固定不变,泥沙交换限制在表层进行,中层和底层暂时不受影响。在时段末,根据床面的冲刷或淤积往下或往上移动表层和中层,保持这两层的厚度不变,而令底层厚度随冲淤厚度的大小而变化。

1.3 方程求解

模型采用水流泥沙非耦合方法求解,计算物理量布置于同位(非交错)网格,微分方程的数值离散采用有限体积法(控制容积法),离散方程的求解基于Pantankar压力校正法原理[10],具体计算步骤如下。①二维水流计算。②根据水深判断各断面是否属于滩岸,并对岸坡进行崩塌稳定性分析。计算土体崩塌结果,包括单位堤防崩塌土体的宽度及崩塌体积。横向移动岸坡崩塌处的计算网格,计算下一步时泥沙侧向输入项及河床地形高程。③悬移质泥沙扩散方程的求解。④根据河床变形方程计算河床变形,并修正河床高程。⑤进行下一时段的计算。

2 算例及结果分析

由于目前缺少河岸冲刷过程以及弯道内悬移质泥沙不平衡输移的实测资料,故以一概化的90°弯道的河床变形作为模拟对象,模拟凹岸的冲刷后退,对结果进行定性分析。

2.1 模型概化

假定有一90°概化弯道,顺直进口段350 m,90°弯道550 m,顺直出口段800 m,河宽100 m。初始河岸高程为8.0 m,出口河底高程为0.0 m,河床初始比降0.000 2,初始糙率0.020。模拟河段共划分网格数为58×21,其中ξ方向主槽平均网格大小约为30 m,η方向主槽平均网格大小约为5 m。如图2、3所示。边界条件:进口断面给定流量600 m3/s,出口水位7.0 m。

计算流场如图2。在弯道顺直进口段,纵向平均流速沿河宽基本呈对称分布。进入弯道,凸岸流速稍增。至弯顶以下,出现相反的调整,流速渐渐移向凹岸。出弯水流,在相当长的距离内,最大流速仍靠近外侧河岸。计算结果与张红武[11]弯道模型流速平面分布结论一致。说明本文中所采用的平面二维水流数学模型计算方法及相关处理是可行的。

2.2 二维泥沙计算

床沙平均粒径为0.087 mm,河岸土体的特性为:容重γs为18 kN/m3、内摩擦角Φ为14°、凝聚力C为20 kN/m2,堤岸土体的起动切应力τc为1.2 N/m2。边界条件用恒定流计算。进口断面给定流量Q=600 m3/s ,出口水位Z为7.0 m。进口来流中悬移质含量为0 kg/m3;挟沙力系数k=0.15,m=0.92;恢复饱和系数冲刷时取α=1.0,淤积时取α=0.25。本模型模拟了凸岸边滩形成后,弯道水流和泥沙特性。模型概化计算地形如图4、5所示。

(1)冲淤变化。图6给出典型断面CS1不同时刻形态图。从图中可以看出,由于横向输沙的作用,导致凹岸冲刷后退,凸岸淤积。随着时间的推移,凹岸不断冲刷,导致河床岸坡失稳坍塌,河道展宽。图7(a)～(d)为T=0 d,10 d,20 d,30 d时刻地形云图。

(2)岸坡崩塌。图8(a)～(b)给出始末时刻崩塌变化下的网格移动的相对位置。从图中看出,凹岸在横向环流的冲刷下逐渐失稳坍塌,最大崩塌后退距离发生在顶冲点附近位置。从图9各时刻崩塌后退位置示意图可看出,随着主槽冲深下切及凹岸的冲刷后退,崩塌位置也逐渐移向下游。

上述计算结果与现有的研究结论是吻合的,说明扩展后的模型具有一定的可行性,能够反映弯道环流横向输沙以及由此产生的凹岸的崩塌后退。移动网格技术的应用,解决了崩岸宽度与计算网格的尺寸存在的矛盾,克服了采用固定网格进行计算时难以精确调整模型计算边界的困难。

下一步将在本文的基础上,结合具体河段的实测资料,将修正后的二维水沙模型应用于实际河道,使其能够较好地模拟弯道水流的水沙运动规律及河岸的崩塌现象。

3 结语

(1)采用移动网格技术,能够较好的解决崩岸宽度与计算网格的尺寸存在的矛盾,并能精确的调整模型计算边界。在保证计算速度的前提下,提高计算精度。

(2)扩展后的平面二维水沙模型引入横向输沙项,近似考虑了弯道环流所引起的横向输沙。使之更好的反映实际河道的冲淤变化。

(3)采用修正后的模型,对一概化弯道进行计算,计算结果与现有的认识一致,证明该模型具有一定的可靠性和实用性。

参考文献

[1]N Nagata,T Hosoda,Y Muramoto.Numerical Analysis of River图8网格位移图图9各时刻崩塌后退位置示意图Channel Processes with Bank Erosion[J].Journal of HydraulicEngineering,2000,126(4):243-252.

[2]Chang-Lae Jang,Yasuyuki Shimizu.Numerical Simulation of Rel-atively Wide,Shallow Channels with Erodible Banks[J].Journalof Hydraulic Engineering,2005,131(7):565-575.

[3]杨建明,吴建华.动网格技术数值模拟挑流冲刷过程[J].水动力学研究与进展A辑,2001,16(2):6.

[4]刘玉玲,刘哲.弯道水流数值模拟[J].应用力学学报,2007,24(2).

[5]钟德钰,张红武.考虑环流横向输沙及河岸边形的平面二维扩展数学模型[J].水利学报,2004,(7).

[6]夏军强.河岸冲刷机理研究及数值模拟[D].北京:清华大学,2002.

[7]谢作涛,张小峰.一般曲线坐标系平面二维水沙数学模型研究与应用[J].泥沙研究,2005,(6):12.

[8]袁晶.移动坐标下平面二维水沙数学模型研究及应用[D].武汉:武汉大学,2006.

[9]张红武,吕昕.弯道水力学[M].北京:水利电力出版社,1993.

[10]帕坦卡S V.传热与流体流动的数值计算[M].张政译.北京:科学出版社,1984.

[11]曾庆华.关于弯道的底沙运动问题[J].泥沙研究,1982,(3).

移动坐标篇2

1 不同坐标系统之间坐标转换的基本原理

1.1 四参数相似变换模型

地方独立坐标系统与国家坐标系统的转换, 实际上就转换成为了平面直角坐标系统之间的转换, 当然它是有一定的适用范围的。四参数模型只顾及小范围之内的平面坐标系统之间的转换。如独立坐标系的原点在西安80坐标系统下的坐标为, 独立坐标系统对西安80坐标系统的旋转角为, 尺度比参数为m, 则相似变换模型如下:

公共点在两个坐标系中的坐标差为:

上式即由公共点的坐标差表示的“观测值”方程, 由于“观测值”带有偶然误差, 设观测值的改正数为测值的改正数为, 则观测值的误差方程如下:

每一对公共点可以列立2个误差方程式, 方程中4个未知参数, 至少需要2个公共点, 在多于2个公共点的情况下, 可以利用最小二乘原理, 平差求解四个未知参数的最或然值。最后按下式对所有控制点批量进行坐标转换:

坐标重合点可采用在两个坐标系下均有坐标成果的点。但最终重合点还需根据所确定的转换参数, 计算重合点坐标残差, 根据其残差值的大小来确定, 若残差大于3倍 (或2倍) 中误差则剔除, 重新计算坐标转换参数, 直到满足精度要求为止。

1.2 多项式拟合模型

多项式拟合模型如下, 其中的 (x1, y1) 表示源坐标系, (x2, y2) 表示目标坐标系, ai (i=0, 1…, 5) 和bi (i=0, 1, …, 5) 称为多项式拟合系数, 也是待求的未知参数。采用间接平差原理, 可以求出未知参数的平差值。

2 阜新独立坐标系与西安80坐标系转换方案

2.1 阜新独立坐标系统概述

阜新城建坐标系的中央子午线121°30′, 西安80椭球, 投影至高斯平面, 起算点 (西安80坐标系) :孙家湾, 起算方位:孙家湾->八家子山。阜新城建坐标系的起算数据采用的是西安80坐标系的一点和一方位;与西安80坐标系相比存在微小的旋转和尺度变化, 另外由于平面四参数相似变换模型的四参数存在相关性, 因此也引起了微小的坐标平移;阜新市平均海拔150米左右, 地面水平长度归算至国家规定的椭球面要加改正。

2.2 外业测量方案

采用CORS系统的RTK控制测量技术, 直接采集实验点位的阜新独立坐标系统和西安80坐标。

(1) 架设三脚架整平对中; (2) 分时段测量, 上午、下午或者换一天, 每一时段测量3个测回共20组坐标, 每个测回测量6-7次坐标, 每个坐标采用5次平滑取平均值。测回间观测前, 关闭数据通讯链路, 取下流动站接收机, 等待5分钟后重新开机测量; (3) 测回内和测回之间的坐标较差不超过4cm; (4) 取各时段坐标平均值作为该点坐标值。

3 坐标转换结果及分析

3.1 基于四参数的阜新独立坐标系控制点转换

使用已有的的121°30′坐标与使用求解之后的转换参数转换原中央子午线为123°的国家坐标, 所求的转换四参数和公共点坐标差的改正数V (平差的残差) 如 (表1) 。

明显可以看出, 转换残差V的x/y分量部分超过5cm, 分析原因可能是两套坐标系统的长度变形不一致, 因此考虑将两套坐标系统的投影中央子午线设定为相同的值, 再进行转换 (表2) 。

3.2 基于多项式拟合的阜新独立坐标系控制点转换

采用一次多项式拟合模型, 共计6个未知参数, 选择四个公共点:孙家沟、扣莫、烟台营子、古香园作为计算转换参数的点。多项式拟合系数如下:X拟合三系数 (A0、A1、A2) , Y拟合三系数 (B0、B1、B2) 。公共点坐标差的不符值V如 (表3) :

4 结语

以上通过分析四参数转换模型及多项式拟合模型的基本原理, 本文结合阜新市的具体地理位置进行了阜新独立坐标系统与国家80坐标系统之间的转换研究。无论在城市建成区范围, 还是阜新市辖区范围内, 控制点坐标转换及公共点检核的精度在3.5cm范围之内, 满足生产的技术要求。中心化测区的转换参数虽然与非中心化测区的转换参数不同, 但二者转换的效果是一样的。因此, 本文充分的验证了四参数模型与多项式拟合模型在阜新独立坐标系统建立中的可行性、准确性。同时也为在此类地区的工程项目建设提供了转换依据, 提高了测量工作的效率。但要注意不存在一套可以通用的转换参数, 在具体应用时应根据作业区域的坐标系统情况进行坐标系之间的分析, 确定坐标转换模型、进行坐标转换精度估计, 并按照坐标转换的实施步骤进行。

参考文献

[1]陈俊勇.对我国建立现代大地坐标系统和高程系统的建议[J].测绘通报, 2002 (8) :1-5.

[2]茹仕高, 李倩霞.面向智慧城市的空间坐标系统维持与转换[J].测绘通报, 2015 (2) :19-22.

移动坐标篇3

随着GPS(Global Positioning System,全球卫星定位系统)技术的发展,如今在许多项目建设中应用了GPS技术。道路工程中的控制测量往往距离长、范围广、常规方法施测难度大、工期长,更可发挥GPS不需通视和点距可长可短的特点。因此,在道路工程中,施工平面控制点的建立主要使用了GPS、RTK技术。

利用GPS技术平差后得到的坐标属WGS- 84坐标系,它是一个以地球质心为坐标原点的地心坐标系。而我国使用的坐标系统主要有北京54坐标系、西安80坐标系,以及2008年才实施的国家2000坐标系。为了充分利用已有的测量成果,如地形图、工程图等,目前坐标成果大多使用的是北京54坐标系。北京54坐标系以克拉索夫斯基椭球为参考椭球,采用Gauss投影的方式,以30°或60°带划分整个中国所在区域。为了避免投影误差过大,很多城市还采用了地方独立坐标系,把中央子午线定在城市的中央,投影平面选择为城市平均高程面。这些原因使得我国的平面坐标系统比较复杂,GPS直接测定的结果需要进行适当的坐标转换。

本文将详细介绍GPS定位结果由WGS-84坐标转换为北京54平面坐标系的算法。

2 平面转换

2.1 坐标转换的步骤

(1)搜集测区内的坐标重合点成果(搜集同时具有WGS-84坐标系下坐标值和北京54坐标系下的坐标值、高程值的点);

(2)选择适合的重合点用于计算坐标转换参数;

(3)将重合点的两坐标系坐标换算成同一投影带的高斯平面坐标;

(4)采用平面四参数转换模型,根据确定的转换方法与转换模型利用最小二乘法初步计算四参数坐标转换参数;

(5)分析重合点坐标转换残差,将存在粗差的公共点剔除(假设在测区内有n个已知公共点,通过优化选择,剔除残差大于2倍残差中误差的公共点),然后重新计算转换参数,直到满足一定的精度要求为止;

(6)坐标转换残差满足精度要求时,计算最终的坐标转换参数并估计坐标转换参数精度;

(7)根据转换参数,将转换点的坐标转换到目标坐标系。

2.2 平面四参数坐标转换模型

在两平面直角坐标系之间进行转换,需要有四个转换参数,其中有两个平移参数(Δx0,Δy0)、一个旋转参数α和一个尺度比因子m,转换公式如下:

$\begin{array}{l} [\begin{matrix} x \\ y \end{matrix}]_{54} = [\begin{matrix} Δ x_{0} \\ Δ y_{0} \end{matrix}] + (1 + m) [\begin{matrix} \cos α & \sin α \\ - \sin α & \cos α \end{matrix}] \\ [\begin{matrix} x \\ y \end{matrix}]_{W G S - 84} (1) \end{array}$

通过选择的坐标重合点(WGS-84的高斯投影坐标(x,y)WGS-84和北京-54的坐标值(x,y)54)计算出相应的四个转换参数(Δx0,Δy0)、α、m。进而可以将需要进行坐标转换的所有GPS点坐标值转换为我们需要的北京54坐标。

可以看出上述平面坐标转换的方法有较高精度,原理简单,数值简单可靠,但由于使用模型为线性模型,用到的高斯投影是非线性模型,因此大面积使用时会产生变形。平面坐标转换的方法适合小面积坐标转换的使用。

3 空间转换

GPS测定的部分点中,应有一部分点平面坐标已知,并且大地高已知。若使用空间转换方法,应先将已知的平面坐标(x,y)54进行高斯坐标反算,得到54坐标系统的大地坐标(B,L)54,再加上大地高,转为54椭球的空间坐标(X,Y,Z)54 ,而GPS测得的坐标表示为(X,Y,Z)84

WGS-84坐标向北京54坐标转换的三维转换公式为:

$\begin{array}{l} [\begin{array}{l} X \\ Y \\ Ζ \end{array}]_{54} = [\begin{array}{l} Δ X_{0} \\ Δ Y_{0} \\ Δ Ζ_{0} \end{array}] + (1 + Δ μ) [\begin{matrix} 1 & ε_{z} & - ε_{y} \\ - ε_{z} & 1 & ε_{x} \\ ε_{y} & - ε_{x} & 1 \end{matrix}] \\ [\begin{array}{l} X \\ Y \\ Ζ \end{array}]_{84} (2) \end{array}$

在式(2)中,(ΔX0,ΔY0,ΔZ0)、Δμ、εx、εy、εz分别为坐标转换的三个平移参数、一个尺度参数、三个旋转参数。这是WGS-84坐标系统到北京54坐标系统进行转换的七参数,因此称为七参数模型。为可靠地求定7个转换参数,至少需要三个公共已知点方可求解。

当式(2)中,不求旋转参数和尺度参数时,式(2)变成三参数模型,见式(3)。

$[\begin{array}{l} X \\ Y \\ Ζ \end{array}]_{54} = [\begin{array}{l} Δ X_{0} \\ Δ Y_{0} \\ Δ Ζ_{0} \end{array}] + [\begin{array}{l} X \\ Y \\ Ζ \end{array}]_{84} (3)$

此时,使用三参数模型时,仅需1个公共已知点。

求得这些参数后(七参数或者三参数),可以进行坐标转换得到每个未知点在54坐标系下的空间坐标,再将这些空间坐标转为大地坐标,并进行高斯转换得到平面坐标。

空间转换模型通常在测区很大的时候使用。空间转换模型中,高程的精度对平面坐标精度的影响很小,由于旋转参数和平移参数较强的相关性,当测区面积较小时,七参数和三参数法的差别不大。

4 总结

通过对平面坐标转换及空间坐标转换的七参数转换和三参数转换的比较,我们得到以下结论:

在工程测量中,如果测区范围不太大,可以使用平面坐标转换的方法进行,或者采用空间三参数坐标转换的方法进行,都可得到足够的精度。如果测区范围较大,则要采用空间七参数转换的方法进行才能得到足够的精度。

施工中采用的是水准高,而非大地高,因此可以采取水准测量的方法得到的水准高来代替大地高,或者采用高程拟合的方法,利用部分点的水准高来拟合其他点的水准高。

摘要：GPS测量得到的是WGS-84的地心空间直角坐标,而工程中通常使用地方独立坐标系。两者的转换方法一般分为平面转换模型和空间转换模型。平面转换模型原理简单,数值稳定可靠,但只适用于小范围的CPS测量;空间转换模型可用于大范围的CPS测量,按实际情况又分为七参数转换和三参数转换两种。对这两种转换模型进行介绍。

关键词：坐标系,GPS,平面转换,空间转换

参考文献

[1]王解先,等.WGS-84与北京54坐标的转换问题[J].大地测量与地球重力学.2003(8).

移动坐标篇4

${\begin{cases} X = (Ν + Η) c o s B c o s L \\ Y = (Ν + Η) c o s B s i n L \\ Ζ = [Ν (1 - e^{2}) + Η] s i n B \end{cases}$

(1)

式中, B代表纬度;L代表经度;H代表高度[1];e代表偏心率。

而地心空间直角坐标转化为大地坐标的反解无法通过上式直接解算。反解算法按照计算方式主要分为迭代法和直接法两类[2]。文中就是对基于Bowring思想推导的直接算法进行分析和改进, 并推导出计算高程的一种高精度的解算算法。

后面的计算全部采用WGS84椭球, 其有关参数, 如表1所示[1]。

1 Bowring算法

为了提高计算的速度, 国内外大地测量学者曾以不同的途径, 相继推导出了多种直接解算公式。经过全面的分析和比较, 发现Bowring研究思路导出的一组转换公式以其计算简捷且精度高的特点, 而备受工程人员的推崇。下面先讲解该算法的思路, 而后对其进行分析和改进。

1.1 经度计算算法

经度可以通过式 (2) 直接解算

$L = \arctan \frac{Y}{X}$ (1)

计算得到的结果取值范围在 $(- \frac{π}{2}, \frac{π}{2})$ 区间, 须将该结果按照表2进行换算才能解得经度。

1.2 纬度计算算法

建立子午面直角坐标系[3], 如图1所示。图中外圆是以地球长轴a为半径的辅助圆。

图中A是地球外部空间某一点, 其地心空间直角坐标 (X, Y, Z) , P为A点在椭球面上的垂足点, Q为A点和地心O的连线OA与椭球面的交点。Φ表示地心纬度[1], U表示规划纬度[1], B为大地纬度, r=X2+Y2。设M点为子午圈在P点的曲率中心, M′点为子午圈在Q点的曲率中心。由解析几何得知, 椭圆上某点的参数方程以及M和M′点的坐标分别如式 (3) ～式 (5) 所示。

${\begin{cases} r = a \cos U \\ Ζ = b \sin U \end{cases}$

(3)

${\begin{cases} r_{Μ} = e^{2} a (\frac{r_{Ρ}}{a})^{3} = e^{2} a \cos^{3} U_{Ρ} \\ Ζ_{Μ} = - e^{'}^{2} b (\frac{Ζ_{Ρ}}{b})^{3} = - e^{'}^{2} b \sin U_{Ρ} \end{cases}$

(4)

${\begin{cases} r_{Μ^{'}} = e^{'}^{2} a (\frac{r_{Q}}{a})^{3} = e^{2} a \cos^{3} U_{Q} \\ Ζ_{Μ} = - e^{'}^{2} b (\frac{Ζ_{Q}}{b})^{3} = - e^{'}^{2} b \sin U_{Q} \end{cases}$

(5)

由图1得出

$\tan B = \frac{Ζ_{A} - Ζ_{Μ}}{r_{A} - r_{Μ}} = \frac{Ζ + e^{'}^{2} b \sin^{3} U_{Ρ}}{\sqrt{x^{2} + y^{2} - e^{2} a \cos^{3} U_{Ρ}}}$ (6)

式 (6) [4]中B和UP都是未知量, 所以不能直接求解。然而UQ却是可以直接计算得到的, 于是Bowring提出了用M′点坐标代替M点坐标的想法, 也就是用UQ代替UP计算B。UQ可用式 (7) 直接计算得到。

$\tan U_{Ρ} = \frac{a}{b} \tan Φ = \frac{a z}{b \sqrt{x^{2} + y^{2}}}$ (7)

计算得到UQ后代入式 (8) 估算纬度B。

$\tan B = \frac{z_{Ρ} - z_{Μ}}{r_{Ρ} - r_{Μ}} = \frac{z + e^{'}^{2} b \sin^{3} U_{Q}}{\sqrt{x^{2} + y^{2} - e^{2} a \cos^{3} U_{Q}}}$ (8)

当H=0时, A, P, Q这3点重合, 式 (8) 严格成立;当H>0, 尤其是H>1 000 km时, 式 (8) 计算精度下降。

1.3 高程计算算法

高程计算公式如式 (9) [5]和式 (10) [6]两种形式。

$Η = \frac{Ζ}{s i n B} - Ν (1 - e^{2})$ (9)

$Η = \frac{\sqrt{x^{2} + y^{2}}}{c o s B} - Ν$ (10)

2 算法改进

可见, 经度算法不需要改进, 下面论述对纬度和高程算法的改进。为了使其在H很大时也能保持算法精度的稳定性, 下面进一步对Bowring算法进行改进。

2.1 纬度算法的改进

对于某个固定的纬度B, 当高程H不一样时, 利用式 (8) 估算得到的纬度记为B*。B*与B之间的差值会随着H的增加变换。另一方面, B在[-90, 90]范围变化也会对B*有影响。因此引入一个修正系数q=B*/B, q是B和H的函数。将B和H在取值范围内进行划分, 由于篇幅关系只取了整个划分中的部分数据点, 如表3所示。因为经度L对计算结果没有任何影响, 所以可以任意设置一个L值, 表3是在L=45°时计算得到的。

从计算结果可以看出, 当高程>1 000 km时, 其计算误差也会跟着加大, 这是因为Bowring算法是只有当空间点在椭球面上时才严格成立。为了保证在高精度定位领域中, 解算算法的精度能稳定在一定的范围之内, 而不会随着高程恶化的问题, 下面说明运用二维线性插值法来修正式 (8) 的计算结果。

该方法主要是估算校正系数q。假设待求的纬度为B, 高度为H, 通过查表3可知B和H所属的区间, Bi-1<B<Bi, Hj-1<H<Hj (i表示行号, j表示列号) 。设由 (Bi-1, Hj-1) , (Bi-1, Hj) , (Bi, Hj-1) , (Bi, Hj) 确定的q分别为qi-1, j-1, qi-1, j, qi, j-1, qi, j。由二维线性插值的性质可推出下列公式

$q_{j - 1} = \frac{(q_{i, j - 1} - q_{i - 1, j - 1}) B + B_{i} q_{i - 1, j - 1} - B_{i - 1} q_{i, j - 1}}{B_{i} - B_{i - 1}}$ (11)

$q_{j} = \frac{(q_{i, j} - q_{i - 1, j}) B + B_{i} q_{i - 1, j} - B_{i - 1} q_{i, j}}{B_{i} - B_{i - 1}}$ (12)

$q = \frac{(q_{j} - q_{j - 1}) Η + Η_{j} q_{j - 1} - Η_{j - 1} q_{j}}{Η_{j} - Η_{j - 1}}$ (13)

上面的推导过程是在假设纬度B和高程H都已知的条件下进行的, 之后这两个量为待求的未知量, 所以在实际运用中需要先用式 (8) 估算纬度, 得到估值B*, 而后代入式 (9) 或式 (10) 估算高程, 得到估值H*。用B*和H*分别代替上面的B和H, 完成查表和线性插值等运算。最后可以计算得到q的估值q*, 代入式 (14) 即可得到校正后的B值。

B=B*/q* (14)

2.2 高程算法的改进

从前面介绍的高程计算式 (9) 和式 (10) 的分母分别包含sinB和cosB可以看出, 当B→0°时, 式 (9) 不再适用, 而当B→90°时, 式 (10) 不适用。为此对式 (9) 和式 (10) 稍加变形和运算, 以得到一个通用公式。

将式 (9) 和式 (10) 分别乘以sin2B和cos2B, 再将两式相加, 整理后得

$Η = \sqrt{x^{2} + y^{2}} \cos B + z \sin B - a \sqrt{1 - e^{2} \sin^{2} B}$ (15)

利用式 (15) 对B求偏导数, 如式 (16) 所示

$\frac{\partial Η}{\partial B} ≜ - \sqrt{x^{2} + y^{2}} \sin B + z \cos B - e^{2} Ν \cdot \sin B \cdot \cos B$ (16)

根据式 (1) 可推导出如下结果

$\sqrt{x^{2} + y^{2}} = (Ν + Η) \cos B$ (17)

将式 (17) 和公式z=[N (1-e2) +H]sinB代入式 (16) 得 $\frac{\partial Η}{\partial B} = 0$ 。

说明式 (16) 不受纬度B计算误差的一阶无穷小项的影响, 且该公式在B的整个取值范围上都适用, 算法稳定性好。

3 仿真计算

文中的仿真计算是对式 (8) 计算得到的纬度B和修正后的结果进行比较, 以验证该二维插值算法的有效性。

由表4中数据可以看出, 直接用式 (8) 计算得到的纬度, 其误差会随着高程的增加明显增加。而经修正后的纬度值其误差被有效地控制在10-5量级以下, 而且不会随着高程增加而恶化, 保证了计算算法的稳定性。但是该二维插值算法会在每次换算中增加12次加法、9次乘法和3次除法, 另外还有查表操作, 给运算带来额外的负担, 所以该算法应作为Bowring算法的补充, 工程人员可在实现算法时对用Bowring算法估算的高程结果进行判断, 当高程>1 000 km时才采用该算法以保证算法稳定性。

4 结束语

文中重点是讲解基于Bowring思路的直接算法, 以及针对该直接算法的修正算法。分析了Bowring算法在高程>1 000 km时, 其计算误差随高程增加而加大, 提出了用二维线性插值的方法来提高该算法的精度和稳定性, 并通过仿真给予验证, 该修正算法可作为对Bowring算法的补充, 建议在高程>100 km时使用。另外, 文中提出的高程计算公式对B的整个取值范围都适用, 且其计算精度受纬度计算误差的影响小, 具有计算简单、适用范围广、精度高的特点。

参考文献

[1]边少锋.大地坐标系与大地基准[M].北京:国防工业出版社, 2005.

[2]徐绍铨.大地测量学[M].武汉:武汉测绘科技大学出版社, 1996.

[3]Bowring B R.Transformation From Spatial to Geographical Coordinates[J].Survey Review, 1976 (23) :323-327.

[4]束蝉方, 李斐, 沈飞.空间直角坐标向大地坐标转换的新算法[J].武汉大学学报:信息科学版, 2009 (5) :561-563.

[5]崔永俊.空间直角坐标与大地坐标之间的变换方法研究[J].华北工学院学报, 2003 (1) :73-75.

【移动坐标】推荐阅读：

大地坐标10-14

坐标测量06-05

坐标平移06-11

图像坐标07-05

坐标投影07-30