CN113009523B - 长时间相干积分捕获的多普勒频率估计与补偿方法及系统 - Google Patents
长时间相干积分捕获的多普勒频率估计与补偿方法及系统 Download PDFInfo
- Publication number
- CN113009523B CN113009523B CN202110198507.2A CN202110198507A CN113009523B CN 113009523 B CN113009523 B CN 113009523B CN 202110198507 A CN202110198507 A CN 202110198507A CN 113009523 B CN113009523 B CN 113009523B
- Authority
- CN
- China
- Prior art keywords
- doppler frequency
- signal
- module
- frequency
- gnss
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 70
- 230000001427 coherent effect Effects 0.000 title claims abstract description 35
- 230000010354 integration Effects 0.000 title claims abstract description 33
- 238000005070 sampling Methods 0.000 claims abstract description 32
- 239000013598 vector Substances 0.000 claims abstract description 31
- 239000011159 matrix material Substances 0.000 claims abstract description 25
- 238000004422 calculation algorithm Methods 0.000 claims description 58
- 238000001228 spectrum Methods 0.000 claims description 21
- 230000008859 change Effects 0.000 claims description 17
- 238000004364 calculation method Methods 0.000 claims description 16
- 238000001514 detection method Methods 0.000 claims description 16
- 238000003860 storage Methods 0.000 claims description 13
- 230000009466 transformation Effects 0.000 claims description 12
- 230000007774 longterm Effects 0.000 claims description 5
- 230000008569 process Effects 0.000 claims description 5
- 238000005316 response function Methods 0.000 claims description 5
- 238000004519 manufacturing process Methods 0.000 claims description 4
- 238000011897 real-time detection Methods 0.000 claims description 4
- 230000003595 spectral effect Effects 0.000 claims description 3
- 241000764238 Isis Species 0.000 claims 1
- 230000035945 sensitivity Effects 0.000 abstract description 12
- 230000007613 environmental effect Effects 0.000 abstract description 3
- 238000009825 accumulation Methods 0.000 description 7
- 230000000737 periodic effect Effects 0.000 description 6
- 238000011160 research Methods 0.000 description 5
- 238000012545 processing Methods 0.000 description 4
- 230000009286 beneficial effect Effects 0.000 description 3
- 239000000969 carrier Substances 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 230000000903 blocking effect Effects 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- 230000006835 compression Effects 0.000 description 1
- 238000007906 compression Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 230000002035 prolonged effect Effects 0.000 description 1
- 230000007480 spreading Effects 0.000 description 1
- 238000003892 spreading Methods 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 239000005436 troposphere Substances 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
- G01S19/13—Receivers
- G01S19/24—Acquisition or tracking or demodulation of signals transmitted by the system
- G01S19/25—Acquisition or tracking or demodulation of signals transmitted by the system involving aiding data received from a cooperating element, e.g. assisted GPS
- G01S19/254—Acquisition or tracking or demodulation of signals transmitted by the system involving aiding data received from a cooperating element, e.g. assisted GPS relating to Doppler shift of satellite signals
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Computer Networks & Wireless Communication (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Position Fixing By Use Of Radio Waves (AREA)
Abstract
本发明公开了一种长时间相干积分捕获的多普勒频率估计与补偿方法及系统。该方法中，首先将GNSS零中频数据IGIFS作为输入信号分别与m条相位补偿序列相乘后，将每个相乘结果划分为M个含N个采样点数据块，将M个数据块对应位相加后得到一N维向量；m个N维向量视为m×N维矩阵；再分别遍历矩阵中的每一列j，将m维列向量中的最大值以及对应的位序号n分别保存于一个2×N维矩阵的第j列的两个元素中；最后搜索并记录该2×N维矩阵第一行N个元素的最大值，并确定该最大值对应的序号j*和对应第二行保存的相位补偿序列的序号n*；计算载波多普勒频率，本发明可以在信号受到遮挡或一定的环境干扰下保持正常接收GNSS信号，提高GNSS信号多普勒频率估计的灵敏度和估计精度。
Description
技术领域
本发明属于卫星信号处理领域，具体涉及一种长时间相干积分捕获的多普勒频率估计与补偿方法及系统。
背景技术
卫星导航定位系统(GNSS)是一种以卫星为基础的无线电导航系统，能为陆、海、空的各类载体提供全天候、不间断、高精度、实时导航定位服务[1，2]，已经应用于国民经济与日常生活的各个领域，如地面交通监管、飞机与船舶导航、精密受时、大地测量等[2]。目前，全球范围内研发最早、应用最早的卫星定位系统GPS系统在我国已得到广泛应用，我国自主研发全球卫星定位系统Compass(北斗二代)，该系统2012年底已经在我国及其周边地区提供定位服务，北斗3号系统已于2020年6月建成并提供全球导航定位服务。因此，研究GNSS系统的高性能接收机技术必将成为国内未来一段时间内的研究重点。
GPS信号(扩频信号)到达地面接收机时已相当微弱，大约为-130dBmW，比接收机内部热噪声低20～30dB。特别的，在室内、城市、森林等复杂环境(本文中统称为室内环境)中，GPS接收信噪比更低，而室内环境恰是人类活动的主要环境之一。
文献调研表明，长时间预相干积分是进一步提高捕获处理增益的首选方法。
当预相干积分累加时间超过20ms时，载波频率搜索步长不能大于50Hz，对于卫星信号中约±10kHz的多普勒频率变化范围来说，搜索时间耗费巨大。为此，本专利将估计并补偿GNSS信号中的多普勒频率，压缩频率搜索范围，降低长时间预相干积分捕获计算的复杂度。
GNSS信号高灵敏度捕获算法的研究是很有必要的。目前，虽然国内外已经有相关文献研究GNSS信号的高灵敏度捕获算法，在高灵敏度捕获、快速捕获算法建模等方面已做了一定研究工作，如有文献提出了二倍分组块补零(DBZP)算法、组合比特翻转估计算法、多级相干累加捕获算法，将累加时间充分的延长，但最大可能预相干积分时间仍没有大于1个导航电文比特长度，平方损耗和多普勒残差仍然是低信噪比环境下影响捕获性能的主要因素，而且，这些算法很难同时兼顾捕获灵敏度和捕获效率这两个方面。为了实现快速高效的GNSS信号高灵敏度捕获，在民航、铁路交通等对可靠性要求较高领域中得到推广应用，亟需研究辅助GNSS信号高灵敏度快速捕获的多普勒频率估计与补偿算法与理论，解决上述关键科学问题。
发明内容
为了克服现有GPS接收机灵敏度低和动态性能低而不能给出定位结果的不足，本发明提供了一种辅助实现GNSS信号高灵敏度快速捕获的多普勒频率估计与补偿算法方法及系统，使用该模块的接收机可提高GNSS信号捕获速度及其捕获灵敏度，提高GNSS弱信号的检测和锁定能力，并且，当信号受到遮挡、环境噪声较大时也可稳定的捕获卫星信号。
本发明所采用的具体技术方案如下：
一种长时间相干积分捕获的多普勒频率估计与补偿方法，其步骤如下：
S1：获取长度为TIms的GNSS零中频数据IGIFS作为输入信号r(·)；
S2：选取一组调整频率δf，得到一组共m条相位补偿序列Ts为采样周期，上标j为虚数单位；对于其中任意第n条相位补偿序列β(n)，n＝1，2，…m，分别执行S21～S23的计算过程：
S21：将所述输入信号r(·)与相位补偿序列值β(n)相乘，得到乘积结果然后根据伪码周期将/>分为M块，每块N个采样数据，M块数据块叠加表示为：
其中，表示/>中第i+1块数据块，r(n)表示/>中第1块数据块，v(n)为噪声，fd表示多普勒频率；HLN(fd)为分块叠加操作的响应函数，计算公式为：
S22：用平方运算去除结果rLN(n)中调制的伪随机码PRN码，得到分块信号数据叠加结果
S23：通过快速傅立叶变换将分块信号数据叠加结果转换到频域上，得到其中f表示自变量为频率；
S3：分别针对m条相位补偿序列执行S21～S23后，得到m个N维向量分别记为/>个N维向量视为一个m×N维矩阵；然后，分别遍历m个向量/>中的每一维，将m个向量/>中第j维即m×N维矩阵第j列的最大值以及该最大值所属向量/>的序号n分别保存于矩阵/>第j列的第一行和第二行中，n∈{1，2，…，m}，j∈{1，2，…，N}；最终得到2×N阶矩阵/>
S4：每颗可见卫星的信号频谱附近都会出现的最大值，搜索并记录/>第一行N个元素的最大值max，并确定该最大值max所在列的序号j和对应第二行保存的相位补偿序列的序号n，分别记为j*和n*，j*∈{1，2，…，N}，n*∈{1，2，…，m}；最后计算载波多普勒频率/>
其中，IF表示GNSS中频信号的标称频率。
作为优选，还包括轨迹观测法作为估计卫星信号多普勒频率的辅助方法；所述轨迹观测法中，基于地球静止接收机24小时内的GNSS每颗卫星信号多普勒频率实时检测数据，预先制作GNSS每颗卫星信号绕地球运行一周期内的多普勒频率曲线，接收机捕获时根据接收机时钟计算当前所处的周期内时间，进而根据所述多普勒频率曲线估计出该颗卫星信号的多普勒频率
将多普勒频率和多普勒频率/>进行融合，得到最终的多普勒频率估计值：
其中：τ＝SNR/SNRB，τ≤1，SNRB为空旷地方GPS接收机采样输出端的信噪比，SNR为实际中频采样信号IGIFS的信噪比。
作为优选，设rn(n)为不含GPS信号时的中频采样信号，r(n)为含GPS信号时的中频采样信号，则信噪比SNR采用如下公式计算：
式中E表示求期望，Var表示求方差。
作为优选，当卫星信号极弱导致块叠加与相位补偿模块无法检测并估计GNSS信号多普勒频率的时候，完全通过所述轨迹观测法估计载波多普勒频率，即设此时的参数τ＝0。
作为优选，所述多普勒频率曲线的制作方法为：从GNSS每颗卫星信号绕地球运行一周期内的多普勒频率变化数据中选取6～12个取值变化最为明显的数据点，做分段线性内插，得到建模形成的多普勒频率曲线。
作为优选，所述多普勒频率的估计结果存储在查找表Acqb中；估计完多普勒频率后，开始执行GNSS信号捕获算法。
作为优选，所述查找表Acqb中共有10条记录。
另一方面，本发明提供了一种长时间相干积分捕获的多普勒频率估计与补偿系统，其包括依次连接的块叠加与相位补偿模块、平方模块、谱峰检测模块、轨迹观测模块、存储模块和查表模块；所述块叠加与相位补偿模块中对输入GNSS的中频数据进行分块和叠加操作，压缩载波频率搜索空间，其模块输出端的信号输入平方模块中去掉导航电文后再输入谱峰检测模块；所述谱峰检测模块中，首先对输入信号进行傅里叶变换将其转换到频域上，然后通过搜索谱峰最大值或极大值的方式检测卫星信号，进而估计出存于卫星信号中的多普勒频率，所述谱峰检测模块的输出端连接存储子模块，用于存储估计得到的多普勒频率；所述存储子模块输出端连接查表子模块的输入端，用于查找估计得到的多普勒频率；所述查表子模块输出端连接本地载波生成器的输入端，根据查表得到的多普勒频率估计值生成本地载波信号。
本发明相对于现有技术而言，具有以下有益效果：
本发明可以在信号受到遮挡或一定的环境干扰下保持正常接收GNSS信号，通过采用分块叠加和相位补偿、轨迹观测法等先进技术，提高GNSS信号多普勒频率估计的灵敏度和估计精度(误差小于50Hz)，提高GNSS信号高灵敏度捕获的效率与灵敏度。其中，1)可在较低信噪比环境下实现多普勒频率估计，当输入80毫秒长度的中频采样数据时，即使信噪比降低15分贝仍可以估计出GNSS卫星信号中的多普勒频率；2)提高GNSS信号捕获灵敏度17分贝以上，当信噪比降低17dB时，长时间预相干积分捕获算法可以捕获到4颗以上的GNSS卫星信号；3)GNSS信号捕获所需时间显著缩短，当参与捕获的数据时长为80毫秒时，长时间预相干积分捕获算法仅需运行750秒左右，而公认高效率的非相干累加捕获算法NCI至少需要1400秒的运行时间。
附图说明
图1为长时间相干积分捕获的多普勒频率估计与补偿方法流程图。
图2为仰角计算示意图。
图3为卫星绕地球一周对应的多普勒频率曲线(北纬30°，东经120°)。
图4为多普勒频率周期变化曲线的建模及应用流程。
图5为多普勒频率估计算法输出结果，4中横轴坐标为多普勒频率，单位为Hz，纵坐标为预相干积分的输出。
图6为多普勒频率估计对捕获灵敏度的影响，图中横轴坐标表示GNSS卫星的PRN编号，纵坐标为捕获检测输出，“CIT＝60ms”表示预相干积分时长是60毫秒，“no DBST”指仅延长预相干积分时间而未考虑导航电文比特翻转影响的情况下进行捕获得到的结果，“withDBST”指消除了导航电文比特翻转影响后进行捕获得到的结果。
图7为捕获算法运行时间，图中横轴坐标为参与捕获的数据时长，单位是毫秒，纵坐标为捕获算法运行时间，单位是秒。
具体实施方式
下面结合附图和具体实施方式对本发明做进一步阐述和说明。本发明中各个实施方式的技术特征在没有相互冲突的前提下，均可进行相应组合。
为了克服现有GPS接收机灵敏度低和动态性能低而不能给出定位结果的不足，本发明提供了一种辅助实现GNSS信号高灵敏度快速捕获的多普勒频率估计与补偿方法和系统，使用该方法和系统的接收机可提高GNSS信号捕获速度及其捕获灵敏度，提高GNSS弱信号的检测和锁定能力，并且，当信号受到遮挡、环境噪声较大时也可稳定的捕获卫星信号。
如图1所示，作为本发明的一种实现形式，提出了一种长时间相干积分捕获的多普勒频率估计与补偿方法，其步骤如下：
S1：获取长度为TIms的GNSS零中频数据IGIFS作为输入信号r(·)。
S2：选取一组调整频率δf，得到一组共m条相位补偿序列Ts为采样周期，上标j为虚数单位；对于其中任意第n条相位补偿序列β(n)，n＝1，2，…m，分别执行S21～S23的计算过程：
S21：将所述输入信号r(·)与相位补偿序列值β(n)相乘，得到乘积结果然后根据伪码周期将/>分为M块，每块N个采样数据，M块数据块叠加表示为：
其中，表示/>中第i+1块数据块，r(n)表示/>中第1块数据块，v(n)为噪声，fd表示多普勒频率；HLN(fd)为分块叠加操作的响应函数，计算公式为：
其中：Ts为采样间隔，上标j为虚数单位；
S22：用平方运算去除结果rLN(n)中调制的伪随机码PRN码，得到分块信号数据叠加结果
S23：通过快速傅立叶变换将分块信号数据叠加结果转换到频域上，得到其中f表示自变量为频率；
S3：分别针对m条相位补偿序列执行S21～S23后，得到m个N维向量分别记为/>m个N维向量视为一个m×N维矩阵；然后，分别遍历m个向量/>中的每一维，将m个向量/>中第j维(即m×N维矩阵的第j列)的最大值以及该最大值所属向量/>的序号n分别保存于矩阵/>第j列的第一行和第二行中，n∈{1，2，…，m}，j∈{1，2，…，N}；最终得到2×N阶矩阵/>
上述过程中，m条相位补偿序列共得到m个N维向量，可以视为一个m×N维矩阵。因此上述获得2×N阶矩阵的过程，等价于分别遍历m行N列矩阵中的每一列j，将m维向量中的最大值以及该最大值所对应的位序号n分别保存于一个2×N维矩阵的第j列的两个元素中。
S4：每颗可见卫星的信号频谱附近都会出现的最大值，搜索并记录/>第一行N个元素的最大值max，并确定该最大值max所在列的序号j和对应第二行保存的相位补偿序列的序号n，分别记为j*和n*，j*∈{1，2，…，N}，n*∈{1，2，…，m}；最后计算载波多普勒频率/>
其中，IF表示GNSS中频信号的标称频率。
多普勒频率的估计结果存储在查找表Acqb中(一般可以假设共有10条记录)。估计完多普勒频率后，开始执行GNSS信号捕获算法。
需注意的是，在GNSS接收机捕获过程中上述长时间相干积分捕获的多普勒频率估计与补偿方法可以只运行一次，结果存在表格中供捕获后续卫星信号时查找，而得到的多普勒频率的估计结果，可以在采用多普勒频率估计-压缩载波频率搜索空间-快速捕获算法的方式下，辅助GNSS信号捕获，其计算流程为：捕获卫星信号时，从表Acqb中取出一条记录，由该记录得到多普勒频率/>(或IGIFS的频率)估计值，辅助快速捕获GNSS信号。
另外，当卫星信号极弱导致块叠加与相位补偿模块无法检测并估计GNSS信号多普勒频率的时候，上述长时间相干积分捕获的多普勒频率估计与补偿方法可能会失效，因此本发明中进一步设计了一种轨迹观测建模法用于作为估计卫星信号多普勒频率的辅助方法。
在该轨迹观测法中，其卫星信号多普勒频率估计方法如下：基于地球静止接收机24小时内的GNSS每颗卫星信号多普勒频率实时检测数据，预先制作GNSS每颗卫星信号绕地球运行一周期内的多普勒频率曲线，接收机捕获时根据接收机时钟计算当前所处的周期内时间，进而根据所述多普勒频率曲线估计出该颗卫星信号的多普勒频率
然后，将多普勒频率和多普勒频率/>进行融合，得到最终的多普勒频率估计值：
其中：τ＝SNR/SNRB，τ≤1，SNRB为空旷地方GPS接收机采样输出端的信噪比，SNR为实际中频采样信号IGIFS的信噪比。
其中，设rn(n)为不含GPS信号时的中频采样信号，r(n)为含GPS信号时的中频采样信号，则信噪比SNR可以采用如下公式计算：
式中E表示求期望，Var表示求方差。
当卫星信号极弱导致块叠加与相位补偿模块无法检测并估计GNSS信号多普勒频率的时候，完全通过所述轨迹观测法估计载波多普勒频率，即设此时的参数τ＝0。
另外，前述的多普勒频率曲线需要预先制作，本发明中给出了一种制作方法如下：从GNSS每颗卫星信号绕地球运行一周期内的多普勒频率变化数据中选取6～12个取值变化最为明显的数据点，做分段线性内插，进而得到建模形成的多普勒频率曲线。该多普勒频率曲线可进行存储，以便于后续调用。
基于上述长时间相干积分捕获的多普勒频率估计与补偿方法，还可以进一步提供一种实现该方法的系统，即一种长时间相干积分捕获的多普勒频率估计与补偿系统，其包括依次连接的块叠加与相位补偿模块、平方模块、谱峰检测模块、轨迹观测模块、存储模块和查表模块；所述块叠加与相位补偿模块中对输入GNSS的中频数据进行分块和叠加操作，压缩载波频率搜索空间，其模块输出端的信号输入平方模块中去掉导航电文后再输入谱峰检测模块；所述谱峰检测模块中，首先对输入信号进行傅里叶变换将其转换到频域上，然后通过搜索谱峰最大值或极大值的方式检测卫星信号，进而估计出存于卫星信号中的多普勒频率，所述谱峰检测模块的输出端连接存储子模块，用于存储估计得到的多普勒频率；所述存储子模块输出端连接查表子模块的输入端，用于查找估计得到的多普勒频率；所述查表子模块输出端连接本地载波生成器的输入端，根据查表得到的多普勒频率估计值生成本地载波信号。
下面将上述方法应用于一个具体实施例中，并通过实施例来展示上述方法和系统的具体实现方式以及技术效果。
实施例
在本实施例中，基于多普勒频率估计与补偿算法模块辅助的软件无线电的GNSS(全球导航卫星系统)接收机长时间相干积分捕获算法模块结构如图1所示。捕获算法模块与GNSS接收机的射频前端电路连接，多普勒频率估计与补偿算法模块位于捕获算法的前端，即载波频率搜索的前面，包括依次连接的块叠加与相位补偿模块、平方模块和谱峰检测模块和存储模块和查表模块、轨迹观测法模块。其中，块叠加与相位补偿模块对输入GNSS中频数据进行分块和叠加操作，压缩载波频率搜索空间，其模块输出端信号输入平方模块，去掉导航电文后输入谱峰检测模块，该模块的主要作用是检测相关峰，估计出载波多普勒频率，首先对平方后的信号进行傅里叶变换转换到频域上，然后通过搜索谱峰最大/极大值的方式检测卫星信号，进而计算出存于卫星信号中的多普勒频率，所述谱峰检测模块的输出端连接存储子模块，用于存储估计的多普勒频率；所述存储子模块输出端连接查表子模块的输入端，所述查表子模块输出端连接本地载波生成器的输入端，根据多普勒频率估计值生成本地载波信号。
轨迹观测法模块的输出端接入查表模块，由捕获算法运行时调用，查找当前时刻的天顶可见卫星及其信号中的多普勒频率估计值，用于捕获算法生成本地载波信号。轨迹观测建模法模块作为估计卫星信号多普勒频率的备选方法，基于地球静止接收机24小时的时时检测GNSS每颗卫星信号多普勒频率的数据，制作GNSS每颗卫星信号绕地球运行一周期内的多普勒频率曲线，接收机捕获时根据接收机时钟计算当前所处的周期内时间，进而估计出该颗卫星信号的多普勒频率。该方法一方面可为捕获算法提供捕获搜索的初始参数赋值，包括当前天顶可见卫星PRN编号及每颗卫星信号的多普勒频率预测值，另一方面，该方法在当卫星信号极弱导致块叠加与相位补偿模块无法检测并估计GNSS信号多普勒频率的时候用于估计载波多普勒频率。
整个GNSS信号捕获算法的处理流程为：以GNSS中频采样信号作为输入信号，利用多普勒频率估计与补偿算法模块的多普勒频率估计值，生成不同相位的两路正交本地载波，与输入信号相乘，产生I支路信号和与其正交的Q支路信号。然后I支路和Q支路合为一路复输入信号并进行傅立叶变换，与本地C/A码经过共轭傅立叶变换结果相乘，其结果通过反傅立叶变换转换到时域，再取绝对值，得到输入信号与本地信号间的相关值，最后通过搜索最大相关值判断是否捕获到信号。
为了进一步阐述上述估计系统内具体的多普勒频率估计与补偿算法，下面结合图3对其算法的具体实现方案和步骤详述如下。
多普勒频率估计与补偿算法的处理流程为：以GNSS中频采样信号IGIFS作为输入信号，利用多普勒频率估计与补偿算法模块的多普勒频率估计值压缩载波频率搜索空间，具体步骤如下。
设输入信号为TIms长的零中频采样数据(表示为IGIFS信号)，并设r(·)表示此输入的中频采样信号IGIFS。选取一组调整频率δf，得到一组相位补偿序列其中，Ts为采样周期。对于其中任意第n条相位补偿序列β(n)，n＝1，2,…m，分别执行1)～3)的计算过程：
1)GNSS中频采样数据IGIFS与序列β(n)相乘，乘积结果用表示。然后，根据卫星信号PRN码周期将结果分块，每块N个采样数据，共分为M块。最后，对数据块进行累加(每个数据块中的每个对应位相加)，得到结果/>
其中，表示/>中第i+1块数据块，r(n)表示/>中第1块数据块，v(n)为噪声，fd表示多普勒频率；HLN(fd)为分块叠加操作的响应函数，计算公式为：
HLN(fd)即为分块叠加操作的响应函数。由上式可知，当fdNTs为整数时，即fd取值1KHz的整数倍时(NTs＝1ms)，HLN(fd)取最大值M。这样，搜索GNSS信号多普勒频率的范围由±10KHz变为1KHz了，而且，M个数据分段累加将可增大信噪比(10lgM)dB。
2)用平方运算去除结果中调制的伪随机码PRN码，得到分块信号数据叠加结果/>
3)通过快速傅立叶变换将分块信号数据叠加结果转换到频域上，得到其中f表示自变量为频率
分别针对m条相位补偿序列执行1)～3)后，得到m个N维向量分别记为这m个向量可以看作一信号矩阵。然后，分别遍历m个向量/>中的每一维，将m个向量/>中第j维(视为m×N维矩阵第i列)的最大值以及该最大值所属向量/>的序号n保存于2×N阶矩阵/>的第j列中，即分别遍历m行N列矩阵中的每一列j，将m维向量中的最大值以及该最大值所对应的位序号n分别保存于一个2×N维矩阵的第j列第一行和第二行的两个元素中；n∈{1，2，…，m}，j∈{1，2，…，N}。
每颗可见卫星的信号频谱附近，第一行都会出现最大值(需要注意的是，该最大值并不一定完全是最大值，还有可能是次大值，因此实质为最大/次大值)，该最大/次大值对应的频率的一半即为载波多普勒频率估计值/> 中具有2个N维维度，其中一个记录了最大向量值，而另一个记录了最大向量值所属向量/>的序号n，对应相位补偿序列的序列号。搜索并记录/>中记录前述最大向量值的维度中N个元素的最大值并确定该最大值/>对应的序号n和所属维度j，分别记为n*和j*，j*∈{1，2，…，N}，n*∈{1，2，…，m}；最后计算载波多普勒频率/>
其中，IF表示GNSS中频信号的标称频率。
上述方法获得的载波多普勒频率记录在查找表Acqb中(一般可以假设共有10条记录)。估计完多普勒频率后，即可设置GNSS信号捕获算法开始。但是需要注意的是，如果信号受到遮挡、环境噪声较大时，上述方法可能无法稳定的捕获卫星信号。因此本实施例中进一步设计了轨迹观测法预测多普勒频率。
由于卫星绕地球的运动具有周期性，因此，当接收机相对静止的条件下，卫星信号中的多普勒频率的变化也具有周期性。而且，多普勒随时间变化的曲线具有单值性，因此，通过建模多普勒随时间周期变化的曲线，只要计算出当前时刻对应的多普勒频率曲线中的具体时刻，就可方便的估计多普勒频率的取值。首先，介绍如何计算卫星运动引起的卫星信号中的多普勒频率，以便建模多普勒频率周期变化曲线。下面阐述其具体原理。
卫星位置计算：
卫星在地心地固坐标系中的位置(xs，ys，zs)，可由下列公式表示：
式中，Ω为升交点赤经，i为轨道倾角，ω为近地点角距，f为卫星的真近点角，r为卫星到地心的距离，由于地球的自转，地心地固坐标系也有一个旋转角，表示为：
ΩG＝ΩG0+ωe(t-t0)
其中，t0为参考时间，ΩG0为参考时间的格林威治子午线的赤经，ωe为地球自转角速度，取ωe＝7.292115147e-5(rad/s)。
卫星速度计算：
速度矢量表达式为：
其中和/>分别为：
可见星预测计算：
设接收机位置为接收机记录的上一次工作时记录的位置，在已知接收机位置的情况下，选择仰角大于00的可见卫星，方法如下所述。设地心o、接收机u、卫星s构成三角形Δous(如图2所示)
三角形Δous周长的一半为：
s＝(ou+us+os)/2 (6)
内切圆半径为：
rn＝[(s-ou)(s-os)(s-us)/s]1/2 (7)
则仰角为：
由于仰角太小时，电离层和对流层延迟等误差源所造成的测量误差相对较大，故取仰角大于50的卫星作为可见星。
多普勒频率的计算，计算公式为：
Δfd＝vscosθ·fs/c (9)
式中，c为光速，fs为GPS卫星的载波频率；vs为卫星运动速度(考虑与地球的相对运动)；θ为卫星与接收机的连线和卫星速度vs的夹角。
由于GPS卫星绕地球旋转的周期为11h58min(11小时58分钟)，决定了每颗GPS卫星相对于地球某点的运行周期约为23h56min。所以估计GPS卫星的多普勒频率特性只需估计某地任一23h56min，该地其他时间GPS卫星的多普勒频率特性可由此推算。
按照上述方法，通过仿真分析得到多普勒频率随时间变化的曲线。另外，接收机相对地球静止且已捕获的条件下，由跟踪环路估计得到多普勒频率，记录这些卫星信号中的多普勒频率，本实施例中经多项式拟合得到多普勒频率随时间变化的曲线如图3。根据图3所示，由卫星运动引起的多普勒频率随时间变化具有明显的规律性，可以采用6-12点数据的分段线性内插对多普勒频率周期变化曲线建模，其存储和计算都很简单，对捕获算法复杂度的影响可以忽略不计。具体方法为：在1周期的多普勒频率变化数据中选取6～12个取值变化最为明显的数据点，做分段线性内插，将建模结果存储起来。
多普勒频率周期变化曲线建模后，用于预测多普勒频率时，首先要确定当前接收机时刻所对应的周期变化曲线上的时刻点te，处理方法如下所述。
每周做一次校准。在正常开机并跟踪到信号的情况下，每隔400秒记录一次跟踪到的6-8颗卫星信号的多普勒频率，记录9点数据后，与时间段[ta，tb]内的频移曲线做比较，确定当前时刻位于多普勒频率周期曲线上的具体时刻te。
预测多普勒频率。当长时间丢星后需要重新捕获时，或关机一段时间重新启动并捕获时，可根据上次记录的时刻te(标记为)、记录的丢星至重新捕获这段时间的时长Te，计算出当前重捕时刻位于多普勒频率周期变化曲线上的具体时刻/>查找所有卫星信号的多普勒频率周期变化曲线簇，得到当前时刻下的天顶可见卫星编号及其对应的多普勒频率。
依据如上所述得到的时刻点te即可计算得到多普勒频率的预测值。设得到的GPS信号多普勒频率为Δfd，则GPS信号捕获时，载波频率的搜索以频点(fs+Δfd)为中心展开，伪码相位的搜索以时延(tc+Δtd+Δfd/a)为中心展开。其中，tc为卫星信号发射信号中的伪码相位，Δtd为GPS信号由卫星传播到接收机的传播时延引入的伪码相位变化，Δfd/a为伪码多普勒频率，α为伪码速率与载波中心频率之比。
可见，采用本发明方法可有效缩小捕获搜索的空间，提高捕获速度。经分析可知，在半径800公里左右的范围内，绝大部分时间内，多普勒频率预测误差小于300Hz。因此，可在较大地理范围内使用相同的初始化多普勒频率周期变化曲线。接收机工作一段时间后，逐渐修正多普勒频率周期变化曲线使其更加精确。修正后的多普勒频率估计误差许多情况下小于50Hz。可见，本发明所提方法将有助于缩小多普勒频率搜索空间，提高捕获速度。
当上述轨迹观测法构建完毕后，即可将分块叠加与相位补偿法-轨迹观测法相结合，实现多普勒频率估计与补偿。
轨迹观测法预测多普勒频率，一方面可为捕获算法提供捕获搜索的初始参数赋值，包括当前天顶可见卫星PRN编号及每颗卫星信号的多普勒频率预测值。这样，GPS信号捕获算法无需搜索所有24颗卫星信号，只需搜索当前时刻预测的天顶可见卫星的信号。另一方面，根据多普勒频率的预测值反推其对应频点的坐标值j，搜索矩阵/>的第j列，选出不同i值下的最大/>并按如下公式计算多普勒频率：
设空旷地方GPS接收机采样输出端的信噪比为SNRB，实际中频采样信号IGIFS的信噪比为SNR，设τ＝SNR/SNRB，并限制τ≤1，则多普勒频率估计值为：
设rn(n)为不含GPS信号时的中频采样信号，r(n)为含GPS信号时的中频采样信号，则信噪比采用如下方法计算：
当卫星信号极弱导致块叠加与相位补偿模块无法检测并估计GNSS信号多普勒频率的时候，轨迹观测建模法用于估计载波多普勒频率，即设此时的参数τ＝0。
同样的，上述融合后的多普勒频率估计结果也存储在查找表Acqb中(一般可以假设共有10条记录)。估计完多普勒频率后，GNSS信号捕获算法开始。
上述方案采用多普勒频率估计-压缩载波频率搜索空间-快速捕获算法的方式下，多普勒频率估计与补偿算法模块辅助GNSS信号捕获的计算流程为：捕获卫星信号时，从表Acqb中取出一条记录，由该记录得到多普勒频率(或IGIFS的频率)估计值，然后，依据如下步骤捕获GNSS信号：
(a)依据生成本地载波，并与输入信号相乘，然后做傅里叶变换FFT；
(b)生成sat号卫星对应的本地伪码信号并做共轭傅里叶变换FFT；
(c)步骤(a)、(b)结果乘积后做逆傅里叶变换IFFT，搜索结果中的峰值，进行门限判决。
除频率不确定度范围由±10KHz降低到50Hz外，其余步骤与无多普勒频率估计与补偿算法辅助的GNSS信号捕获算法一样。因此，通过多普勒频率估计与补偿算法估计多普勒频率取值，捕获算法的复杂度可明显降低。GNSS信号捕获算法的功能结构如图1所示。
本专利的有益效果是，可以在信号受到遮挡或一定的环境干扰下保持正常接收GNSS信号，通过采用分块叠加和相位补偿等先进技术，提高GNSS信号多普勒频率估计的灵敏度和估计精度(误差小于50Hz)，提高GNSS信号高灵敏度捕获的效率与灵敏度。其中，1)可在较低信噪比环境下实现多普勒频率估计，当输入80毫秒长度的中频采样数据时，即使信噪比降低15分贝仍可以估计出GNSS卫星信号中的多普勒频率，如图5所示；2)提高GNSS信号捕获灵敏度17分贝以上，如图6所示，当信噪比降低17dB时，长时间预相干积分捕获算法可以捕获到4颗以上的GNSS卫星信号；3)GNSS信号捕获所需时间显著缩短，如图7所示，当参与捕获的数据时长为80毫秒时，长时间预相干积分捕获算法仅需运行750秒左右，而公认高效率的非相干累加捕获算法NCI至少需要1400秒的运行时间。
使用本发明上述频率压缩和估计的方法缩小捕获算法频率搜索空间，进而提高GNSS信号捕获算法的运行效率，当预相干积分时长高于10毫秒时，捕获所需时长一般低于原算法运行时间的五分之一，接近公认高效率算法-非相干累加捕获算法NCI-运行时间的二分之一。
在matlab环境下运行GNSS信号捕获算法，仿真分析结果如下表所述。
表1捕获算法运行时间
备注：在Darius Plausinaitis和Dennis M.Akos等人开发的SoftGNSS v3.0版GPS软件接收机平台上，基于matlab 2014b编程环境。仿真结果中改进前的算法即为SoftGNSSv3.0中的捕获算法，采用基于FFT的频域并行捕获算法。处理的GPS信号来自网上的开源GPS接收机采样数据。捕获算法输入的中频GPS卫星信号的频率为9.548MHz，采样频率为38.192MHz，设置软件接收机的信号通道个数为10，多普勒频率的搜索范围为±7kHz。
以上所述的实施例只是本发明的一种较佳的方案，然其并非用以限制本发明。有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型。因此凡采取等同替换或等效变换的方式所获得的技术方案，均落在本发明的保护范围内。
Claims (8)
1.一种长时间相干积分捕获的多普勒频率估计与补偿方法，其特征在于，步骤如下：
S1：获取长度为TIms的GNSS零中频数据IGIFS作为输入信号r(·)；
S2：选取一组调整频率δf，得到一组共m条相位补偿序列Ts为采样周期，上标j为虚数单位；对于其中任意第n条相位补偿序列β(n)，n＝1，2，…m，分别执行S21～S23的计算过程：
S21：将所述输入信号r(·)与相位补偿序列值β(n)相乘，得到乘积结果然后根据伪码周期将/>分为M块，每块N个采样数据，M块数据块叠加表示为：
其中，表示/>中第i+1块数据块，r(n)表示/>中第1块数据块，v(n)为噪声，fd表示多普勒频率；HLN(fd)为分块叠加操作的响应函数，计算公式为：
S22：用平方运算去除结果rLN(n)中调制的伪随机码PRN码，得到分块信号数据叠加结果
S23：通过快速傅立叶变换将分块信号数据叠加结果转换到频域上，得到其中f表示自变量为频率；
S3：分别针对m条相位补偿序列执行S21～S23后，得到m个N维向量分别记为m个N维向量视为一个m×N维矩阵；然后，分别遍历m个向量/>中的每一维，将m个向量/>中第j维即m×N维矩阵第j列的最大值以及该最大值所属向量的序号n分别保存于矩阵/>第j列的第一行和第二行中，n∈{1,2，…，m}，j∈{1，2，…，N}；最终得到2×N阶矩阵/>
S4：每颗可见卫星的信号频谱附近都会出现的最大值，搜索并记录/>第一行N个元素的最大值max，并确定该最大值max所在列的序号j和对应第二行保存的相位补偿序列的序号n，分别记为j*和n*，j*∈{1，2，…，N}，n*∈{1,2,…,m}；最后计算载波多普勒频率/>
其中，IF表示GNSS中频信号的标称频率。
2.如权利要求1所述的长时间相干积分捕获的多普勒频率估计与补偿方法，其特征在于，还包括轨迹观测法作为估计卫星信号多普勒频率的辅助方法；所述轨迹观测法中，基于地球静止接收机24小时内的GNSS每颗卫星信号多普勒频率实时检测数据，预先制作GNSS每颗卫星信号绕地球运行一周期内的多普勒频率曲线，接收机捕获时根据接收机时钟计算当前所处的周期内时间，进而根据所述多普勒频率曲线估计出该颗卫星信号的多普勒频率
将多普勒频率和多普勒频率/>进行融合，得到最终的多普勒频率估计值：
其中：τ＝min{SNR/SNRB，1}，τ≤1，SNRB为空旷地方GPS接收机采样输出端的信噪比，SNR为实际中频采样信号IGIFS的信噪比。
3.如权利要求2所述的长时间相干积分捕获的多普勒频率估计与补偿方法，其特征在于，设rn(n)为不含GPS信号时的中频采样信号，r(n)为含GPS信号时的中频采样信号，则信噪比SNR采用如下公式计算：
式中E表示求期望，Var表示求方差。
4.如权利要求3所述的长时间相干积分捕获的多普勒频率估计与补偿方法，其特征在于，当卫星信号极弱导致块叠加与相位补偿模块无法检测并估计GNSS信号多普勒频率的时候，完全通过所述轨迹观测法估计载波多普勒频率，即设此时的参数τ＝0。
5.如权利要求2所述的长时间相干积分捕获的多普勒频率估计与补偿方法，其特征在于，所述轨迹观测法中的多普勒频率曲线的制作方法为：从GNSS每颗卫星信号绕地球运行一周期内的多普勒频率变化数据中选取6～12个取值变化最为明显的数据点，做分段线性内插，得到建模形成的多普勒频率曲线。
6.如权利要求1所述的长时间相干积分捕获的多普勒频率估计与补偿方法，其特征在于，所述多普勒频率的估计结果存储在查找表Acqb中；估计完多普勒频率后，开始执行GNSS信号捕获算法。
7.如权利要求6所述的长时间相干积分捕获的多普勒频率估计与补偿方法，其特征在于，所述查找表Acqb中共有10条记录。
8.一种长时间相干积分捕获的多普勒频率估计与补偿系统，其特征在于，包括依次连接的块叠加与相位补偿模块、平方模块、谱峰检测模块、轨迹观测模块、存储模块和查表模块；所述块叠加与相位补偿模块中对输入GNSS的中频数据进行分块和叠加操作，压缩载波频率搜索空间，其模块输出端的信号输入平方模块中去掉导航电文后再输入谱峰检测模块；所述谱峰检测模块中，首先对输入信号进行傅里叶变换将其转换到频域上，然后通过搜索谱峰最大值或极大值的方式检测卫星信号，进而估计出存于卫星信号中的多普勒频率，所述谱峰检测模块的输出端连接存储子模块，用于存储估计得到的多普勒频率；所述存储子模块输出端连接查表子模块的输入端，用于查找估计得到的多普勒频率；所述查表子模块输出端连接本地载波生成器的输入端，根据查表得到的多普勒频率估计值生成本地载波信号。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110198507.2A CN113009523B (zh) | 2021-02-22 | 2021-02-22 | 长时间相干积分捕获的多普勒频率估计与补偿方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110198507.2A CN113009523B (zh) | 2021-02-22 | 2021-02-22 | 长时间相干积分捕获的多普勒频率估计与补偿方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113009523A CN113009523A (zh) | 2021-06-22 |
| CN113009523B true CN113009523B (zh) | 2023-09-29 |
Family
ID=76406400
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110198507.2A Active CN113009523B (zh) | 2021-02-22 | 2021-02-22 | 长时间相干积分捕获的多普勒频率估计与补偿方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113009523B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114465691A (zh) * | 2022-02-15 | 2022-05-10 | 上海兆煊微电子有限公司 | 一种低复杂度的恒包络相位调制信号采样偏差估计和补偿方法及系统 |
| CN115811355B (zh) * | 2023-02-09 | 2023-05-23 | 中国人民解放军战略支援部队航天工程大学 | 一种高动态载波捕获方法 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2009108915A2 (en) * | 2008-02-28 | 2009-09-03 | Magellan Systems Japan, Inc. | Method and apparatus for acquisition, tracking, and sub-microsecond time transfer using weak gps/gnss signals |
| CN103809192A (zh) * | 2014-02-25 | 2014-05-21 | 浙江理工大学 | 一种gnss接收机的动态校正算法 |
| CN104459731A (zh) * | 2014-11-27 | 2015-03-25 | 上海交通大学 | 一种接收极弱gnss信号的高轨卫星定轨方法 |
| US10514466B1 (en) * | 2015-12-07 | 2019-12-24 | Marvell International Ltd. | Method and apparatus for demodulating GNSS navigation data bits under poor clock condition |
| CN111399004A (zh) * | 2020-04-07 | 2020-07-10 | 北京理工大学 | 一种高动态高灵敏度gnss信号捕获方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7477189B2 (en) * | 2007-01-30 | 2009-01-13 | Sirf Technology Holdings, Inc. | Methods and systems for acquisition, reacquisiton and tracking of weak navigational signals |
-
2021
- 2021-02-22 CN CN202110198507.2A patent/CN113009523B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2009108915A2 (en) * | 2008-02-28 | 2009-09-03 | Magellan Systems Japan, Inc. | Method and apparatus for acquisition, tracking, and sub-microsecond time transfer using weak gps/gnss signals |
| CN103809192A (zh) * | 2014-02-25 | 2014-05-21 | 浙江理工大学 | 一种gnss接收机的动态校正算法 |
| CN104459731A (zh) * | 2014-11-27 | 2015-03-25 | 上海交通大学 | 一种接收极弱gnss信号的高轨卫星定轨方法 |
| US10514466B1 (en) * | 2015-12-07 | 2019-12-24 | Marvell International Ltd. | Method and apparatus for demodulating GNSS navigation data bits under poor clock condition |
| CN111399004A (zh) * | 2020-04-07 | 2020-07-10 | 北京理工大学 | 一种高动态高灵敏度gnss信号捕获方法 |
Non-Patent Citations (3)
| Title |
|---|
| Fast GNSS signal acquisition with Doppler frequency estimation algorithm;Faqin Gao et al.;GPS SOLUTION;全文 * |
| GNSS信号中多普勒频移的估计与补偿技术;高法钦 等;无线电工程;第47卷(第04期);全文 * |
| GPS弱信号捕获算法研究;张风国;刘承禹;张红波;;全球定位系统(05);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113009523A (zh) | 2021-06-22 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104076371B (zh) | 利用卫星系统的因特网热点定位 | |
| US8542718B2 (en) | Method and apparatus for acquisition, tracking, and sub-microsecond time transfer using weak GPS/GNSS signals | |
| US8120529B2 (en) | Method and apparatus for autonomous, in-receiver prediction of GNSS ephemerides | |
| US8259012B2 (en) | Software GNSS receiver for high-altitude spacecraft applications | |
| JP3548853B2 (ja) | マルチビット相関器を備えたスペクトル拡散受信機 | |
| US5774829A (en) | Navigation and positioning system and method using uncoordinated beacon signals in conjunction with an absolute positioning system | |
| CN113009523B (zh) | 长时间相干积分捕获的多普勒频率估计与补偿方法及系统 | |
| FI110290B (fi) | Menetelmä informaation vaiheen selvittämiseksi ja elektroniikkalaite | |
| CN108871348B (zh) | 一种利用天基可见光相机的低轨卫星自主定轨方法 | |
| WO2014089988A1 (zh) | 一种高灵敏度北斗辅助授时装置和授时接收机及授时方法 | |
| CN113009524B (zh) | 长时间相干积分捕获的导航电文比特翻转估计方法及系统 | |
| TWI425238B (zh) | 全球導航衛星系統（ｇｎｓｓ）接收器的定位方法 | |
| KR20050087784A (ko) | 위성 기반 위치 지정 시스템 | |
| US20100201569A1 (en) | System and method for detecting location using data communication network | |
| CN102033236A (zh) | 一种卫星导航位置速度联合估计方法 | |
| Altadill et al. | A method for real-time identification and tracking of traveling ionospheric disturbances using ionosonde data: First results | |
| CN104049262A (zh) | 一种基于矢量跟踪环路的波束形成抗干扰方法 | |
| Soloviev et al. | Extending GPS carrier phase availability indoors with a deeply integrated receiver architecture | |
| CN113009522B (zh) | 多普勒频率残差校正的长时间相干积分捕获算法模块 | |
| JPH11183586A (ja) | 特にｇｐｓ型の、完全高周波航行用受信装置 | |
| US20090168851A1 (en) | A correlator sum method for spread spectrum signal receivers | |
| CN105204050B (zh) | 一种惯性辅助的多通道混合型矢量跟踪方法 | |
| WO2018014980A1 (en) | Global navigation satellite system (gnss) signal tracking | |
| Su et al. | A novel GNSS single-frequency PPP approach to estimate the ionospheric TEC and satellite pseudorange observable-specific signal bias | |
| Yang et al. | An adaptive inter-frequency aiding carrier tracking algorithm for the mountain-top GPS radio occultation signal |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |