CN101079687A - 一种基于最小二乘法模型的时钟调整算法 - Google Patents

Abstract

本发明涉及一种基于最小二乘法模型的时钟调整算法，其特征在于，建立相位差的数学模型，每秒对高稳晶体的频率控制采用前N秒的鉴相结果进行预估，并将本次输出的鉴相结果保存作为下一秒预估的参数；并以2σ以外的鉴相结果对累积相位差进行补偿，随着样点增多，最后将高稳晶体的频率波动控制在一个很小的范围之内，本发明的优点是能保证时钟输出准确性。

Description

一种基于最小二乘法模型的时钟调整算法

技术领域

本发明涉及一种基于最小二乘法模型的时钟调整算法，可用于采用修正高稳晶体方法来提供高精度时钟源的时钟设备，属于通信技术领域。

背景技术

数字同步网是数字通信网正常运行的基础，也是保障各种业务网运行质量的重要手段。它与电信管理网、信令网一起并列为电信网的三大支撑网，在电信网中具有举足轻重的地位。

对于任何通信设备，都需要时钟为其提供工作频率，所以时钟性能是影响设备性能的一个重要方面。时钟常被称为设备的心脏。时钟工作时的性能主要由两个方面决定：自身性能和外同步信号的质量。而外同步信号的质量就是由数字同步网来保证的。当设备组成***和网络后，数字同步网必须为***和网络提供精确的定时，以保障其正常运行，网内各节点时钟的精度影响一个数字通信网工作是否正常。

我国的数字通信网规模庞大，分布范围广，所以数字同步网一般要接收几个基准主时钟共同控制。如果采取定时链路来传输定时信号，那么随着数字传输距离的增长，传输损伤逐渐增大、可靠性逐渐降低。而利用装配在基准钟上的GPS接收机跟踪UTC(世界协调时)，来实现对基准钟的不断调整，使之与UTC保持一致的长期频率准确度，从而达到各个基准钟同步运行和全网高度同步的目的是切实可行的，也是方便实用的。并且，在数字同步网中采用GPS配置基准钟，实现方法简单，同步时间精度高，提高了全网性能，成本却相对低廉，并且便于维护管理，所以GPS时钟在基准钟中得到广泛使用。

GPS是英文Navigation Satellite Timing and Ranging/Global PositionSystem的字头缩写词NAVSTAR/GPS的简称。它的含义是，利用导航卫星进行测时和测距，构成全球定位***，这一全球卫星定位***简称为GPS。GPS***的空间部分由24颗卫星组成，均匀分布在6个仰角为55度的轨道面上。GPS***的利用者接收卫星发送的扩频信号，测量电波传播时间求出卫星到接收机天线的距离，利用空间三球相交一点的原理，解算以接收机位置为未知数的方程，从而确切知道接收机的位置。它可以为全球用户提供连续、实时、高精度的三维位置、速度和时间信息，可以满足各种不同用户的需要。GPS***由空中卫星、地面跟踪监测站、地面卫星数据注入站、地面数据处理中心和数据通讯网络等部分组成。用户只需购买GPS接收机，就可享受免费的导航、授时和定位服务。当前，GPS时钟已成为世界上传播范围最广、精度最高的时间发布***之一，民用接收机接收到的GPS时钟与国际标准时间UCT(Universal Co-ordinated Time)保持高度同步，最高精度可达20ns。但是，民用GPS时钟的可靠性并没有得到保障，其所有者并不保证GPS的精度和可靠性，对民用用户不承担责任。而且GPS接收机接收到的GPS时钟或多或少存在以下一些误差：(1)星历误差：GPS信号中予报的卫星位置的误差，(2)卫星钟差：GPS信号中予报的卫星原子钟的偏差，(3)电离层误差：由于大气电离层效应引起的GPS信号接收的误差，(4)对流层误差：由于大气对流层效应引起的GPS信号接收的误差，(5)多径误差：由于反射信号进入接收机天线引起的GPS信号接收的误差，(6)接收机误差：由于热噪声、软件和各通道之间的硬件偏差等引起的测量值误差，(7)跟踪卫星过少误差：在某些条件下GPS接收机锁定的卫星小于4颗而产生定时误差。因此，在实际应用中，GPS接收机产生的时钟信号的精度和稳定性难以得到保证。通常GPS接收机给出的时钟精度以概率指标表示，接收机产生秒脉冲(1PPS)的误差服从正态分布。例如MOTOROLAUT ONCORE型接收机，统计精度为50ns(1σ)。表示该接收机的秒脉冲偏差服从正态分布，GPS时钟误差落于1σ范围(50ns)内的慨率为0.6828，落于2σ范围(100ns)内的慨率为0.9546；落于3σ范围(150ns)内的慨率为0.9974。但在卫星失锁或卫星时钟实验跳变的条件下，GPS时钟误差甚至达几十上百毫秒。

为了与GPS***抗衡，俄罗斯构造了GLONASS(Global Orbiting NavigationSatellite System)定位***。该***提供给民用时钟信号的精度较高，但接收机价格较贵，同样俄罗斯也不保证GLONASS的精度和可靠性，对民用用户也不承担责任，且造价高，难以得到广泛推广应用。

除了GPS时钟和GLONASS时钟，其它常规时钟频率产生方法可以是晶体、铷钟等，这一类时钟稳定性较高，单个时间间隔的随机漂移非常小，但长时间运行的累计误差较大。晶体会老化，易受外界环境变化影响，长期的精度漂移影响；原子钟长期使用后也会产生偏差，需要定时校准。而GPS***由于其工作特性的需要，定期对自身时钟***进行修正，所以其自身时钟***长期稳定，具有对外界物理因素变化不敏感特性。晶体或铷钟以GPS为长期参考，可以获得低成本、高性能的基准时钟。

为推广GPS时钟的应用，需要解决两方面问题：

(1)对GPS时钟实时监测，对秒脉冲序列精度进行评估，保证输出的准确性；

(2)对GPS时钟进行误差补偿，以提高输出时钟的精度。

为解决以上2个问题，一种方法是采用几个不同的卫星时钟***如：GPS***和GLONASS***进行相互比较校验，以提高时钟的精度和可靠性，该方法实现复杂，成本高，另一种方法是采用GPS时钟同步守时钟(原子钟、晶振时钟或CPU内部时钟)的方法，正常运行时由GPS时钟校正守时钟，在GPS失步状态下，由守时钟代替GPS时钟。当采用晶振时钟或CPU内部时钟作为守时钟时，该方法只能消除GPS短时失锁产生的较大偏差，不能消除SA干扰等产生的较小偏差。如要提高时钟精度，需采用比GPS时钟精度更高的原子钟对GPS时钟进行同步比较监测，该方法造价高，很难推广实现。

发明内容

本发明的目的是提供一种保证时钟输出准确性的基于最小二乘法模型的时钟调整算法。

为实现以上目的，本发明的技术方案是提供一种基于最小二乘法模型的时钟调整算法，其特征在于，采用C语言编制的程序，运行于时钟同步设备中，其方法为：

一种基于最小二乘法模型的GPS时钟调整算法，其特征在于，采用C语言编制的程序，运行于时钟同步设备中，其方法为：

第一步：根据压控晶体的参数计算D/A转换芯片的控制初值DA，以此初值上电控制晶体；

第二步：将恒温高稳晶体的秒脉冲与GPS接收机的秒脉冲对齐，开始跟踪；

第三步：根据相差计算对应的晶体电压调整值调整晶体，相位落后加快频率，相位超前降低频率，开始进入快速跟踪GPS的过程，此过程持续半小时，

记录下最后四百次的相位差数据，即公式9的偏差序列Y；半小时后设此时晶体控制电压D/A值为DA，记录s时间序列X即公式2值取四百；

第四步：每一秒由前X-1秒的相位差数据由线形回归算法预估本秒的相位差即公式14，根据相位和频率的对应关系以及高稳晶体频率和控制电压的计算关系计算本次对晶体的控制电压值补偿值ΔDA1即公式15，以此补偿值修正晶体的电压的控制参数为DA+ΔDA1，调整完毕后X递加；调整完毕后DA值更新为DA＝DA+ΔDA1；

第五步：如果本秒相差大于100ns(使用MOTOROLA VP ONCORE时为2σ单位，此时秒鉴相结果(超前/滞后)准确性的概率等于GPS秒脉冲落于2σ范围内的慨率为95.46％)时可认为相位差过大，必须进行一定相位补偿，补偿的相位为超出100ns的部分，根据相位和频率的对应关系以及高稳晶体频率和控制电压的计算关系计算本次对晶体的控制电压值补偿值ΔDA2(见公式15)，以此补偿值修正晶体电压的控制参数为DA+ΔDA2；调整完毕后DA值不更新；

第六步：重复第四步反复修正，时间越长，样点数越多，此时压控晶体的每秒钟的电压控制参数DA将沿着某个中心值上下小幅波动，总体趋近于一条水平直线，通过线性回归克服了GPS秒脉冲抖动带来的影响。

本发明采用高精度晶振对GPS时钟进行监测与校正的简便实用方法，建立GPS时钟误差的测量模型，提出一种高精度时钟的产生方法，并成功地将其应用于时钟同步设备中。

GPS接收机在正常工作条件下，其时钟的误差服从正态分布，只存在单个秒脉冲的抖动，但从长时间来看，GPS时钟并不存在累计误差。而高稳晶体短时间内的抖动很小，但存在较大的累计误差。GPS时钟与晶振时钟的精度是互补的，如果把二者进行比较分析，使二者互为参考，利用GPS的长期精度特性对高稳晶体进行校正，采用数理统计的方法估计出二者的误差，进而对误差进行在线主动补偿，可以实现高精度时钟。

在具体实现时采用：对高精度晶振进行分频，产生秒脉冲时钟信号；晶振秒脉冲时钟与GPS的秒脉冲进行相位比较，产生偏差序列，该偏差包括GPS时钟的左右抖动相差和相位的累计误差；采用数学回归对两种误差进行估计，从而分离出各自的误差，并对晶体累计误差进行修正，由此构造一种简便的高精度时钟发生装置。具体为通过分频电路，将高精度时钟的输出时钟进行分频，分频后得到秒脉冲信号经过数字鉴相滤器和GPS的秒脉冲信号进行相位比较，得到相位偏差序列样本，晶振的累积误差由CPU设置补偿值修正。每秒调整一次，第N秒时的补偿值由前N-1秒的GPS时钟误差和前N-1秒的补偿值的历史数据采用线性回归的算法得到预估的补偿值。

本发明的优点是能保证GPS时钟输出准确性。

附图说明

图1为时钟同步设备的结构示意图；

图2为数字鉴相滤波器的结构示意框图；

图3为数字鉴相滤波器内部的顶层电原理图；

图4为“DP_CTRL”模块中的电原理图；

图5为最小二乘法模型的时钟调整算法流程框图。

具体实施方式

以下结合附图和实施例对本发明作进一步说明。

实施例

如图1所示，为时钟同步设备为结构示意图，所述的时钟同步设备由数字鉴相滤器、GPS接收机、CPU和恒温高稳晶体组成，数字鉴相滤器与GPS接收机连接，同时数字鉴相滤器与CPU和恒温高稳晶体连接组成回路，CPU：Rabbit2000，性能约等同于采用16M晶振的8051系列。

如图2所示，为数字鉴相滤波器的结构示意框图，本发明针对这种低成本的GPS接收模块所提供的1pps的秒脉冲信号(或等效1pps信号)，作为外参考比较信号源，利用数模混和的锁相环路所控制的本地的压控石英晶体振荡器的输出等效1pps信号作为本地比较信号源，通过用FPGA进行数字鉴相及滤波技术，可灵活调整输出信号的特性及精度，以满足各种需要场合。本发明中，使用了两个本地高频信号源，1)，由16.384MHZ的压控晶体振荡器组成的可控本地标准高频时钟源，通过分频，可产生2.048MHZ的输出时钟及本地1pps信号。2)，由16.384MHZ的4倍频既65.536MHZ信号，作为数字鉴相及滤波工作时钟。数字鉴相器比较本地秒脉冲信号与外参考秒脉冲信号的相位差，再通过数字滤波，及外部配套软件的算法分析，可解决GPS模块的1pps输出扰动问题，以保证本地输出信号精度，由此可产生出与输入信号源精度相当的多种信号输出。

本发明实施例所提供的一种数字鉴相滤波器，以及分别与数字鉴相滤波器连接的GPS模块、2M频率信号源模块、CPU和16.384MHZ压控频率单元；

由于数字鉴相及数字滤波还有一些辅助单元均是在一块FPGA(可编程逻辑器件)中实现，所以在这里就FPGA内部的设计实现做详细的说明。

如图3所示，数字鉴相滤波器内部的顶层电原理如下：

顶层包括3个模块(模块1～模块3)，每个模块的左边引脚均为输入，右边均为输出。但是标注“双向数据总线”的引脚为输入输出双向。

1)接口模块1：“IO_INTERFACE”，负责与外部CPU进行数据交换，接收配置指令，上报本地秒信号与外参考秒信号的相位相对关系等；接口模块1的输入端分别连接GPS模块和2M频率信号源模块，经数据总线连接CPU；根据配置，选出一个外参考源，并将其转为秒脉冲输出EXT_1PPS，作为鉴相的外参考依据。例如，CPU通过地址线ADDRES[31..24]，/CS，ALE，/WR配合，将地址为x“00”所对应的数据线DB[7..0]上的内容写入本模块，若数据内容为x“00”，则表明选“GPS_1PPS”为外参考源，接口模块可直接将其输出到模块输出端口“SEL_1PPS”，供外部模块使用；若数据内容为x“01”，则表明选“2MHZ[4..1]”中的第一路(2MHZ[1])为外参考源，由于此时是2.048MHZ的频率输入，故模块对其做分频，产生1pps的秒脉冲，并将其输出到模块输出端口“SEL_1PPS”，供外部模块使用；其他类推，这样“外参考”秒脉冲就产生了。

同理，CPU也将输出信号的占空比参数，滤波器的参数写入模块，分别输出到接口模块“PW[7..0]”和“FILT_D[7..0]”端口上，供外部模块使用。

另外，来自外模块的“UP”、“DOWN”、“OVER_DIR”相位比较输出信号也通过同样的方式，由本接口模块提供给CPU访问读取。

这样，CPU就可不间断的通过数据线掌握外参考源与本地源之间的相位相对关系，从而对压控晶体振荡器组成的可控本地标准高频时种源进行调整，这样持续反复的进行，通过CPU的特定算法，就最终达到闭环反馈环路的稳定状态，本地输出信号也就精密的保持对外参考源的跟踪。

2)本地模块(LOCAL_CQ)2，用于完成对本地晶体振荡器输入频率的分频，按16.384MHZ转到1hz的分频系数进行分频，并抽出其他所需频率，同时还要按CPU配置的脉冲宽度数据调整输出秒脉冲的占空比；本地模块2的输入端分别连接接口模块1的外参考秒脉冲输出端、滤波器参数输出端、重置位输出端和16.384MHZ频率信号源模块，本地模块2分别设有其他所需频率如秒脉冲、2.048MHZ、200HZ、8KMHZ等信号的输出端；另外，当接收到CPU的重行同步置位指令后，要将分频器状态与外参考秒信号的相位重新对齐一次。所产生的秒脉冲L1PPS是鉴相的本地参考依据。

这里的核心就是一个可异步置位的同步时钟定模计数器，模值为16.384×10⁶，每当计数值归零时，既为本地秒脉冲的起始边沿，而高电平的宽度，即由“PW_D[7..0]”端口数值决定，当计数器值小于等于PW_D[7..0]的值时，输出维持高电平，反之为低，若要扩充调整范围，可通过增加脉冲配置数据“PW_D[7..0]”的位数来实现，由此，“本地”秒脉冲就产生了。另外，当RESYN端收到重置位信号时，就产生一个窄脉冲，并且保证使之与外参考的秒脉冲上升沿同步，由这个窄脉冲对计数器进行异步清零，这样，重置位之后，本地计数器的状态就与外参考源保持一致，这样处理的目的，是为了加快对外参考源的跟踪锁定进程，而且还可以使后继模块的数字鉴相及滤波电路得以简化，因为重置位之后，本地秒脉冲的相位已与外参考源的秒脉冲相位近乎一致，在短时间段内，两者之间的相差漂移不会太大，所以滤波电路的观察范围就可以缩小。

3)数据处理模块(DP_CTRL)3，用于内外两个脉冲之间的相位检测，还包括可调整数字滤波器。数据处理模块3的输入端分别连接接口模块1的输出端、本地模块2的输出端和16.384MHZ频率信号源模块的4倍频输出端，数据处理模块3的输出端连接接口模块1的输入端；

本数据处理模块的相位检测，就是检测两个脉冲之间的相位相对关系，这两个输入脉冲就是先前描述过的外部模块已处理过的“外参考”秒脉冲，“本地”秒脉冲，鉴相结果可反映出“外参考”秒脉冲与“本地”秒脉冲之间的相对位置和变化趋势，鉴相结果是每秒种更新一次。

可调整数字滤波器的引入，一方面是可以对偶发的电路干扰做屏蔽，最主要的作用在于，当“外参考”秒脉冲是GPS接收模块所产生时，其固有的秒脉冲“扰动”相对较大，以高精度的时钟***来衡量，其相位的跳变是不符合要求的，同时还会对鉴相器的鉴相结果造成干扰，尤其是当本地的时钟源事实上已与外参考源保持高度一致时，这一现象会更明显，为此，滤波器在这里的作用就是尽量将这种相位“跳变”过滤掉，或者说，就是将“突发”的“大相位差”上报给CPU，CPU会根据一定的算法，判断出当前的实际状态，从而避免不必要的调整或者“过调整”而导致的输出信号性能指标下降。但是，一个实际使用上的问题需要考虑，因为个体电路元器件之间的潜在差异，尤其是可能的GPS接收模块的电性能的差异(例如：生产制造商的不同，生产的批次不同等等)，都有可能让使用固定参数的的滤波器失去效能，为此，这里的滤波器参数是可设的，并且实时的由CPU来控制，由此，通过“过滤”的相位差宽度就可以改变，大大提高了电路效能的可靠性和实用性。

另外，由于“速度”和“精度”向来是一对矛盾，滤波范围变小，***跟踪精度变高，但***进入稳定态就需时更长，反之，滤波范围变大，***跟踪精度变低，但***进入稳定态就更快，为此，针对不同的应用场合，适当调整滤波器的参数设置，可改变***的整体响应速度，提高工作效率。

如图4所示，为“DP_CTRL”模块中的电原理图，所述数据处理模块内设有鉴相器单元4和滤波器单元5、6；

鉴相器单元(V_DPD)4：为用VHDL硬件描述语言构建的数字鉴相器，鉴相器单元4的输入端分别连接接口模块1的外参考秒脉冲输出端和本地模块2的本地参考秒脉冲输出端，用于识别出两个输入脉冲的相位关系，利用脉冲的边沿来鉴别，和输入的脉冲宽度无关，有很高的分辨率，当输入的本地秒信号频率高于或相位超前外参考秒信号时，输出“DOWNDIR”就变‘1’，反之“UPDIR”就变‘1’。均为‘0’时，表示两个输入脉冲完全一致，即同频同相。“DOWNDIR”为‘1’，指示要将本地压控晶振的频率降下来，而“UPDIR”为‘1’时，表示要将本地压控晶振的频率上调加快。鉴相器单元4的加频指示、减频指示输出端分别经与门连接接口模块1的加频控制、减频控制输入端。

a)超前滤波器单元(V_OVER_CTRL)5、滞后滤波器单元(V_OVER_CTRL)6，为用VHDL硬件描述语言构建的数字滤波器，用于溢出指示；当两个输入脉冲相差超出FID[7..0]所设定的值时，“OVER_ACT”输出‘1’电平，这里使用了两个这样的模块，就是为了区分两种相对状态：a。本地秒信号超前外参考秒信号溢出，b。本地秒信号滞后外参考秒信号溢出，并合二为一成“OVER_DIR”输出，供CPU检测，FID[7..0]是由CPU所设，可动态调整，它反映的是“clk×4”输入时钟(既65.536MHZ)的周期数，改变其值，可改变输出时钟的特性参数。

超前滤波器单元5的输入端分别连接鉴相器单元4的加频指示输出端、16.384MHZ频率信号源模块的4倍频输出端、接口模块1的滤波器参数输出端，滞后滤波器单元6的输入端分别连接鉴相器单元4的减频指示输出端、16.384MHZ频率信号源模块的4倍频输出端、接口模块1的滤波器参数输出端，超前滤波器单元5和滞后滤波器单元6的溢出指示输出端经或门连接接口模块1的溢出控制输入端。

采用数学回归对两种误差进行估计，从而分离出各自的误差，并对晶体累计误差进行修正，由此构造一种简便的高精度时钟发生装置。具体为通过分频电路，将高精度时钟的输出时钟进行分频，分频后得到秒脉冲信号经过数字鉴相滤器和GPS的秒脉冲信号进行相位比较，得到相位偏差序列样本，晶振的累积误差由CPU设置补偿值修正。每秒调整一次，第N秒时的补偿值由前N-1秒的GPS时钟误差和前N-1秒的补偿值的历史数据采用线性回归的算法得到预估的补偿值。算法如下：

第一步.建立GPS时钟同步晶振时钟的数学模型

GPS输出时钟与国际标准时间存在一定误差e，e服从正态分布：

ε～N(0，σ²)。    (1)

s时间序列X：

1，2，3，4，…，x，…，n    (2)

GPS输出的s时钟序列对应的国际标准时间可记为

1-e₁，2-e₂，3-e₃，4-e₄，…，x-e_x，…，n-e_n    (3)

通用公式表示为y_x′＝x-ε_x x∈N ε_x～(0，σ²)  (4)

式中y_x′为GPS输出第x个s时钟对应的国际标准时间，时间误差为e_x。

设晶振分频产生的s时钟序列第0个s时钟与国际标准时间的初始偏差为a；时间间隔误差为b；由于高精度晶振的随机误差远小于GPS的s时钟的随机误差(如精度为10^-9s的晶振随机误差＜1ns)，因此不考虑晶振s时钟的随机误差，晶振分频输出s时钟产生的时钟序列对应的国际标准时间可记为

1+a+b，2+a+2b，3+a+3b，4+a+4b，…，x+a+bx，…，n+a+bn    (5)

通用公式表示为y_x″＝x+a+bx  x∈N   (6)

式中y_x″为晶振分频输出第x个s时钟对应的国际标准时间，时间误差为

m(x)＝a+bc        (7)

则晶振分频s时钟(简称晶振s时钟)与GPS的s时钟的偏差为

y_x＝y_x″-y_x′＝a+bx+ε_x  x∈N      (8)

偏差序列Y可表示为

y₁，y₂，y₃，y₄，y₅，y₆，…，y_x，…，y_n    (9)

由Y的时间序列可以对a，b进行估计，从而计算晶振s时钟误差m(x)；同时由Y的时间序列可以对GPS的误差e_x进行估计，从而可以衡量GPS接收机的精度水平。

第二步.晶振s时钟的误差估计

采用式(8)的一元线性回归模型，即采用式(7)的回归方程，分析时间序列X与偏差序列Y的相关特性，对a，b进行最小二乘估计，估计值分别为

$\hat{b}=\frac{\underset{x=1}{\overset{n}{\mathrm{\Σ}}}(x-\stackrel{\‾}{x})(y_x-\stackrel{\‾}{y})}{\underset{x=1}{\overset{n}{\mathrm{\Σ}}}{(x-\stackrel{\‾}{x})}^2}---\left(10\right)$

$\hat{a}=\stackrel{\‾}{y}-\hat{b}\stackrel{\‾}{x}---\left(11\right)$

式中 y， x分别为偏差序列Y及时间序列X的平均值，

$\stackrel{\‾}{x}=\frac{1}{n}\underset{x=1}{\overset{n}{\mathrm{\Σ}}}x=\frac{n+1}{2}---\left(12\right)$

$\stackrel{\‾}{y}=\frac{1}{n}\underset{x=1}{\overset{n}{\mathrm{\Σ}}}{y}_x---\left(13\right)$

则晶振输出的第x个s时钟与国际标准时间之间的误差估计值

为

$\widehat{\mathrm{\μ}}\left(x\right)=\hat{a}+\hat{b}x---\left(14\right)$

得到的

是一个时间量，单位为秒，还需将其和高稳晶体的控制量电压(单位为V)进行转换。控制电压由D/A芯片的控制量决定，转换公式如下：

$\mathrm{\ΔDA}=\frac{\widehat{\mathrm{\μ}}\left(x\right)*D_{\mathrm{MAX}}}{\frac{f_{\mathrm{Max}}}{f_{\mathrm{MIN}}}-1}---\left(15\right)$

式中D_MAX取决于DA转换芯片的位数，如DA为十六位转换芯片，则D_MAX＝65536。如DA为八位转换芯片，则D_MAX＝256。f_MIN为压控高稳晶体调节到的最小频率，f_MAX为压控高稳晶体可调节到的最大频率。

根据误差估计值

得到的ΔDA对晶振s时钟进行调整，即可产生高精度时钟输出。

第三步.GPS接收机输出的s时钟误差估计

GPS接收机输出的秒脉冲误差的估计，实际上是对式(1)中σ²的估计，由于σ²＝D(ε)＝E(ε²)是e的二阶原点距，按距估计法，可用样本二阶原点距作为它的估计值

${\widehat{\mathrm{\σ}}}^2=\frac{1}{n}\underset{x=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ϵ}}_x^2---\left(16\right)$

而由式(8)得

ε_x＝y_x-a-bx          (17)

用 分别替代a，b，即得σ²的估计量

${\widehat{\mathrm{\σ}}}^2=\frac{1}{n}\underset{x=1}{\overset{n}{\mathrm{\Σ}}}{(y_x-\hat{a}-\hat{b}x)}^2---\left(18\right)$

的大小反映接收机输出秒脉冲的精度。当误差估计值

时，认为GPS接收机失步或故障。

如图5所示，为最小二乘法模型的GPS时钟调整算法流程框图，本发明通过分频电路，将高精度时钟的输出时钟进行分频，分频后得到秒脉冲信号经过数字鉴相滤器和GPS的秒脉冲信号进行相位比较，得到相位偏差序列样本，晶振的累积误差由CPU设置补偿值修正。每秒调整一次，第N秒时的补偿值由前N-1秒的GPS时钟误差和前N-1秒的补偿值的历史数据采用线性回归的算法得到预估的补偿值。以此补偿值作为修正对高稳晶体频率进行调整，并由本次的鉴相结果计算出下一秒的补偿值，其步骤为：

第一步：根据压控晶体的参数计算D/A转换芯片的控制初值DA，以此初值上电控制晶体。具体方法如下：

本例高稳压控晶体的频率调节范围为16.384M±15Hz，DA转换芯片为16位，参考电压5V，DA赋值为0时输出电压为0，输出频率16.384M-15Hz。DA赋值为0时输出电压为5V，输出频率16.384M+15Hz。为了保证上电频率初始值大致为16.384M，DA取值32768即可；

第二步：将晶体的秒脉冲与GPS的秒脉冲对齐，开始跟踪；

第三步：根据相差计算对应的晶体电压调整值调整晶体，相位落后加快频率，相位超前降低频率，开始进入快速跟踪GPS的过程。此过程的目的是为了使晶体的输出频率的中心值能进入一个比较准确的范围，并在此范围内小幅波动。此过程持续半小时，记录下最后四百次的相位差数据(即公式9的偏差序列Y)。此时晶体控制电压D/A值为DA，本例DA为32210。记录s时间序列X(公式2)值取四百；

第四步：每一秒由前X-1秒的相位差数据由线形回归算法预估本秒的相位差

(设本例为+5ns)(见公式14)，根据相位和频率的对应关系以及高稳晶体频率和控制电压的计算关系计算本次对晶体的控制电压值补偿值ΔDA1＝-179(见公式15)。以此补偿值修正晶体的电压的控制参数为为DA+ΔDA1＝32210-179＝32031，调整完毕后X递加为401；调整完毕后DA值更新为DA＝DA+ΔDA1＝32031；

第五步：如果本秒相差大于100ns(使用MOTOROLA VP ONCORE时为2σ单位，此时秒鉴相结果(超前/滞后)准确性的概率等于GPS秒脉冲落于2σ范围内的慨率为95.46％)如本例为+105ns时可认为相位差过大，必须进行一定相位补偿，补偿的相位为超出100ns的部分，根据相位和频率的对应关系以及高稳晶体频率和控制电压的计算关系计算本次对晶体的控制电压值补偿值ΔDA2＝-179(见公式15)，以此补偿值修正晶体电压的控制参数为DA+ΔDA2＝32031-179＝31852；调整完毕后DA值不更新；

第六步：重复第四步反复修正，时间越长，样点数越多，此时压控晶体的每秒钟的电压控制参数DA将沿着某个中心值上下小幅波动，总体趋近于一条水平直线。通过线性回归克服了秒脉冲抖动带来的影响。

Claims (9)

1.一种基于最小二乘法模型的时钟调整算法，其特征在于，采用C语言编制的程序，运行于时钟同步设备中，其方法为：

第一步：根据压控晶体的参数计算D/A转换芯片的控制初值DA，以此初值上电控制晶体；

第二步：将恒温高稳晶体的秒脉冲与GPS接收机的秒脉冲对齐，开始跟踪；

第三步：根据相差计算对应的晶体电压调整值调整晶体，相位落后加快频率，相位超前降低频率，开始进入快速跟踪GPS的过程，此过程持续半小时，记录下最后四百次的相位差数据，即公式9的偏差序列Y；半小时后设此时晶体控制电压D/A值为DA，记录s时间序列X即公式2值取四百；

第四步：每一秒由前X-1秒的相位差数据由线形回归算法预估本秒的相位差

即公式14，根据相位和频率的对应关系以及高稳晶体频率和控制电压的计算关系计算本次对晶体的控制电压值补偿值ΔDA1即公式15，以此补偿值修正晶体的电压的控制参数为DA+ΔDA1，调整完毕后X递加；调整完毕后DA值更新为DA＝DA+ΔDA1；

第五步：如果本秒相差大于100ns时可认为相位差过大，必须进行一定相位补偿，补偿的相位为超出100ns的部分，根据相位和频率的对应关系以及高稳晶体频率和控制电压的计算关系计算本次对晶体的控制电压值补偿值ΔDA2即公式15，以此补偿值修正晶体电压的控制参数为DA+ΔDA2；调整完毕后DA值不更新；

第六步：重复第四步反复修正，时间越长，样点数越多，此时压控晶体的每秒钟的电压控制参数DA将沿着某个中心值上下小幅波动，总体趋近于一条水平直线，通过线性回归克服了GPS秒脉冲抖动带来的影响。

2.根据权利要求1所述的一种基于最小二乘法模型的时钟调整算法，其特征在于，所述的步骤(3)-步骤(6)中的计算方法为：

第一步：建立GPS时钟同步晶振时钟的数学模型

GPS输出时钟与国际标准时间存在一定误差e，e服从正态分布：

ε～N(0，σ²)；    (1)

s时间序列X：

1，2，3，4，...，x，...，n    (2)

GPS输出的s时钟序列  对应的国际标准时间可记为

1-e₁，2-e₂，3-e₃，4-e₄，...，x-e_x，...，n-e_n      (3)

通用公式表示为y_x′＝x-ε_x   x∈N   ε_x～(0，σ²)  (4)

式中y_x′为GPS输出第x个s时钟对应的国际标准时间，时间误差为e_x，

设晶振分频产生的s时钟序列第0个s时钟与国际标准时间的初始偏差为a；时间间隔误差为b；不考虑晶振s时钟的随机误差，晶振分频输出s时钟产生的时钟序列对应的国际标准时间可记为

1+a+b，2+a+2b，3+a+3b，4+a+4b，...，x+a+bx，...，n+a+bn    (5)

通用公式表示为y_x″＝x+a+bx   x∈N         (6)

式中y_x″为晶振分频输出第x个s时钟对应的国际标准时间，时间误差为

m(x)＝a+bx                                (7)

则晶振分频s时钟(简称晶振s时钟)与GPS的s时钟的偏差为

y_x＝y_x″-y_x′＝a+bx+ε_x     x∈N          (8)

偏差序列Y可表示为

y₁，y₂，y₃，y₄，y₅，y₆，...，y_x，...，y_n  (9)

由Y的时间序列可以对a，b进行估计，从而计算晶振s时钟误差m(x)；同时由Y的时间序列可以对GPS的误差e_x进行估计，从而可以衡量GPS接收机的精度水平；

第二步：晶振s时钟的误差估计

采用式(8)的一元线性回归模型，即采用式(7)的回归方程，分析时间序列X与偏差序列Y的相关特性，对a，b进行最小二乘估计，估计值分别为

$\hat{b}=\frac{\underset{x=1}{\overset{n}{\mathrm{\Σ}}}(x-\stackrel{\‾}{x})(y_x-\stackrel{\‾}{y})}{\underset{x=1}{\overset{n}{\mathrm{\Σ}}}{(x-\stackrel{\‾}{x})}^2}---\left(10\right)$

$\hat{a}=\stackrel{\‾}{y}-\hat{b}\stackrel{\‾}{x}---\left(11\right)$

式中 y， x分别为偏差序列Y及时间序列X的平均值， $\stackrel{\‾}{x}=\frac{1}{n}\underset{x=1}{\overset{n}{\mathrm{\Σ}}}x=\frac{n+1}{2}---\left(12\right)$

$\stackrel{\‾}{y}=\frac{1}{n}\stackrel{n}{\underset{x=1}{\mathrm{\Σ}}}{y}_x---\left(13\right)$

则晶振输出的第x个s时钟与国际标准时间之间的误差估计值 为

$\widehat{\mathrm{\μ}}\left(x\right)=\hat{a}+\hat{b}x---\left(14\right)$

得到的

是一个时间量，单位为秒，还需将其和高稳晶体的控制量电压(单位为V)进行转换，控制电压由D/A芯片的控制量决定，转换公式如下：

$\mathrm{\ΔDA}=\frac{\widehat{\mathrm{\μ}}\left(x\right)*D_{\mathrm{MAX}}}{\frac{f_{\mathrm{Max}}}{f_{\mathrm{MIN}}}-1}---\left(15\right)$

根据误差估计值

得到的ΔDA对晶振s时钟进行调整，即可产生高精度时钟输出；

第三步：GPS接收机输出的s时钟误差估计

GPS接收机输出的秒脉冲误差的估计，实际上是对式(1)中σ²的估计，由于σ²＝D(ε)＝E(ε²)是e的二阶原点距，按距估计法，可用样本二阶原点距作为它的估计值

${\widehat{\mathrm{\σ}}}^2=\frac{1}{n}\underset{x=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ϵ}}_x^2---\left(16\right)$

而由式(8)得

ε_x＝y_x-a-bx         (17)

用

分别替代a，b，即得σ²的估计量

${\widehat{\mathrm{\σ}}}^2=\frac{1}{n}\underset{x=1}{\overset{n}{\mathrm{\Σ}}}{(y_x-\hat{a}-\hat{b}x)}^2---\left(18\right)$

的大小反映接收机输出秒脉冲的精度，当误差估计值

时，认为GPS接收机失步或故障。

3.根据权利要求1所述的一种基于最小二乘法模型的时钟调整算法，其特征在于，所述的时钟同步设备由数字鉴相滤波器、GPS接收机、CPU和恒温高稳晶体组成，数字鉴相滤器与GPS接收机连接，同时数字鉴相滤波器与CPU和恒温高稳晶体连接组成回路。

4.根据权利要求3所述的一种基于最小二乘法模型的时钟调整算法，其特征在于，所述的数字鉴相滤波器在FPGA(可编程逻辑器件)中实现，分别与数字鉴相滤波器连接的有GPS模块、标准信号源模块、CPU和可控本地标准压控频率单元，构成一高精度的时钟装置，其特征在于，所述数字鉴相滤波器包括：

一接口模块，用于与包括CPU在内的FPGA的外部模块进行信息交换，即将内外两个脉冲之间的相位检测信息输出给CPU，将CPU的占空比参数、滤波器参数、重置位等输入信号传输给FPGA内相应模块；其输入端分别连接所述GPS模块和所述标准信号源模块，经数据总线连接所述CPU；

一本地模块，用于完成对本地晶体振荡器输入频率的分频，及本地秒脉冲的重置位，并输出其他所需频率；其输入端分别连接所述接口模块的外参考秒脉冲输出端、滤波器参数输出端、重置位输出端和可控本地标准压控频率单元模块，所述本地模块分别设有其他所需频率的输出端；

一数据处理模块，包括调整由CPU设置滤波器参数的可调整数字滤波器，用于内外两个脉冲之间的相位检测；其输入端分别连接所述接口模块的输出端、所述本地模块的输出端和可控本地标准压控频率单元模块的倍频输出端，其输出端连接所述接口模块的输入端；

所述CPU依据内外两个脉冲之间的相位检测信息，控制本地模块的重置位、滤波器的参数、输出信号的占空比。

5.根据权利要求4所述的一种基于最小二乘法模型的时钟调整算法，其特征在于，所述的数据处理模块内设有：

一鉴相器单元，用于识别出两个输入脉冲的相位关系，其输入端分别连接接口模块的外参考秒脉冲输出端和本地模块的本地参考秒脉冲输出端，其加频指示、减频指示输出端分别经与门连接接口模块的加频控制、减频控制输入端；

一超前滤波器单元和一滞后滤波器单元，用于溢出指示；超前滤波器单元的输入端分别连接鉴相器单元的加频指示输出端、可控本地标准压控频率单元模块的倍频输出端、接口模块的滤波器参数输出端；滞后滤波器单元的输入端分别连接鉴相器单元的减频指示输出端、可控本地标准压控频率单元模块的倍频输出端、接口模块的滤波器参数输出端，超前滤波器单元和滞后滤波器单元的溢出指示输出端经或门连接接口模块的溢出控制输入端。

6.根据权利要求4所述的一种基于最小二乘法模型的时钟调整算法，其特征在于，所述的标准信号源模块为2M频率信号源模块。

7.根据权利要求4所述的一种基于最小二乘法模型的时钟调整算法，其特征在于，所述的可控本地标准压控频率单元为16.384MHZ的压控频率单元。

8.根据权利要求4所述的一种基于最小二乘法模型的时钟调整算法，其特征在于，所述的本地模块输出端的其他所需频率为秒脉冲、2.048MHZ、200HZ、8KMHZ。

9.根据权利要求4所述的数字鉴相滤波器，其特征是，所述可控本地标准压控频率单元模块的倍频输出端为16.384MHZ压控频率单元模块的4倍频输出端。

Images