实时相位（通用3篇）
实时相位 篇1
0引言
连续波测控体制具有高精度特点。为了实现远作用距离, 必须提高连续波发射功率。而高可靠性要求发射机尽量采用固态功放模块, 一般采用功率合成技术以获得较大规模的发射功率。
某X波段连续波发射系统要求具备高可靠性和机动性。天线采用平板阵列天线以降低设备的体积重量。为提高发射效率, 减少发射馈线的损耗, 发射机要求尽可能靠近天线单元。此外, 由于波长较短, 设备在加工、使用中会出现相位不一致的情况。如发射机工作时, 温度不均匀性造成各组件的相位不一致, 由于器件的制造工艺影响, 系统相位实际测试时, 会出现变化规律相反的现象。采用自由空间功率合成, 可以满足系统要求。
1自由空间功率合成技术
1.1自由空间功率合成原理
采用自由空间功率合成技术, 合成方式类似于有源相控阵天线的模式。自由空间功率合成原理如图1所示。
对于Nx×Ny的相似天线组成的平面阵, 可以把平面阵分成Nx个Y方向的直线阵 (Y阵) , 则此平面阵可看成以Y阵为天线元的Nx元x方向的直线阵, 则Nx×Ny平面阵的方向函数为:
f (θ, ϕ) =f1 (θ, ϕ) ·fNx (θ, ϕ) ·fNy (θ, ϕ) 。 (1)
式中, f1 (θ, ϕ) 为单个天线的方向函数。
对于均匀 (间隔D) 分布的直线天线阵, 假设每个单元发射功率、天线增益都相同, 则在远端P (r, θ, ϕ) 点的合成场强EΣ为:
式中,
因P点位于天线阵的远区, 即r>>nD, 故可认为各天线元至P点的射线是相互平行的, Nx个辐射场强极化单位矢量相同, 即为θ, 则在远场条件下, 平面阵等效辐射功率ERP。
由于采用等幅全分布体制, 因此在天线阵正前方 (θ=0°) 等幅同相的天线阵对应着远场叠加合成的最大场强, 此即为远场控制算法的基本依据。
1.2合成效率及其影响仿真
各单元相位控制精度直接影响了天线阵合成效率, 当各单元同相位发射信号时, 单元数越多, 远场ERP值就越大。单元数与ERP值的关系示意如图2所示。
空间功率合成效率 (合成效率) 为:
EIRP=PGNxNyη, (3)
式中, Δθ21、…、ΔθNx1、Δθ22、…、ΔθNxNy分别为第2个、第3个、…、第NxNy个单元天线相对于第1个单元天线辐射信号到达接收点的相位差。显然η≤1, 当NxNy个信号同相叠加, 可获得最大的空间合成功率。
由式 (3) 、式 (4) 和图2可见, 等效功率与阵元数的平方成正比, 随着阵元数的增加, 等效功率将迅速增大。
相位差的变化直接影响了功率合成的效率。当相位控制精度在20°以内, 其合成效率能达到90%以上。当相位控制精度在10°以内, 其合成效率能达到97%以上。
此外, 组成天线阵列的单元数也影响着合成效率。当各单元相位差的均方根值一定时, 由于存在各项异同性, 天线阵越庞大, 各单元出现相位差抵消的几率也越大, 其合成效率将有所下降。
同时, 天线阵互藕、温度升高导致发射机功率下降等影响, 大规模天线阵列的功率合成效率较小阵列要低。
2X波段空间功率合成方案
2.1系统组成
为了克服相位变化对合成效率的影响, 采用实时相位校准方案构建验证系统, 其组成框图如图3所示。
验证系统由4个单元组成, 每个单元发射功率约为5 W左右, 天线增益约为24.4 dB。
信号源产生的X频段信号等幅度、等相位分布到各发射模块。经过移相放大后, 等相位地发射到自由空间。
在发射端与天线之间加装定向耦合器, 可随时获得各通道的发射信号。通过开关选择, 可实现对每一发射通道信号的选取, 作为待修正信号, 传送至幅相校准单元, 进行幅度、相位进行测量与校准。
信号源输出一路信号, 作为发射信号的校准参考信号进行比对。
移相器是实现相位校准的修正部件, 该方案采用6位移相器。
控制器根据幅相校准单元工作流程, 控制功放单元、多路开关选通和移相器修正等。
2.2实时相位校准
为保证合成功率最大, 各通道输出相位控制必须保证一致。由于功放为发热部件, 其通道受热后容易引起相位的温度漂移变化。而天线阵组成, 特别是大型天线阵的热分布是不平衡的, 各通道相位也随着温度变化而各不相同, 需要及时进行相位校准。
如图3所示, 双路AD同时对校准参考信号与待修正信号进行采集, 在FPGA进行FFT运算。按式 (5) 计算出信号的相对相位:
式中, Im、Rm分别为所求多普勒频率的虚部值与实部值。
将待修正信号的相对相位与校准信号的相对相位进行比较, 即可得到该路信号的相位差。分别选择不同发射支路作为待修正信号进行幅相校准, 即可完成全阵列单元的幅度相位测量。
通过控制各通道移相器, 即可实现对系统的相位校准。在系统工作间隙, 可实现实时的相位校准。
2.3相位校准误差分析
造成发射信号相位误差的主要有:天线单元相位误差 (包括馈线误差) 、信道相位不一致性测量误差、信道相位不一致修正误差和检测信道不一致性误差等, 各误差源独立不相关。
2.3.1 天馈误差
天线单元间相位误差可达40°, 可通过远场测量获得相位误差。在整机调试中, 通过修正馈线相位的方法进行天线相位修正。其修正误差约为10°。
天线在温度变化和机械应力等条件下还会出现随机的相位漂移, 其误差约为2°左右。
2.3.2 信道不一致性误差
信道不一致性误差是相位修正的主要误差。主要包括各无源微波器件本身故有的相位不一致性误差、功放相位漂移误差、移相器精度引起的修正误差和微波链路相位误差等。
无源器件与链路相位误差属于系统误差, 可以测量并修正。有源器件往往会产生一定的相位抖动, 该相位误差属于随机误差, 约为3°左右。
功放相位漂移误差变化过程相对较慢, 可以通过定期测量得到其变化。一般功放稳定前相位变化较大, 而且各器件见相位变化趋势也不相同, 此时可采用缩短检测间隔时间, 实时对各通道进行检测方式, 获取相位变化的细节信息。
相位测量主要采用比对方法, 根据仿真结果, 当校准信号信噪比在30 dB以上, 采用FFT计算方法, 相位差计算精度在5°以内。
X频段移相器精度一般为4° (均方根误差) 以内, 属于随机误差, 在不同修正值条件下误差各不相同。系统中采用6位移相器, 其量化误差精度为5.625°。
2.3.3 检测电路误差
相位测量除了计算误差外, 还包括检测电路误差。检测电路误差主要包括耦合器、馈线和多路开关等引起的相位误差, 各误差分配如表1所示。
而中频幅相校准单元产生的相位误差, 由于各通道分时共用同一信道, 只需要考虑不同时刻的相位漂移即可, 在50 ms间隔内, 该误差较小, 可按1°考虑。
耦合器为无源器件可计入链路误差, 该误差可通过仪器测量并校准, 最大误差按5°考虑。
分路器误差包括各端口相位不一致性误差 (约为3°) 和开关造成的相位突变。此外在大规模天线阵, 其开关及分路器需要级联, 误差会发生累计。一般开关引起的相位误差约为10°左右。
由表1可见, 总的校准误差可以控制在20°以内, 实际系统测试结果与分析相近。
3功率合成效率测试
3.1测试场要求
由空间功率合成的原理可知, 相控阵天线阵等效辐射功率的测试必须在空间辐射场进行。因此, 要准确测量其等效辐射功率, 测试场地必须满足自由空间传播条件或者可知地面反射等因素对测量结果的影响。
满足或近似满足自由空间传播条件的测试场地称为自由空间测试场。理想的自由空间测试场为电波暗室。其中, 发射天线和接收天线的架设高度应保持一致。测试场必须满足电波传播的远场条件。因此。发射天线和接收天线的距离应同时满足式 (6) 和式 (7) 的要求:
R>2 (d+D) 2/λ, (6)
R>10λ。 (7)
式中, R为发射天线至接收天线的距离 (m) ;D为发射天线的尺寸 (m) ;d为接收天线的尺寸 (m) ;λ为波长 (m) 。
实际应用时, 由于R、D往往较大, 一般电波暗室无法满足测试需要。通常采用高架场或斜式测试场。
3.2测试方法
采用通过引入已知的参考系统测试天线阵等效辐射功率的比较方法进行ERP值测试。再通过测试结果反算出功率合成效率。其等效辐射功率计算方法为:
ERP=Es-Ex+Ps+Gx-Ax-Gs+As。
式中, Es为接收被测系统发射信号电平 (dBm) ;Ex为接收参考系统发射的信号电平 (dBm) ;Ps为参考系统输出功率 (dBW) ;Gx为被测系统发射天线增益 (dBi) ;Gs为参考系统天线线增益 (dBi) ;As为被测系统发射时测试系统的衰减值 (dB) ;Ax为参考系统发射时测试系统的衰减值 (dB) 。
测试中, 衰减值可以提前校准, 并扣除。此外, 为消除布站造成的方向图对准、测试距离变化、接插头操作失误等对测量精度的影响, 采用被测系统中的一个单元作为参考系统。
获得的ERP值与理论计算得到的ERP值进行比较, 即获得系统的合成效率。
3.3测试结果
采用上节方案, 对2*2阵列功率合成效率进行了测试。先对各功放单元输出功率、天线增益进行测试, 测试结果如表2所示。
采用远场测量单路A1发射单元的Ex值作为参考接收电平为-17.9 dBm。
进行A1、B1、C1、D1的合成, 与相位修正。相位修正值如表2所示。其中, 天线相位偏差已通过远场修正。测量功率合成后的Es值为-6.38 dBm。
理论上合成功率应为18.87 W, 4个天线与4个功放所达到的ERP值, 与单路A1接收到Ex值相比较, 其增益应为11.83 dB。
实际增益测量得到的增益为:
-6.38- (-17.9) =11.52 dB。
假设天线合成效率为100%, 理论上, 合成的天线阵增益为30.42 dB。合成所获得的功率增益为5.5 dB, 其合成功率约为17.74 W;合成效率约94%。
实际上天线合成效率不可能达到100%。相位修正的影响, 同时作用于天线与功放, 天线增益合成效率应与功率合成效率相同。
根据以上测试结果, 在相位控制精度小于10°时, 空间功率的合成效率应为96%。与理论分析结果相吻合。
4结束语
某X波段发射设备采用12组发射模块, 在静态测量中, 各单元相位最大偏差超过100°, 如不进行修正, 必然会影响合成效率。
该系统采用实时相位校准方法, 在X波段实现了合成100 W以上的连续波功率发射, 合成效率达到90%以上。其技术已应用于某靶场X波段测量系统, 效果良好。
摘要：介绍了自由空间功率合成技术基本原理和主要技术, 重点分析了X波段功率合成的难点以及影响合成效率的主要因素, 并对合成单元数量及各单元相位偏差对功率合成的影响进行了仿真。在此基础上, 提出一种基于实时相位校准技术的X波段连续波功率合成的方案。采用空间辐射功率的比对测试方法, 对自由空间功率合成效率、合成功率进行了测试, 并给出测试结果, 其最大合成效率达到90%以上。
关键词：空间功率合成,连续波,相位校准,合成效率
参考文献
[1]王芳.Ka频段功率合成放大技术研究[D].成都:电子科技大学, 2005:2-4.
[2]张凤伟, 田智辉, 刘苏, 等.效辐射功率测试方法研究[J].无线电工程, 2006, 36 (4) :47-49.
[3]江志浩, 蔡德荣, 王孜.空间功率合成技术的合成效率问题研究[J].工程实践及应用技术, 2008 (2) :54-56.
实时相位 篇2
实时全息干涉法是全息干涉技术中的主要方法之一。它能够对任意形状、任意材料和粗糙表面的三维物体进行测量,具有全场、实时、非接触、条纹对比度好、测量精度高等特点。通过实时地观察干涉条纹的变化,可以定性地分析出物体表面的形变及其应力等问题[1]。若要自动、定量地进行测量,关键是测量出反映物光场变化的干涉条纹在不同时刻的相位分布,目前对干涉条纹进行自动相位测量的方法有傅里叶变换法、相移法等[2,3]。其中傅里叶变换法适合于单幅条纹的相位分析,只能得到一幅图像中的相对变形量,不适合大尺度物体变形场的测量,而且所提取的相位值为相对相位值(即视场中一点相对于另一点的相位值),不能得出测试过程中该点的相位变化。而相移法需要参考光路中加入相移器,且需要进行精确的控制,一般情况难以实现。这些方法一般都只适合于静态干涉条纹图相位的提取。
希尔伯特变换以其适合于大范围的测量、计算时间短、高分辨力、算法简单等特点在时域信号移相、相位提取等方面得到了广泛的应用[4,5,6]。
本文利用希尔伯特变换法在时间域内对摄像机拍摄的实时全息条纹进行处理,得到视场中每一点在整个测量过程中的相位变化,从而得到全场的相位变化分布情况。
1 实时全息干涉条纹
被测物体在静止状态下完成全息图记录,全息干版H在经过冲洗处理后精确复位。如图1所示,在参考光R与物光O同时照射下,在再现光路中既有由参考光再现产生的包含原物光O′再现像的再现光场UR1,同时又有由物体直接透过全息图形成的实时光场UO0,这两个光场之间具有很好的相干性,形成实时干涉条纹。如果物体发生变化,原参考光再现的原来的物光与变化后的物光叠加的结果,将引起干涉条纹的实时变化,可表示为
式中:φ(x,y,t)=Δφ(x,y,t)-π/2,J为贝塞尔函数,K为与平均曝光量、记录材料等有关系数。当实验条件确定时,系数K、J均为常数[3]。如果物光强O和再现参考光强R近似保持不变,则叠加干涉条纹强度的变化就只取决于物光的相位变化Δφ。如果能够测量出条纹的实时相位φ,就可以根据它与引起相位变化的物理量之间的关系,求出该物理量的变化[1]。
从式(1)可以看出,实时全息干涉条纹是一种动态条纹,在某一点(即x,y确定)的光强将随时间t呈余弦关系变化,因此我们可以在时域内对实时全息干涉条纹进行处理。
2 基于希尔伯特变换的相位提取
希尔伯特变换是完全在时域中进行的一种特殊的正交变换。函数u(t)的希尔伯特变换定义为
其中Hi{u(t)}表示对u(t)作希尔伯特变换,是u(t)的线性函数,它可看作是函数u(t)与(-πt)-1的卷积:
希尔伯特变换等同于一个滤波器,一个信号经希尔伯特变换后,其幅度频谱、功率(能量)谱以及自相关函数和功率(能量)均不变。在频域内,希尔特变换的关系可表示为
式中:V(f)是函数v(t)的傅里叶变换,U(f)是函数u(t)的傅里叶变换,j是虚数符号。由式(4)可知,求一个实函数u(t)的希尔伯特变换,就是将u(t)的所有正频率成份相位移动-π/2,对所有负频率成份相位移动π/2。因此,对f(φ(t))=cos(φ(t))作希尔伯特变换得到sin(φ(t)),利用这一特性,可以得到信号的实时相位:
在Matlab信号处理工具箱中提供了现成的希尔伯特变换函数[7]。
由式(5)直接得到的相位值是随时间变化从-π/2到π/2周期性变化的实时相位,则从i时刻到p时刻的相位变化为
式中:n为条纹变化的级数,即φ(t)从π/2跳变到-π/2的次数;φ(p)和φ(i)分别表示p时刻和i时刻的实时相位。
3 实验与结果分析
本文对模拟动态条纹和实测实时全息干涉条纹进行了相位提取。在模拟实验中,视场大小为200 pixels×200 pixels,初始时刻(t=0 s)物光相位变化为0,物光的相位在41 s内按高斯曲面线性增加到如图2所示的相位分布。初始时刻视场无条纹,如图3(a)所示。在相位增加过程中,视场中间有干涉条纹冒出并不断向外移动,在t=20.48 s、t=40.96 s时刻干涉条纹如图3(b)、(c)所示,点(70,100)的灰度随时间的变化的时域信号如图4(a)所示。
为使背景强度在相位计算中带来的影响最小,本文在进行希尔伯特变换之前,先将信号经过一个高通滤波器以便去除背景强度和噪声的影响。从而得到如图4(b)所示的去除背景后的时域信号。再进行希尔伯特变换之后得到图4(c)所示的波形图,可以看出其波形的信号相位与变换前的相位相差90°。
由式(5)进行相位提取,得到如图4(d)所示相位图。此图为相位值从-π/2到π/2周期性变化的实时相位图。在程序中采用逐级计数的方法来确定条纹变化级数n,即:先设定n=0,再将前一时刻的相位值与后一时刻的相位值进行对比,当发现前一时刻的相位值减后一时刻的相位值大于阈值(实验中取2 rad)时,n加1,从而确定两个时刻间条纹变化的级数。利用式(6)可计算出该点在任意两个时刻间发生的相位变化值。
对整个观察范围内的点进行扫描,求出各点的相位变化值。从t=0到t=40.96 s过程中,各点的相位变化三维图如图5所示。从图中可以看出其相位提取结果与标准值完全相符。
通过大量实验发现,当灰度时域信号噪声增加时,相位提取结果也将出现较大误差,如图6所示为信噪比SNR=20时所得实验结果,从图中可以看出提取结果不仅存在很多随机变化的起伏,而且所得结果一般比标准值要偏大。这是因为当存在较大噪声时,程序在检测式(6)中的n(即φ(t)从π/2跳变到-π/2的次数)时,将可能受噪声影响而大于真实值。
在对实测实时全息干涉条纹提取相位时,作者利用数码摄像机拍摄了铝片受力变形时的实时全息干涉条纹的变化过程,以视频文件保存,利用Matlab软件对实时全息干涉条纹进行处理。因为数码摄像机得到的是彩色图像,在Matlab中每一帧彩色图像用3个二维矩阵分别表示红、绿、蓝三种色彩的图像数据,要处理视频文件就需要建立一个四维矩阵。考滤到实时全息干涉实验中一般采用波长为632.8 nm的He-Ne激光器作为光源,在处理时只取其中的红色图像数据,这样只要一个三维矩阵就可表示出视频文件,大大减小数据处理量。从t=0时刻到t=15.6 s的记录过程中,观察到条纹从上部向下移动,图7为t=0、t=6.4 s和t=15.6 s时刻的实时全息干涉条纹图,其尺寸大小为270×250,像素点(75,25)(以图片左下角顶点作为原坐标原点,以下边缘和左边缘分别作为x轴和y轴)的灰度随时间的变化的时域信号如图8所示。
用以上算法求得从t=0到t=6.4 s、t=11.2 s、t=15.6 s过程中,各点的相位变化三维图如图9所示。从图中可以看出:在没有条纹的部分,得到的相位变化值为0;随着时间的变化,相位变化量增加;在条纹边缘部分从现了一些毛刺现像,这是因为在条纹边缘上条纹的移动不均匀、不连续而引起的误差;在条纹的上部相位变化值比下部的相位变化值要小。这与传统的实时全息干涉条纹人工分析结果是一致的:条纹越密的地方,物光相位变化越大,条纹越疏则相位变化越小。
4 结论
利用希尔伯特变换法能有效提取低噪声情况下动态干涉条纹的相位值及其变化,在时域内对实时全息干涉条纹进行处理,可以得到视场中每一点每一时刻的相位值。进而可求出任意两个时刻间各点发生的相位变化,将为实现全自动、高精度实时全息干涉测量提供有效手段。
参考文献
[1]谢锦平,刘冬梅,王凤鹏,等.用实时全息干涉法分析岩石的侧向受力问题[J].激光杂志,2004,25(3):59-60.XIE Jing-ping,LIU Dong-mei,WANG Feng-peng,et al.Analyzing the problem of rock receive the side directions strength by the way of holographic at any time[J].Laser Journal,2004,25(3):59-60.
[2]Creach K,Schmit J.N-point spatial phase measurement techniques for non-destructive testing[J].Optics and Lasers in Engineering(S0143-8166),1996,24(6):365-379.
[3]吕晓旭,张以谟.实时全息干涉计量中时间序列条纹位相的测量方法[J].中国激光,2004,31(2):219-222.LüXiao-xu,ZHANG Yi-mo.A Novel Phase Measurement Method in Real-Time Holographic Interference Time Sequence Fringes[J].Chinese Journal of Lasers,2004,31(2):219-222.
[4]Violeta D M,Hirofumi K.Dynamic electronic speckle pattern interferometry(DESPI)phase analyses with temporal Hilbert transform[J].Optics Express(S1094-4087),2003,11(6):617-623.
[5]吕捷,王鸣,宦海,等.希尔伯特变换条纹分析法及其在非球面镜测量上的应用[J].光学学报,2005,25(6):781-785.LüJie,WANG Ming,HUAN Hai,et al.Fringe Analysis with Hilbert Transform and Its Application in the Measurement of Aspheric Mirror[J].Acta Optica Sinica,2005,25(6):781-785.
[6]高成勇,周灿林.一种基于希尔伯特变换实施相移的相位解调算法[J].光学技术,2004,30(4):508-512.GAO Cheng-yong,ZHOU Can-lin.A kind of calibration algorithm of phase-shifting based on Hilbert transform[J].Optical Technique,2004,30(4):508-512.
实时相位 篇3
北斗卫星定位系统是我国自主建设的全球卫星定位系统,到目前为止已经初步具备了向亚太地区提供实时定位、导航和授时服务的能力[1]。周跳探测是北斗载波相位数据预处理中的重要部分,是确保获得“干净”数据的重要保障。目前已有很多学者对周跳探测方面的问题进行了研究,并提出了很多可以实际应用的解决方案。其中比较普遍的是多项式拟合法[2]、高次差法[3]、电离层残差法[4,5]、小波滤波法[6]等,但是不同的方法适用情况不一样,也不够完善。多项式拟合法对于大周跳的探测效果较好,但其探测效果受到多项式类型、拟合次数以及截断误差等因素的影响;高次差法需要逐级做差,周跳探测无法实时进行,而且在探测的过程中放大了噪声,不适用于探测小周跳甚至近百周的周跳值;电离层残差法推导过程比较严密,解决了周跳解的多值性问题,适用于小周跳的探测,且是目前探测效果最好的方法之一,但需要利用多频载波相位观测数据,当2个载波相位观测值同时发生周跳且周跳之比与频率之比接近时,该方法无法探测周跳[7],同时该方法只针对多频接收机。小波滤波法可以方便地探测到发生周跳历元的突变点,但周跳的具体幅度还需通过其他方式配合探测[8]。
由于多频接收机价格昂贵,所以在实际生产中普遍使用单频GNSS接收机。在应用北斗单频接收机对矿区重点形变区域进行实时变形监测时,为提高北斗单频载波相位观测值预处理中周跳探测的精度和效率,本文提出了一种基于Kalman滤波的单频北斗载波相位实时周跳探测方法,给出了相应的数学模型,论述了周跳与粗差的识别方法,并利用淮北市某矿区沉降变形监测中实测北斗单频载波相位数据对该方法进行了验证。实验结果表明,该方法能准确区分粗差与周跳,同时能较为精确地给出周跳大小,提高了北斗单频载波相位观测数据预处理的精度。
1 周跳探测原理
1.1 Kalman滤波模型
Kalman滤波的数学模型由状态方程和观测方程组成,对于线性系统,其离散模型为[9]
式中:Xk为tk时刻的状态向量;A为tk-1时刻到tk时刻的状态转移矩阵;Wk-1为tk-1时刻的预报模型噪声,其值取决于状态方程表述的好坏,表述得越精确取值越小[10];B为tk-1时刻的动态噪声系数矩阵;φk为tk时刻的载波相位观测向量;C为tk时刻的状态向量系数矩阵;Vk为tk时刻的观测噪声,其方差为σ2,在静态载波相位测量中σ=λ/100[11],即波长的分辨率。
用于周跳探测的Kalman滤波模型是建立在三阶多项式的基础上,记载波相位观测值的采样间隔为T,则
状态向量Xk由φk及其对时间t的一阶导φ′k、二阶导φ″k、 三阶导φ'''k组成, 即
根据最小二乘原理可得Kalman滤波递推步骤:
(1)状态向量的一步预测值为
(2)状态向量一步预测值方差矩阵为
式中:Q为状态噪声的协方差阵,Q=E(W2)。
(3)滤波增益矩阵为
(4)状态向量估计值为
(5)状态向量估计值方差矩阵为
式中:I为单位矩阵。
1.2滤波初始值确定
为了启动Kalman滤波程序,必须要有一个初始值X0和P0,它一般可由前N点通过三阶多项式平滑取得,本文取N=12,但N并不是一个固定数值,N的选取应视观测数据质量情况而定,若N值选取不当,容易造成滤波发散,不能很好地探测周跳。
用N个连续的载波相位观测值进行三阶多项式平滑,求出初始值X0和P0:
式中:L为载波相位观测值向量。
X0的协方差矩阵为
在选取N个点进行三阶多项式平滑时,应确保这N个观测值中不含有周跳,这可以通过残差平方和来检验:
单位权中误差估计值为
检验时通常以2.5倍的中误差作为限差:
当式(12)成立时,便可开始进行滤波计算。
1.3周跳探测具体实现
在进行滤波计算时,先计算得到每个历元的一步预测值估计值与实际观测值之间会存在一定的误差若其不满足预设限差要求,即|V(k)| >u(u为预设的周跳检验阀值),则表明该历元时刻的载波相位观测值含有周跳或粗差,针对此种情况需要做进一步处理。
当|V(k)| >u成立时,需进一步计算第k+1个历元时刻的载波相位观测值与其对应时刻的一步预测值之间的误差并判定V(k+1)与检验阈值u′=2.5u间的关系,u′=2.5u是因为V(k+1)是第k+1个历元时刻的一步预测误差,在计算过程中第k个历元时刻的观测值并没有参加滤波计算。若V(k+1) >u′成立,则表明在第k个历元时刻载波相位观测值发生了周跳,此时应重新进行初始化;若不成立,则表明在第k个历元时刻载波相位观测值中含有野值或者粗差,此时应令第k个历元时刻的增益矩阵Jk=0,即该时刻的观测值不参与滤波计算。
具体的周跳探测流程如图1所示。
2 实验结果分析
2.1 实验设计
实验数据来源:采集于淮北市某矿区,采样间隔为1s,经过其他数据预处理软件检验为不含有周跳的“干净”的北斗单频载波相位观测数据。设计2组实验来验证本文所提出的方法在矿区变形监测数据预处理中的可行性以及周跳探测结果的精确度。实验1分别利用多项式拟合法和Kalman滤波法来检验实验数据的噪声水平,结果如图2所示。
实验2在上述一组干净数据中第100历元处加入大小为1的粗差,在第400,700,1 000,1 300历元处加入大小为2,3,4,5周的周跳后,分别利用多项式拟合法、Kalman滤波法进行周跳探测,结果如图3及表1所示。
2.2 实验结果分析
从图2可以看出,对于采样间隔为1s的北斗单频实测数据,多项式拟合法和Kalman滤波法探测出的无周跳噪声水平相差不多,相比之下Kalman滤波的效果要好一些。在实际观测时,受接收机所处环境的影响,使得原始数据含有较大观测噪声,Kalman滤波预报残差序列在[-0.497,0.461]周范围内波动,在一定程度上会影响该方法探测周跳的精度。
图3直观地反映出Kalman滤波可以更加灵敏、准确、实时地探测出周跳或粗差发生的位置,较好地满足实时变形监测的时效性。从表1中可以清晰地看到,本文所提出的适用于北斗单频载波相位周跳的Kalman滤波法能够准确地区分出滤波曲线发生波动时观测值的变动类型。总的来说,在观测噪声较大时,Kalman滤波依然能够灵敏、精确地探测出周跳发生的位置,准确地区分周跳与粗差,同时也能够较为精确地确定发生周跳的估值大小。
3 结论
针对在矿区变形监测中如何提高北斗单频载波相位观测数据预处理的精度问题,提出了相应的Kalman滤波方法,并利用淮北市某矿区的北斗单频载波相位实测数据,从理论分析与实例验证2个方面证明了该方法的可行性及探测结果的可靠性,并得出以下结论:
(1)提出的Kalman滤波法能较为准确、灵敏地探测和区分出北斗单频载波相位观测值中的粗差或周跳,能近似实时地探测出载波相位观测值中的周跳。
(2)在启动Kalman滤波探测周跳前,需要N个干净的历元数据进行滤波初始化,同样在探测到发生周跳的位置后仍需要进行初始化,这就要求相邻2个发生周跳位置间的历元间隔不得小于N,因此,这在某种程度上限制了Kalman滤波法探测周跳的准实时性。
(3)当观测数据本身带有较大噪声时,周跳探测结果与真实值存在一定误差。因此,如何提高该情况下的周跳探测精度,将是下一步的研究重点。
参考文献
[1]杨元喜,李金龙,王爱兵,等.北斗区域卫星导航系统基本导航定位性能初步评估[J].中国科学:地球科学,2014,44(1):72-81.
[2]蔡诗响,张小红,李星星,等.一种基于多项式拟合的单频周跳探测改进方法[J].测绘信息与工程,2009,34(5):1-3.
[3]DAI Zhen.MATLAB software for GPS cycle-slip processing[J].GPS Solution,2012,16(2):267-272.
[4]CAI Changsheng,LIU Zhizhao,XIA Pengfei,et al.Cycle slip detection and repair for undifferenced GPS observations under high ionospheric activity[J].GPS Solution,2013,17(2):247-260.
[5]王维,王解先,高俊强.GPS周跳探测的方法研究[J].武汉大学学报:信息科学版,2010,35(6):687-690.
[6]蔡昌盛,高井祥.GPS周跳探测及修复的小波变换法[J].武汉大学学报:信息科学版,2007,32(1):39-42.
[7]BANVILLE S,LANGLEY R B.Mitigating the impact of ionospheric cycle slips in GNSS observations[J].Journal of Geodesy,2013,87(2):179-193.
[8]黄兵杰,柳林涛,高光星,等.基于小波变换的GPS精密单点定位中的周跳探测[J].武汉大学学报:信息科学版,2006,31(6):512-515.
[9]GENG Yanrui,WANG Jinling.Adaptive estimation of multiple fading factors in Kalman filter for navigation applications[J].GPS Solution,2008,12(4):273-279.
[10]宋迎春,朱建军,陈正阳.动态定位中测量噪声时间相关的Kalman滤波[J].测绘学报,2006,35(4):328-331.
[11]李猛,文援兰,梁加红,等.基于卡尔曼滤波的周跳探测算法[J].测绘科学技术学报,2010,27(6):407-411.