由GPS基线向量解算地面形变

收稿日期:2004-04-29作者简介:王解先(1963-),男,江苏常州人,教授,博士生导师,工学博士.E 2mail :jiexian @由GPS 基线向量解算地面形变王解先1,2(1.同济大学测量与国土信息工程系,上海　200092;2.现代工程测量国家测绘局重点实验室,上海　200092)摘要:通常采用比较各期坐标变化的方法来分析全球定位系统(G PS )形变网的数据.提出一种整体平差模型,将各期G PS 基线向量放在一起平差,定义各点在任意时刻的坐标为其在参考时刻的坐标和形变速度与形变时间的乘积之和.为了更具直观性、更容易加上先验已知信息,将空间坐标参数转换为平面坐标和大地高参数,将空间坐标形变速度转换为水平位移和沉降速度参数.整体平差同时计算出参考时刻的站坐标和形变速度,还推导了相应的协方差矩阵,算例表明模型正确.关键词:全球定位系统;沉降;形变;整体平差中图分类号:P 226 文献标识码:A 文章编号:0253-374X (2005)07-0967-04Application of GPS Technique on Surface Deformation MonitoringWANG Jie 2xian1,2(1.Depantment of Surveying and G eoin formatics ,T ongji University ,Shanghai 200092,China ;2.K ey Laboratory ofM odern Engineering Surveying ,S tete Bureau of Surveying M apping ,Shanghai 200092,China )Abstract :Usually the deformation velocities are determined by the com paris on of coordinate changes am ong different epochs.In this paper ,a whole adjustment m odel is introduced in processing G PS deformation m onitor 2ing netw ork.Baselines of all epochs are adjusted together.The coordinates of each point are divided to coordi 2nates at reference time and multiplication of deformation velocity and duration interval.T o be m ore clearly the space coordinate parameters are replaced by plan coordinate and ellips oid height ,and the space deformation ve 2locities are replaced by horizontal and vertical deformation velocities.The coordinates of each station at refer 2ence epoch are calculated simultaneously with deformation velocities.The variance and co 2variance matrix is deduced als o.The exam ple shows the m odel w orks well.K ey words :global positioning system (G PS );subsidence ;deformation ;whole adjustment m odel 通常采用水准测量方法来精确测定地面沉降,但水准测量工作量大,对于大型城市,1mm ・km -1的精密水准观测误差,其累积仍然较大.目前G PS (global positioning system )测量精度已经达到较高的水平,而且G PS 观测简单,有很多国际G PS 服务组织(International G PS Service ,IG S )站的资料可以随时取用,与观测数据一起解算,并可以直接得到三维形变数据.因此G PS 越来越多地用在城市形变观测中.通常做法是,在监测区域内布设一个G PS 网,每隔一段时间观测一期,每期数据独立平差第33卷第7期2005年7月同济大学学报(自然科学版)JOURNA L OF T ONG J I UNI VERSITY (NAT URA L SCIE NCE )V ol.33N o.7　Jul.2005得到各点的三维坐标,比较各点各期坐标的变化,得到水平位移和沉降量.笔者试图将各期G PS基线结果采用统一模型,以参考某时刻的三维坐标和形变速度为参数,整体平差直接得到各点的三维坐标和三维形变速度.1　数据处理模型1.1　误差方程以ΔR ij=(ΔX ij　ΔY ij　ΔZ ij)T表示点i至j在t时刻测得的基线k,t以儒略日为单位,误差方程为εk+ΔR ij=R j-R i(1)式中:R i,R j是t时刻i和j点的空间坐标;εk是基线残差.若以R im表示i点在整个网观测中间时刻t m时的坐标,v XYZi表示i点的三维位移速度:v XYZi=v Xi v Yiv Zi 则,R i=R im+v XYZi d tR j=R jm+v XYZj d t(2)式中:d t=t-t mΔ,t m=t l+t f2,Δ=t l-t f2,t f,t l是该网最早观测和最后观测的时间,以儒略日为单位. v XYZi的单位是m・(Δd)-1,d t的取值范围为[-1, 1],有利于法方程的性能.将式(2)代入式(1),误差方程变为εk+ΔR ij=R jm+v j d t-R im-v i d t(3) 对各点在t m时刻的坐标取近似值R0im,将速度取近似值v0i,误差方程写为εk =-EδXiδYiδZi+EδXjδYjδZj-E d tδVXiδVYiδVZi+E d tδVXjδVYjδVZj-l k(4)式中:E为三维单位阵;(δX i　δY i　δZ i)T是参考时刻坐标的改正数;(δv Xi　δv Yi　δv Zi)T是形变速度参数的改正数,常数项为l k=ΔR ij+R0im-R0jm+v0i d t-v0j d t(5)式(5)的权取为G PS基线向量解算得到的协方差阵的逆阵p k.形变速度通常表示为上下(沉降)、东向和北向的值v UE N i=(v U i　v E i　v N i)T,它与空间直角坐标速度v XYZi的关系为v Xiv Yiv Zi=R3(-L i)R2(Bi)v U iv E iv N i=R iv U iv E iv N i式中:B i,L i为该点的经纬度;R3(-L i)与R2(B i)是旋转矩阵.R3(-L)=cos(-L)sin(-L)0-sin(-L)cos(-L)0001,R2(B)=cos B0-sin B010sin B0cos B(6) 测站坐标表示为平面坐标和大地高更为直观,更容易利用平面已知坐标等信息,按文献[1]δXδYδZ=9(X Y Z)9(X g　Y g　h)δXgδYgδh=DδXgδYgδh(7)式中:(X g　Y g　h)T是平面坐标和大地高,由空间坐标化为大地坐标,再由高斯投影求得;D是由数值导数方法求出的9(X Y Z)9(X g　Y g　h)矩阵[1].误差方程(4)变为εk=-D iδXg iδYg iδhi+D jδXg jδYg jδhj- d t R iδvU iδvE iδvN i+d t R jδvU jδvE jδvN j-l k(8)　基线k对法方程的贡献为869 同济大学学报(自然科学版)第33卷　N k =0000000000D Ti p k D i0-D Ti p k D j0d t D T i p k R i0-d t D Ti p k R j00000000000-D Tj p k D i0D Tj p k D j0-d t D Tj p k R i0d t D T j p k R j000000000d t R T i p k D i0-d t R Ti p k D j0d t 2R T i p k R i0-d t 2R Ti p k R j0000000000-d t R Tj p k D i0d t R T j p k D j0-d t 2R Tj p k R i0d t 2R T j p k R j00C k =-D T i p k l kD Tj p k l k-d t R T i p k l k0d t R T j p k l k 总法方程写为:N k δx =C k1.2　起算条件若将k 点在t R 时刻(儒略日)的坐标固定为R k ,相当于加上条件:R 0km+δX kδY k δZ k+(t R -t m )Δv X k v Y k v Z k=Rk即R 0km +D k δX g kδY g k δh k+(tR -t m )ΔR k v 0U k v 0E k v 0N k +(t R -t m )ΔR k δv U kδv E kδv N k=R k(9) (1)固定i点的沉降速度值,相当于不解δv U i .(2)固定i 点的水平位移速度值,相当于不解δv E i 和δvN i .(3)固定i 点平面位置相当于不解δX g i 与δY g i .(4)固定i 点的高程相当于不解δh i .在实际使用中,这些基准条件往往不能准确已知,如已知的坐标或速度值存在误差,处理中,可以采用拟稳平差或对条件方程加权等方法处理[2].1.3　协因数阵按式(8)组成法方程求解,迭代至收敛,得到各点在t m 时刻的站坐标(X g i Y g i h i )T 和单位为Δd 的站漂移速度v Δi =(v U i v E i v N i )T ,及其协因数阵Q ′=Q 11Q 12Q T12Q 22(10)式中:Q 11与Q 22分别对应于求解的站坐标和速度;Q 12为协方差.各站t m 时刻平面坐标与大地高(X g i Y g i h i )T 很容易换算成t m 时刻的空间坐标R im =(X i Y i Z i )T ,R im 与V Δi 的协因数为Q RRQ Rv Q vRQ vv=PQ ′P T(11)其中P =D -11ωD -1n EωE式中:n 为总点数.若站坐标参考时刻改为t 0,速度单位改为m ・年-1,分别以R 0和v 表示,则R 0v=E t 0-t mΔE365ΔER m v Δ(12)其协因数阵Q 为969　第7期王解先:由G PS 基线向量解算地面形变Q =Et 0-t mΔE365ΔEQ RR QRv Q RvTQvvE 0t 0-t mΔE365ΔE=Q RR +t 0-t mΔQ TRv +t 0-t mΔQ TRv +t 0-t mΔ2Q vv365ΔQ Rv +365Δt 0-t mΔQvv365ΔQ I Rv +365Δt 0-t m ΔQ vv365Δ2Q vv(13)1.4　基线残差与验后方差i 至j 在t 时刻的基线向量估值ΔR ^ij 为ΔR ^ij =R ^jm +R j v ^j d t -R ^im -R i v ^i d t(14) 基线残差为εk =-ΔR ij +ΔR ^ij(15)验后误差协因数阵为Q ij =AQAT(16)式中:系数A 的行数为3,列数为总点数乘6,在对应于i 和j 的坐标处为三维单位阵,在对应于i 和j 的速度处为三维单位阵乘d t.A =(-E E -d t R i d t R j )(17) 基线长度的中误差为σ2k =σ2X j -X iρij Y j -Y iρijZ j -Z iρijQ ij ・X j -X i ρijY j -Y i ρij Z j -Z iρijT(18)式中:ρij 为用点i 和j 的坐标计算得到的距离.式(17),(18)可以用于基线剔除.2　算例某大城市采用G PS 技术监测城市形变,全网共有50个点,观测了4期,8台双频G PS 接收机同时观测,数据采样间隔为30s ,观测时段长度为8h ,每期约观测7d.计算时,引入IG S 的佘山(SH AO )站数据一起处理,采用IG S 精密星历,基线处理软件为G AMIT.G AMIT 的基线结果采用以上模型处理,计算时,将其中一个IG S 点的平面坐标固定,将其中一个基岩点的沉降速度固定为零,另外将几个有水准沉降资料的点的沉降速度固定.平差后基线残差的均值为(0.0001　-0.0001　0.0000)m ,水平位移与沉降速度结果也很合理,说明模型正确.3　结论G PS 已经在城市形变方面得到较广泛的应用,笔者提出的整体平差模型,将各期G PS 基线放在一起平差,平差参数为平面坐标、大地高和沉降速度,使得已知条件很容易加,先将坐标参考时刻定在观测中间时刻,将形变速度的单位定为m ・(Δd )-1,使得法方程性能良好,然后将坐标参考时刻从中间时刻改化为任意时刻,并将速度单位改化为惯用的m ・年-1,并推导了方差协方差阵,实际数据计算结果表明,模型正确.参考文献:[1]　W ang J ,Iz H B.A practical G PS netw ork adjustment method[J ].Sur 2vey Review ,1999,35(272):127-133.[2]　樊功瑜.误差理论与测量平差[M].上海,同济大学出版社,1998. FAN G ong 2yu.Error theory and surveying adjustment [M].Shanghai :T ongji University Press ,1998.(编辑:张　弘)79 同济大学学报(自然科学版)第33卷　。