CN109725290B - Error extraction method and device, electronic equipment and readable storage medium - Google Patents
Error extraction method and device, electronic equipment and readable storage medium Download PDFInfo
- Publication number
- CN109725290B CN109725290B CN201910096259.3A CN201910096259A CN109725290B CN 109725290 B CN109725290 B CN 109725290B CN 201910096259 A CN201910096259 A CN 201910096259A CN 109725290 B CN109725290 B CN 109725290B
- Authority
- CN
- China
- Prior art keywords
- epoch
- coordinate
- day
- satellite
- formula
- 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
- 238000000605 extraction Methods 0.000 title claims abstract description 38
- 238000000034 method Methods 0.000 claims abstract description 59
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 46
- 238000012545 processing Methods 0.000 claims abstract description 45
- 238000000354 decomposition reaction Methods 0.000 claims description 79
- 230000006870 function Effects 0.000 claims description 27
- 238000004364 calculation method Methods 0.000 claims description 20
- 238000004891 communication Methods 0.000 claims description 17
- 238000009499 grossing Methods 0.000 claims description 17
- 238000005259 measurement Methods 0.000 claims description 10
- 238000004590 computer program Methods 0.000 claims description 8
- 230000003203 everyday effect Effects 0.000 claims description 2
- 230000003111 delayed effect Effects 0.000 claims 2
- 125000006850 spacer group Chemical group 0.000 claims 2
- 230000008569 process Effects 0.000 abstract description 15
- 238000005516 engineering process Methods 0.000 abstract description 11
- 239000000284 extract Substances 0.000 abstract description 11
- 230000000694 effects Effects 0.000 description 13
- 230000008859 change Effects 0.000 description 8
- 238000004458 analytical method Methods 0.000 description 6
- 230000002354 daily effect Effects 0.000 description 6
- 238000010586 diagram Methods 0.000 description 6
- 230000000875 corresponding effect Effects 0.000 description 5
- 238000004088 simulation Methods 0.000 description 5
- 238000001914 filtration Methods 0.000 description 4
- 238000004321 preservation Methods 0.000 description 3
- 230000002596 correlated effect Effects 0.000 description 2
- 238000013461 design Methods 0.000 description 2
- 230000002093 peripheral effect Effects 0.000 description 2
- 230000003252 repetitive effect Effects 0.000 description 2
- 101100257981 Arabidopsis thaliana S-ACP-DES3 gene Proteins 0.000 description 1
- 238000012935 Averaging Methods 0.000 description 1
- 101150057677 DAD3 gene Proteins 0.000 description 1
- 102100034004 Gamma-adducin Human genes 0.000 description 1
- 101000799011 Homo sapiens Gamma-adducin Proteins 0.000 description 1
- 101001080624 Homo sapiens Proline/serine-rich coiled-coil protein 1 Proteins 0.000 description 1
- 101000652684 Homo sapiens Transcriptional adapter 3 Proteins 0.000 description 1
- 102100027427 Proline/serine-rich coiled-coil protein 1 Human genes 0.000 description 1
- 101100433696 Saccharomyces cerevisiae (strain ATCC 204508 / S288c) AAD3 gene Proteins 0.000 description 1
- 102100030836 Transcriptional adapter 3 Human genes 0.000 description 1
- 238000010219 correlation analysis Methods 0.000 description 1
- 208000036552 dowling-degos disease 3 Diseases 0.000 description 1
- 208000002153 familial abdominal 3 aortic aneurysm Diseases 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
- 230000001151 other effect Effects 0.000 description 1
- 230000000737 periodic effect Effects 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 238000004445 quantitative analysis Methods 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 230000003313 weakening effect Effects 0.000 description 1
Images
Landscapes
- Position Fixing By Use Of Radio Waves (AREA)
Abstract
本发明实施例提供了一种误差提取方法、装置、电子设备及可读存储介质,应用于无线定位技术领域,所述方法包括:分别获取基准站和观测站在各天预设时间段内的接收信号,对接收信号的伪距及载波相位进行差分处理,根据得到的双差伪距、双差载波相位和基准站的位置坐标,得到观测站在各天的差分坐标序列;通过二维移动加权平均法对差分坐标序列进行处理,得到平滑坐标序列,通过小波包算法对任一参考天以及与该参考天互相关系数最大的邻近天的平滑坐标序列进行M层分解,判断各分解层中低频部分的互相关系数中最大值所在的目标层J是否小于M;如果是,提取目标层J中的低频部分进行小波重构,得到多路径误差。本发明提高了误差提取的准确性。
Embodiments of the present invention provide an error extraction method, device, electronic device, and readable storage medium, which are applied to the field of wireless positioning technology. Receive the signal, perform differential processing on the pseudorange and carrier phase of the received signal, and obtain the differential coordinate sequence of the observation station on each day according to the obtained double-difference pseudorange, double-difference carrier phase and the position coordinates of the base station; The weighted average method processes the difference coordinate sequence to obtain a smooth coordinate sequence. The wavelet packet algorithm is used to decompose the smooth coordinate sequence of any reference day and the adjacent day with the largest cross-correlation coefficient with the reference day. Whether the target layer J where the maximum value of the cross-correlation coefficient of the low-frequency part is located is smaller than M; if so, extract the low-frequency part in the target layer J for wavelet reconstruction to obtain the multipath error. The invention improves the accuracy of error extraction.
Description
技术领域technical field
本发明涉及无线定位技术领域,特别是涉及一种误差提取方法、装置、电子设备及可读存储介质。The present invention relates to the technical field of wireless positioning, and in particular, to an error extraction method, apparatus, electronic device and readable storage medium.
背景技术Background technique
在短基线高精度定位中,通过差分技术可以消除或减弱电离层延迟、对流层延迟、卫星轨道误差、接收机钟差和卫星钟差等误差,但多路径效应在基线两端不具有相关性,无法通过差分消除,因此多路径误差成为高精度定位的主要误差源。In short-baseline high-precision positioning, the ionospheric delay, tropospheric delay, satellite orbit error, receiver clock error, satellite clock error and other errors can be eliminated or weakened by differential technology, but the multipath effect is not correlated at both ends of the baseline. It cannot be eliminated by differential, so multipath error becomes the main error source for high precision positioning.
为了检测、抑制和消除多路径误差,相关学者从卫星信号设计、接收天线的设计与选址、数字信号处理和定位导航计算等方面提出了多种策略。其中,定位导航计算无需在射频前端和基带数字信号处理模块内对接收信号的处理做任何改动,只需对测后数据进行处理。已经提出的测后数据处理方法包括仰角加权技术、扩展卡尔曼滤波技术、利用信噪比技术、经验模式分解技术、Vondark滤波技术以及小波技术等。In order to detect, suppress and eliminate multipath errors, relevant scholars have proposed various strategies from the aspects of satellite signal design, receiving antenna design and location selection, digital signal processing, and positioning and navigation calculations. Among them, the positioning and navigation calculation does not need to make any changes to the processing of the received signal in the RF front-end and the baseband digital signal processing module, but only needs to process the post-measurement data. The post-measurement data processing methods that have been proposed include elevation angle weighting technology, extended Kalman filter technology, signal-to-noise ratio technology, empirical mode decomposition technology, Vondark filter technology and wavelet technology.
其中,小波阈值去噪技术利用小波变换中的变尺度特性,对确定信号有一种“集中”能力。如果一个信号的能量集中于小波变换域少数小波系数上,则它们的取值必然大于在小波变换域内能量分散的大量信号和噪声的小波系数,这时可用阈值去噪的方法。给定一个阈值δ,所有绝对值小于δ的小波系数化为“噪声”,它们的值用0代替;而超过阈值δ的小波系数的数值,被缩减后再重新取值。经小波重构后,得到所需要的信号。但是,利用小波变换的多分辨分析只是对信号的低频部分进行分解,因此提取出的多路径只有低频部分的多路径,而信号高频部分的多路径误差没有进行提取。在利用小波变化对信号的多路径误差进行多分辨分析时,忽略了观测噪声的影响,使得信号滤波后的效果并不理想。因此,现有的小波技术中,提取的多路径误差的准确性较低,滤波后的信号准确性较低。Among them, the wavelet threshold denoising technology utilizes the variable scale characteristic of wavelet transform, and has a "concentration" ability to determine the signal. If the energy of a signal is concentrated on a small number of wavelet coefficients in the wavelet transform domain, their values must be larger than the wavelet coefficients of a large number of signals and noises with energy dispersed in the wavelet transform domain. At this time, the method of threshold denoising can be used. Given a threshold δ, all wavelet coefficients whose absolute value is less than δ are converted into "noise", and their values are replaced by 0; while the values of wavelet coefficients exceeding the threshold δ are reduced and then re-valued. After wavelet reconstruction, the required signal is obtained. However, the multi-resolution analysis using wavelet transform only decomposes the low-frequency part of the signal, so the extracted multipath is only the multi-path of the low-frequency part, and the multi-path error of the high-frequency part of the signal is not extracted. When the multi-resolution analysis of the multi-path error of the signal is carried out by using the wavelet change, the influence of the observation noise is ignored, so that the effect of the signal after filtering is not ideal. Therefore, in the existing wavelet technology, the accuracy of the extracted multipath errors is low, and the accuracy of the filtered signal is low.
发明内容SUMMARY OF THE INVENTION
本发明实施例的目的在于提供一种误差提取方法、装置、电子设备及可读存储介质,以提高误差提取的准确性。具体技术方案如下:The purpose of the embodiments of the present invention is to provide an error extraction method, apparatus, electronic device, and readable storage medium, so as to improve the accuracy of error extraction. The specific technical solutions are as follows:
本发明实施例提供了一种误差提取方法,所述方法包括:An embodiment of the present invention provides an error extraction method, and the method includes:
分别获取基准站和观测站在各天预设时间段内的接收信号,针对所述基准站和所述观测站每天的接收信号,对所述基准站和所述观测站在所述预设时间段内各历元的接收信号的伪距及载波相位进行差分处理,得到双差伪距和双差载波相位;Respectively obtain the received signals of the reference station and the observation station within the preset time period of each day, and for the daily received signals of the reference station and the observation station, the reference station and the observation station for the preset time The pseudorange and carrier phase of the received signal of each epoch in the segment are differentially processed to obtain double-difference pseudorange and double-difference carrier phase;
根据所述基准站的位置坐标、各天的双差伪距和双差载波相位进行位置解算,得到所述观测站在各天的差分坐标序列;Perform position calculation according to the position coordinates of the reference station, the double-difference pseudorange and the double-difference carrier phase of each day, and obtain the differential coordinate sequence of the observation station for each day;
通过二维移动加权平均法对所述观测站在各天的差分坐标序列进行处理,得到所述观测站在各天的平滑坐标序列,所述二维移动加权平均法基于时域空间和值域空间进行移动加权平均;The two-dimensional moving weighted average method is used to process the differential coordinate sequence of the observation station on each day to obtain a smooth coordinate sequence of the observation station on each day. The two-dimensional moving weighted average method is based on the time domain space and the value domain. space for moving weighted average;
针对任一参考天,计算所述观测站该参考天与其他邻近天的平滑坐标序列的互相关系数,选取该参考天以及与该参考天互相关系数最大的邻近天的平滑坐标序列;For any reference day, calculate the cross-correlation coefficient of the smoothed coordinate sequence of the reference day and other adjacent days of the observation station, and select the reference day and the smoothed coordinate sequence of the adjacent day with the largest cross-correlation coefficient with the reference day;
通过小波包算法对该参考天以及与该参考天互相关系数最大的邻近天的平滑坐标序列进行M层分解,计算各分解层中低频部分的互相关系数,M为预先设置的分解层数,判断目标层J是否小于M,所述目标层J为所述各分解层中低频部分的互相关系数中最大值所在的层数;The wavelet packet algorithm is used to decompose the smooth coordinate sequence of the reference day and the adjacent days with the largest cross-correlation coefficient with the reference day in M layers, and the cross-correlation coefficient of the low-frequency part in each decomposition layer is calculated. Judging whether the target layer J is less than M, the target layer J is the layer number where the maximum value is located in the cross-correlation coefficient of the low-frequency part in each decomposition layer;
如果是,提取所述目标层J中的低频部分进行小波重构,得到多路径误差;If so, extract the low-frequency part in the target layer J for wavelet reconstruction to obtain the multi-path error;
如果否,将M的值加1,返回所述通过小波包算法对该参考天以及与该参考天互相关系数最大的邻近天的平滑坐标序列进行M层分解的步骤,直至所述目标层J小于M。If not, add 1 to the value of M, and return to the step of performing M-level decomposition on the smooth coordinate sequence of the reference day and the adjacent day with the largest cross-correlation coefficient with the reference day through the wavelet packet algorithm, until the target layer J less than M.
可选的,所述通过二维移动加权平均法对所述观测站在各天的差分坐标序列进行处理,得到所述观测站在各天的平滑坐标序列,包括:Optionally, the two-dimensional moving weighted average method is used to process the differential coordinate sequence of the observation station on each day to obtain a smooth coordinate sequence of the observation station on each day, including:
若所述样本历元坐标的个数为I,确定所述时域空间的移动平均周期T为不大于I的奇数;If the number of the sample epoch coordinates is 1, it is determined that the moving average period T of the time domain space is an odd number not greater than 1;
若第i个样本历元坐标为xi, 表示值域空间,i为1~I的整数,根据公式:计算所述各样本历元坐标的平均值 If the coordinate of the ith sample epoch is x i , Represents the value range space, i is an integer from 1 to I, according to the formula: Calculate the average of the coordinates of each sample epoch
根据公式:计算所述各样本历元坐标的标准差σI;According to the formula: Calculate the standard deviation σ I of the coordinates of each sample epoch;
将所述标准差σI作为选取值域空间历元坐标的阈值θI,则θI=σI,Taking the standard deviation σ I as the threshold θ I for selecting the spatial epoch coordinates of the range, then θ I =σ I ,
若对第r个样本历元坐标进行平滑处理,r为1~I的整数;If smoothing is performed on the coordinates of the rth sample epoch, r is an integer from 1 to I;
根据公式:计算值域空间内各样本历元坐标与当前样本历元坐标xr的标准差σr,According to the formula: Calculate the standard deviation σ r of each sample epoch coordinate and the current sample epoch coordinate x r in the range space,
将标准差σr不大于阈值θI的样本历元坐标作为分权历元坐标xk,所述分权历元坐标xk用于确定在值域空间第r个样本历元坐标平滑后的值,分权历元坐标xk的个数为第r个样本历元坐标在值域空间的移动平均周期Nr,Nr≤I;对于第r个样本历元坐标,若αj表示时域空间第j个延迟历元坐标的加权系数,j为的整数,[]表示取整数, The sample epoch coordinates whose standard deviation σ r is not greater than the threshold θ I are taken as the decentralized epoch coordinates x k , and the decentralized epoch coordinates x k are used to determine the smoothed r-th sample epoch coordinates in the value domain space. value, the number of weighted epoch coordinates x k is the moving average period N r of the r-th sample epoch coordinates in the range space, N r ≤I; for the r-th sample epoch coordinates, if α j represents The weighting coefficient of the coordinate of the jth delay epoch in the domain space, j is Integer of , [] means to take an integer,
βk表示值域空间第k个分权历元坐标的加权系数,k为1~Nr的整数, β k represents the weighting coefficient of the k-th decentralized epoch coordinate in the range space, k is an integer from 1 to N r ,
根据公式:计算所述时域空间的移动平均周期T内各样本历元坐标相对于第r个样本历元坐标的平均延迟历元坐标为Δτj为第j个样本历元坐标相对于第r个样本历元坐标xr的延迟历元,According to the formula: Calculate the average delay epoch coordinate of each sample epoch coordinate relative to the rth sample epoch coordinate within the moving average period T of the time domain space as Δτ j is the delay epoch of the jth sample epoch coordinate relative to the rth sample epoch coordinate x r ,
根据公式:计算时域空间各延迟历元坐标的异化系数为 根据公式:计算时域空间第j个延迟历元坐标的加权系数αj;According to the formula: The dissimilation coefficient of each delay epoch coordinate in the time domain space is calculated as According to the formula: Calculate the weighting coefficient α j of the coordinate of the jth delay epoch in the time domain space;
根据公式:计算值域空间内分权历元坐标xk的平均值Δμk为第k个分权历元坐标相对于第r个样本历元坐标的差值;According to the formula: Calculate the mean of the weighted epoch coordinates x k in the range space Δμ k is the difference between the coordinate of the k-th decentralized epoch relative to the coordinate of the r-th sample epoch;
根据公式:计算值域空间第k个分权历元坐标的异化系数 According to the formula: Calculate the dissimilation coefficient of the k-th weighted epoch coordinate in the range space
根据公式:计算值域空间第k个分权历元坐标的加权系数βk;According to the formula: Calculate the weighting coefficient β k of the coordinate of the k-th decentralized epoch in the value domain space;
根据公式:According to the formula:
将第r个样本历元坐标xr分别在时域空间和值域空间以xr为中心进行移动加权平均处理,得到第r个样本历元坐标的时域空间移动加权平均值以及值域空间移动平均加权值 Perform the moving weighted average processing of the rth sample epoch coordinate x r in the time domain space and the value domain space with x r as the center, and obtain the time domain space moving weighted average of the rth sample epoch coordinates and the range space moving average weighted value
根据公式:计算xr在二维空间下的平均值 According to the formula: Calculate the mean of x r in 2D space
根据公式:计算时域空间的异化系数 According to the formula: Calculate the dissimilation coefficient in the time-domain space
根据公式:计算值域空间的异化系数 According to the formula: Calculate the dissimilation coefficient of the range space
根据公式:计算时域空间的加权系数γT,根据公式:计算值域空间的加权系数 According to the formula: Calculate the weighting coefficient γ T of the time domain space, according to the formula: Calculate the weighting factor of the range space
根据公式:将和进行加权处理,得到经过平滑处理后的坐标 According to the formula: Will and Perform weighting processing to obtain smoothed coordinates
可选的,所述对所述基准站和所述观测站在所述预设时间段内各历元的接收信号的伪距及载波相位进行差分处理,得到双差伪距和双差载波相位,包括:Optionally, performing differential processing on the pseudorange and carrier phase of the received signals of the reference station and the observation station in each epoch within the preset time period to obtain a double-difference pseudorange and a double-difference carrier phase. ,include:
根据公式:According to the formula:
计算观测站r和基准站b在同一历元观测卫星n和卫星m的双差伪距以及观测站r和基准站b在同一历元观测卫星n和卫星m的双差载波相位 Calculate the double-difference pseudoranges of satellite n and satellite m observed by observation station r and base station b in the same epoch and the double-difference carrier phase of satellite n and satellite m observed by observation station r and reference station b in the same epoch
其中,为双差算子,S表示GPS系统,卫星n和卫星m为该GPS系统内的卫星,表示观测站r和基准站b在同一历元观测卫星n和卫星m时卫星中心到接收机相位中心间的几何距离,λ为载波波长;为观测站r和基准站b在同一历元观测卫星n和卫星m的载波相位未知整周模糊度,表示基准站b和观测站r在同一历元观测卫星n和卫星m的伪距多路径误差,表示观测站r和基准站b在同一历元观测卫星n和卫星m的载波相位多路径误差,表示观测站r和基准站b在同一历元观测卫星n和卫星m的伪码,表示观测站r和基准站b在同一历元观测卫星n和卫星m的载波相位测量噪声。in, is the double difference operator, S represents the GPS system, satellite n and satellite m are the satellites in the GPS system, Represents the geometric distance between the satellite center and the receiver phase center when the observation station r and the base station b observe satellite n and satellite m in the same epoch, and λ is the carrier wavelength; is the unknown integer ambiguity of the carrier phase of satellite n and satellite m observed at the same epoch by observation station r and base station b, represents the pseudorange multipath error of satellite n and satellite m observed by base station b and observation station r in the same epoch, represents the carrier phase multipath error of satellite n and satellite m observed by observation station r and reference station b in the same epoch, Pseudocode indicating that observation station r and base station b observe satellite n and satellite m in the same epoch, It represents the carrier phase measurement noise of satellite n and satellite m observed by observation station r and reference station b in the same epoch.
可选的,所述通过小波包算法对该参考天以及与该参考天互相关系数最大的邻近天的平滑坐标序列进行M层分解,包括:Optionally, performing M-level decomposition on the reference day and the smooth coordinate sequence of the adjacent day with the largest cross-correlation coefficient with the reference day by using the wavelet packet algorithm, including:
若该参考天以及与该参考天互相关系数最大的邻近天的平滑坐标序列为xi,根据公式:If the smooth coordinate sequence of the reference day and the adjacent day with the largest cross-correlation coefficient with the reference day is x i , according to the formula:
对xi进行小波包分解, Perform wavelet packet decomposition on x i ,
i为1~I的整数,ω为1~2M的整数,q∈I,q为参考历元;i is an integer from 1 to I, ω is an integer from 1 to 2M, q∈I, q is a reference epoch;
表示未分解时的小波包,表示第M层上的第2ω-1个小波包系数,表示第M-1层上的第ω个小波包系数,表示第M层上的第2ω个小波包系数,G为尺度函数分解滤波器,H为小波函数分解滤波器, represents the wavelet packet when it is not decomposed, represents the 2ω-1 wavelet packet coefficient on the Mth layer, represents the ωth wavelet packet coefficient on the M-1th layer, represents the second ω wavelet packet coefficient on the Mth layer, G is the scale function decomposition filter, H is the wavelet function decomposition filter,
所述提取所述目标层J中的低频部分进行小波重构,得到多路径误差,包括:The extraction of the low-frequency part in the target layer J for wavelet reconstruction to obtain a multi-path error, including:
根据公式:对xi进行小波包重构;According to the formula: Perform wavelet packet reconstruction on xi ;
表示第J+1层上的第2ω-1个小波包系数,表示第J+1层上的第2ω个小波包系数,g为尺度函数重构滤波器,h为小波函数重构滤波器。 represents the 2ω-1 wavelet packet coefficient on the J+1th layer, represents the second ω wavelet packet coefficient on the J+1th layer, g is the scale function reconstruction filter, and h is the wavelet function reconstruction filter.
本发明实施例提供了一种误差提取装置,所述装置包括:An embodiment of the present invention provides an error extraction device, and the device includes:
双差计算模块,用于分别获取基准站和观测站在各天预设时间段内的接收信号,针对所述基准站和所述观测站每天的接收信号,对所述基准站和所述观测站在所述预设时间段内各历元的接收信号的伪距及载波相位进行差分处理,得到双差伪距和双差载波相位;The double-difference calculation module is used to obtain the received signals of the reference station and the observation station within a preset time period of each day respectively, and for the daily received signals of the reference station and the observation station, the reference station and the observation station Perform differential processing on the pseudorange and carrier phase of the received signal of each epoch in the preset time period to obtain a double-difference pseudorange and a double-difference carrier phase;
差分坐标序列计算模块,用于根据所述基准站的位置坐标、各天的双差伪距和双差载波相位进行位置解算,得到所述观测站在各天的差分坐标序列;A differential coordinate sequence calculation module, configured to perform position calculation according to the position coordinates of the reference station, the double-difference pseudorange and the double-difference carrier phase of each day, and obtain the differential coordinate sequence of the observation station for each day;
平滑处理模块,用于通过二维移动加权平均法对所述观测站在各天的差分坐标序列进行处理,得到所述观测站在各天的平滑坐标序列,所述二维移动加权平均法基于时域空间和值域空间进行移动加权平均;The smoothing processing module is used to process the differential coordinate sequence of the observation station on each day by a two-dimensional moving weighted average method, and obtain the smoothed coordinate sequence of the observation station on each day, and the two-dimensional moving weighted average method is based on Moving weighted average in time domain space and value domain space;
平滑坐标序列选取模块,用于针对任一参考天,计算所述观测站该参考天与其他邻近天的平滑坐标序列的互相关系数,选取该参考天以及与该参考天互相关系数最大的邻近天的平滑坐标序列;The smooth coordinate sequence selection module is used to calculate the cross-correlation coefficient of the smooth coordinate sequence of the reference day and other adjacent days of the observation station for any reference day, and select the reference day and the neighbor with the largest cross-correlation coefficient with the reference day. smooth coordinate sequence of days;
分解模块,用于通过小波包算法对该参考天以及与该参考天互相关系数最大的邻近天的平滑坐标序列进行M层分解,计算各分解层中低频部分的互相关系数,M为预先设置的分解层数;The decomposition module is used to decompose the smooth coordinate sequence of the reference day and the adjacent day with the largest cross-correlation coefficient with the reference day through the wavelet packet algorithm, and calculate the cross-correlation coefficient of the low-frequency part in each decomposition layer. M is a preset The number of decomposition layers;
判断模块,用于判断目标层J是否小于M,所述目标层J为所述各分解层中低频部分的互相关系数中最大值所在的层数;Judging module, for judging whether the target layer J is less than M, and the target layer J is the layer number where the maximum value is located in the cross-correlation coefficient of the low-frequency part in each decomposition layer;
重构模块,用于在所述判断模块的判断结果为是时,提取所述目标层J中的低频部分进行小波重构,得到多路径误差;a reconstruction module, configured to extract the low-frequency part in the target layer J for wavelet reconstruction when the judgment result of the judgment module is yes, to obtain a multi-path error;
循环模块,用于在所述判断模块的判断结果为否时,将M的值加1,返回所述通过小波包算法对该参考天以及与该参考天互相关系数最大的邻近天的平滑坐标序列进行M层分解的步骤,直至所述目标层J小于M。The loop module is used to add 1 to the value of M when the judgment result of the judgment module is no, and return the smooth coordinates of the reference day and the adjacent day with the largest cross-correlation coefficient with the reference day through the wavelet packet algorithm. The sequence performs the steps of M layer decomposition until the target layer J is smaller than M.
可选的,所述平滑处理模块具体用于,若所述样本历元坐标的个数为I,确定所述时域空间的移动平均周期T为不大于I的奇数;Optionally, the smoothing processing module is specifically configured to, if the number of the sample epoch coordinates is 1, determine that the moving average period T of the time domain space is an odd number not greater than 1;
若第i个样本历元坐标为xi, 表示值域空间,i为1~I的整数,根据公式:计算所述各样本历元坐标的平均值 If the coordinate of the ith sample epoch is x i , Represents the value range space, i is an integer from 1 to I, according to the formula: Calculate the average of the coordinates of each sample epoch
根据公式:计算所述各样本历元坐标的标准差σI;According to the formula: Calculate the standard deviation σ I of the coordinates of each sample epoch;
将所述标准差σI作为选取值域空间历元坐标的阈值θI,则θI=σI,Taking the standard deviation σ I as the threshold θ I for selecting the spatial epoch coordinates of the range, then θ I =σ I ,
若对第r个样本历元坐标进行平滑处理,r为1~I的整数;If smoothing is performed on the coordinates of the rth sample epoch, r is an integer from 1 to I;
根据公式:计算值域空间内各样本历元坐标与当前样本历元坐标xr的标准差σr,According to the formula: Calculate the standard deviation σ r of each sample epoch coordinate and the current sample epoch coordinate x r in the range space,
将标准差σr不大于阈值θI的样本历元坐标作为分权历元坐标xk,所述分权历元坐标xk用于确定在值域空间第r个样本历元坐标平滑后的值,分权历元坐标xk的个数为第r个样本历元坐标在值域空间的移动平均周期Nr,Nr≤I;对于第r个样本历元坐标,若αj表示时域空间第j个延迟历元坐标的加权系数,j为的整数,[]表示取整数, The sample epoch coordinates whose standard deviation σ r is not greater than the threshold θ I are taken as the decentralized epoch coordinates x k , and the decentralized epoch coordinates x k are used to determine the smoothed r-th sample epoch coordinates in the value domain space. value, the number of weighted epoch coordinates x k is the moving average period N r of the r-th sample epoch coordinates in the range space, N r ≤I; for the r-th sample epoch coordinates, if α j represents The weighting coefficient of the coordinate of the jth delay epoch in the domain space, j is Integer of , [] means to take an integer,
βk表示值域空间第k个分权历元坐标的加权系数,k为1~Nr的整数, β k represents the weighting coefficient of the k-th decentralized epoch coordinate in the range space, k is an integer from 1 to N r ,
根据公式:计算所述时域空间的移动平均周期T内各样本历元坐标相对于第r个样本历元坐标的平均延迟历元坐标为Δτj为第j个样本历元坐标相对于第r个样本历元坐标xr的延迟历元,According to the formula: Calculate the average delay epoch coordinate of each sample epoch coordinate relative to the rth sample epoch coordinate within the moving average period T of the time domain space as Δτ j is the delay epoch of the jth sample epoch coordinate relative to the rth sample epoch coordinate x r ,
根据公式:计算时域空间各延迟历元坐标的异化系数为 根据公式:计算时域空间第j个延迟历元坐标的加权系数αj;According to the formula: The dissimilation coefficient of each delay epoch coordinate in the time domain space is calculated as According to the formula: Calculate the weighting coefficient α j of the coordinate of the jth delay epoch in the time domain space;
根据公式:计算值域空间内分权历元坐标xk的平均值Δμk为第k个分权历元坐标相对于第r个样本历元坐标的差值;According to the formula: Calculate the mean of the weighted epoch coordinates x k in the range space Δμ k is the difference between the coordinate of the k-th decentralized epoch relative to the coordinate of the r-th sample epoch;
根据公式:计算值域空间第k个分权历元坐标的异化系数 According to the formula: Calculate the dissimilation coefficient of the k-th weighted epoch coordinate in the range space
根据公式:计算值域空间第k个分权历元坐标的加权系数βk;According to the formula: Calculate the weighting coefficient β k of the coordinate of the k-th decentralized epoch in the value domain space;
根据公式:According to the formula:
将第r个样本历元坐标xr分别在时域空间和值域空间以xr为中心进行移动加权平均处理,得到第r个样本历元坐标的时域空间移动加权平均值以及值域空间移动平均加权值 Perform the moving weighted average processing of the rth sample epoch coordinate x r in the time domain space and the value domain space with x r as the center, and obtain the time domain space moving weighted average of the rth sample epoch coordinates and the range space moving average weighted value
根据公式:计算xr在二维空间下的平均值 According to the formula: Calculate the mean of x r in 2D space
根据公式:计算时域空间的异化系数 According to the formula: Calculate the dissimilation coefficient in the time-domain space
根据公式:计算值域空间的异化系数 According to the formula: Calculate the dissimilation coefficient of the range space
根据公式:计算时域空间的加权系数γT,根据公式:计算值域空间的加权系数 According to the formula: Calculate the weighting coefficient γ T of the time domain space, according to the formula: Calculate the weighting factor of the range space
根据公式:将和进行加权处理,得到经过平滑处理后的坐标 According to the formula: Will and Perform weighting processing to obtain smoothed coordinates
可选的,所述双差计算模块,具体用于根据公式:Optionally, the double-difference calculation module is specifically used according to the formula:
计算观测站r和基准站b在同一历元观测卫星n和卫星m的双差伪距以及观测站r和基准站b在同一历元观测卫星n和卫星m的双差载波相位 Calculate the double-difference pseudoranges of satellite n and satellite m observed by observation station r and base station b in the same epoch and the double-difference carrier phase of satellite n and satellite m observed by observation station r and reference station b in the same epoch
其中,为双差算子,S表示GPS系统,卫星n和卫星m为该GPS系统内的卫星,表示观测站r和基准站b在同一历元观测卫星n和卫星m时卫星中心到接收机相位中心间的几何距离,λ为载波波长,为观测站r和基准站b在同一历元观测卫星n和卫星m的载波相位未知整周模糊度,表示基准站b和观测站r在同一历元观测卫星n和卫星m的伪距多路径误差,表示观测站r和基准站b在同一历元观测卫星n和卫星m的载波相位多路径误差,表示观测站r和基准站b在同一历元观测卫星n和卫星m的伪码,表示观测站r和基准站b在同一历元观测卫星n和卫星m的载波相位测量噪声。in, is the double difference operator, S represents the GPS system, satellite n and satellite m are the satellites in the GPS system, represents the geometric distance between the satellite center and the receiver phase center when the observation station r and the base station b observe satellite n and satellite m in the same epoch, λ is the carrier wavelength, is the unknown integer ambiguity of the carrier phase of satellite n and satellite m observed at the same epoch by observation station r and base station b, represents the pseudorange multipath error of satellite n and satellite m observed by base station b and observation station r in the same epoch, represents the carrier phase multipath error of satellite n and satellite m observed by observation station r and reference station b in the same epoch, Pseudocode indicating that observation station r and base station b observe satellite n and satellite m in the same epoch, It represents the carrier phase measurement noise of satellite n and satellite m observed by observation station r and reference station b in the same epoch.
可选的,所述分解模块,具体用于若该参考天以及与该参考天互相关系数最大的邻近天的平滑坐标序列为xi,根据公式:Optionally, the decomposition module is specifically used if the smooth coordinate sequence of the reference day and the adjacent day with the largest cross-correlation coefficient with the reference day is x i , according to the formula:
对xi进行小波包分解, Perform wavelet packet decomposition on x i ,
i为1~I的整数,ω为1~2M的整数,q∈I,q为参考历元;i is an integer from 1 to I, ω is an integer from 1 to 2M, q∈I, q is a reference epoch;
表示未分解时的小波包,表示第M层上的第2ω-1个小波包系数,表示第M-1层上的第ω个小波包系数,表示第M层上的第2ω个小波包系数,G为尺度函数分解滤波器,H为小波函数分解滤波器, represents the wavelet packet when it is not decomposed, represents the 2ω-1 wavelet packet coefficient on the Mth layer, represents the ωth wavelet packet coefficient on the M-1th layer, represents the second ω wavelet packet coefficient on the Mth layer, G is the scale function decomposition filter, H is the wavelet function decomposition filter,
所述重构模块,具体用于根据公式:The reconstruction module is specifically used according to the formula:
对xi进行小波包重构; Perform wavelet packet reconstruction on xi ;
表示第J+1层上的第2ω-1个小波包系数,表示第J+1层上的第2ω个小波包系数,g为尺度函数重构滤波器,h为小波函数重构滤波器。 represents the 2ω-1 wavelet packet coefficient on the J+1th layer, represents the second ω wavelet packet coefficient on the J+1th layer, g is the scale function reconstruction filter, and h is the wavelet function reconstruction filter.
本发明实施例提供了一种电子设备,包括:处理器、通信接口、存储器和通信总线,其中,所述处理器、所述通信接口、所述存储器通过所述通信总线完成相互间的通信;An embodiment of the present invention provides an electronic device, including: a processor, a communication interface, a memory, and a communication bus, wherein the processor, the communication interface, and the memory communicate with each other through the communication bus;
所述存储器,用于存放计算机程序;the memory for storing computer programs;
所述处理器,用于执行所述存储器上所存放的程序时,实现上述任一所述的误差提取方法的步骤。The processor is configured to implement the steps of any one of the above error extraction methods when executing the program stored in the memory.
本发明实施例提供了一种计算机可读存储介质,所述计算机可读存储介质内存储有计算机程序,所述计算机程序被处理器执行时,实现上述任一所述的误差提取方法的步骤。An embodiment of the present invention provides a computer-readable storage medium, where a computer program is stored in the computer-readable storage medium, and when the computer program is executed by a processor, the steps of any of the above-described error extraction methods are implemented.
本发明实施例提供的误差提取方法、装置、电子设备及可读存储介质,分别获取基准站和观测站在各天预设时间段内的接收信号,针对基准站和观测站每天的接收信号,对基准站和观测站在预设时间段内各历元的接收信号的伪距及载波相位进行差分处理,并根据得到的各天的双差伪距、双差载波相位以及基准站的位置坐标进行位置解算,得到观测站在各天的差分坐标序列;通过二维移动加权平均法对差分坐标序列进行处理,得到观测站在各天的平滑坐标序列,二维移动加权平均法基于时域空间和值域空间进行移动加权平均;通过小波包算法对互相关系数最大的参考天以及邻近天的平滑坐标序列进行M层分解,计算各分解层中低频部分的互相关系数,判断目标层J是否小于M,目标层J为各分解层中低频部分的互相关系数中最大值所在的层数;如果是,提取目标层J中的低频部分进行小波重构,得到多路径误差;如果否,将M的值加1,返回通过小波包算法对互相关系数最大的参考天以及邻近天的平滑坐标序列进行M层分解的步骤,直至目标层J小于M;提取目标层J中的低频部分进行小波重构,得到多路径误差。本发明实施例中,通过二维移动加权平均法可以削弱低频部分观测噪声对信号的影响,并解决信号平滑处理中的端部效应,达到保边去噪的目的。小波包算法可以将信号频带进行多层次划分,对高频部分进一步分解,消除高频部分中的观测噪声,从而提高多路径误差提取的准确性。当然,实施本发明的任一产品或方法并不一定需要同时达到以上所述的所有优点。The error extraction method, device, electronic device, and readable storage medium provided by the embodiments of the present invention respectively obtain the received signals of the reference station and the observation station within the preset time period of each day, and for the daily received signals of the reference station and the observation station, Perform differential processing on the pseudorange and carrier phase of the received signal at each epoch of the base station and the observation station within a preset time period, and obtain the double-difference pseudorange and double-difference carrier phase of each day and the position coordinates of the base station. Perform position calculation to obtain the differential coordinate sequence of the observation station on each day; process the differential coordinate sequence through the two-dimensional moving weighted average method to obtain the smooth coordinate sequence of the observation station on each day. The two-dimensional moving weighted average method is based on the time domain. The moving weighted average of space and range space is carried out; M-level decomposition is carried out on the smooth coordinate sequence of the reference day with the largest cross-correlation coefficient and the adjacent days through the wavelet packet algorithm, and the cross-correlation coefficient of the low-frequency part in each decomposition layer is calculated, and the target layer J is determined. Whether it is less than M, the target layer J is the number of layers where the maximum value of the cross-correlation coefficient of the low-frequency part in each decomposition layer is located; if so, extract the low-frequency part in the target layer J for wavelet reconstruction to obtain the multi-path error; if not, Add 1 to the value of M, and return to the step of performing M-level decomposition on the smooth coordinate sequence of the reference day with the largest cross-correlation coefficient and the adjacent days through the wavelet packet algorithm, until the target layer J is less than M; extract the low-frequency part in the target layer J for Wavelet reconstruction to get multipath error. In the embodiment of the present invention, the influence of low-frequency part observation noise on the signal can be weakened by the two-dimensional moving weighted average method, and the end effect in the signal smoothing process can be solved, so as to achieve the purpose of edge preservation and denoising. The wavelet packet algorithm can divide the signal frequency band into multiple layers, further decompose the high-frequency part, and eliminate the observation noise in the high-frequency part, thereby improving the accuracy of multi-path error extraction. Of course, it is not necessary for any product or method of the present invention to achieve all of the advantages described above at the same time.
附图说明Description of drawings
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其他的附图。In order to explain the embodiments of the present invention or the technical solutions in the prior art more clearly, the following briefly introduces the accompanying drawings that need to be used in the description of the embodiments or the prior art. Obviously, the accompanying drawings in the following description are only These are some embodiments of the present invention. For those of ordinary skill in the art, other drawings can also be obtained according to these drawings without creative efforts.
图1为本发明实施例的误差提取方法的流程图;1 is a flowchart of an error extraction method according to an embodiment of the present invention;
图2为本发明实施例中3层分解后的小波包分解树结构图;Fig. 2 is the wavelet packet decomposition tree structure diagram after 3-layer decomposition in the embodiment of the present invention;
图3为通过小波阈值去噪算法、二维移动加权平均算法以及基于二维移动加权平均处理的小波包算法进行滤波去噪后的信号;Fig. 3 is the signal after filtering and denoising by wavelet threshold denoising algorithm, two-dimensional moving weighted average algorithm and wavelet packet algorithm based on two-dimensional moving weighted average processing;
图4为通过小波阈值去噪算法、二维移动加权平均算法以及基于二维移动加权平均处理的小波包算法进行滤波后的残差信号;Fig. 4 is the residual signal after filtering by the wavelet threshold denoising algorithm, the two-dimensional moving weighted average algorithm and the wavelet packet algorithm based on the two-dimensional moving weighted average processing;
图5为本发明实施例的误差提取装置的结构图;5 is a structural diagram of an error extraction device according to an embodiment of the present invention;
图6为本发明实施例的电子设备的结构图。FIG. 6 is a structural diagram of an electronic device according to an embodiment of the present invention.
具体实施方式Detailed ways
下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅仅是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, but not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.
在使用高精度定位技术进行变形观测时,由于观测站一般固定且周围环境基本保持不变,GPS(Global Positioning System,全球定位系统)卫星与观测站以及周围反射面的几何构型具有重复性,以一个恒星日(11小时58分)为周期。在存在多路径效应的情况下,其重复性的特性就使得临近天之间的坐标序列存在相关性。基于此特性,本发明实施例提供了一种误差提取方法、装置、电子设备及可读存储介质,以提高误差提取的准确性,从而削弱多路径效应。When using high-precision positioning technology for deformation observation, since the observation station is generally fixed and the surrounding environment remains basically unchanged, the geometry of the GPS (Global Positioning System, Global Positioning System) satellite, the observation station and the surrounding reflector is repetitive. With a sidereal day (11 hours 58 minutes) as the cycle. In the presence of multipath effects, its repetitive nature makes the coordinate sequence between adjacent days correlated. Based on this characteristic, embodiments of the present invention provide an error extraction method, apparatus, electronic device, and readable storage medium, so as to improve the accuracy of error extraction, thereby weakening the multipath effect.
下面首先对本发明实施例所提供的误差提取方法进行详细介绍。The error extraction method provided by the embodiment of the present invention is first introduced in detail below.
参见图1,图1为本发明实施例的误差提取方法的流程图,包括以下步骤:Referring to FIG. 1, FIG. 1 is a flowchart of an error extraction method according to an embodiment of the present invention, including the following steps:
S101,分别获取基准站和观测站在各天预设时间段内的接收信号,针对基准站和观测站每天的接收信号,对基准站和观测站在预设时间段内各历元的接收信号的伪距及载波相位进行差分处理,得到双差伪距和双差载波相位。S101: Acquire the received signals of the base station and the observation station in the preset time period of each day respectively, and for the daily received signals of the base station and the observation station, obtain the received signals of the base station and the observation station in each epoch in the preset time period. The pseudorange and carrier phase are differentially processed to obtain double-difference pseudorange and double-difference carrier phase.
本发明实施例中,由于观测站天线与接收信号反射面距离微小的变化将导致多路径效应发生很大的变化,且这种变化表现为周期性的变化,将使多路径效应重复的周期发生变化,因此观测站在尽可能小的范围内变化且趋于不变。当反射信号的入射角过大时,多路径效应的重复性将变低甚至消失。因此,反射信号的入射角尽可能的小且小于15°。In the embodiment of the present invention, due to the slight change in the distance between the observation station antenna and the receiving signal reflecting surface, the multipath effect will change greatly, and this change is manifested as a periodic change, which will cause the repeated cycle of the multipath effect to occur. changes, so the observatory changes in as small a range as possible and tends to be constant. When the incident angle of the reflected signal is too large, the repeatability of the multipath effect will decrease or even disappear. Therefore, the incident angle of the reflected signal is as small as possible and less than 15°.
基准站的位置坐标是确定的,观测站的位置坐标是未知的,可以通过基准站和观测站的接收信号确定观测站的位置。具体的,由于临近天之间的坐标序列存在相关性,可以获取基准站和观测站在各天预设时间段内的接收信号,预设时间段可以是一天中任意的时间段,例如8:00-12:00,在此不做限定。对于基准站和观测站每天在预设时间段内的接收信号,可以将预设时间段分为多个历元,历元在天文学是一些天文变数作为参考的时刻点,那么,可以对基准站和观测站在各历元的接收信号的伪距及载波相位进行差分处理,得到双差伪距和双差载波相位。The position coordinates of the base station are determined, and the position coordinates of the observation station are unknown. The position of the observation station can be determined through the received signals of the base station and the observation station. Specifically, due to the correlation between the coordinate sequences between adjacent days, the received signals of the reference station and the observation station in the preset time period of each day can be obtained, and the preset time period can be any time period of the day, such as 8: 00-12:00, not limited here. For the signals received by the base station and the observation station within a preset time period every day, the preset time period can be divided into multiple epochs. In astronomy, an epoch is a point in time when some astronomical variables are used as a reference. Differential processing is performed with the pseudorange and carrier phase of the received signal at each epoch of the observation station, and the double-difference pseudo-range and double-difference carrier phase are obtained.
本发明的一种实现方式中,根据公式:In an implementation manner of the present invention, according to the formula:
计算观测站r和基准站b在同一历元观测卫星n和卫星m的双差伪距以及观测站r和基准站b在同一历元观测卫星n和卫星m的双差载波相位 Calculate the double-difference pseudoranges of satellite n and satellite m observed by observation station r and base station b in the same epoch and the double-difference carrier phase of satellite n and satellite m observed by observation station r and reference station b in the same epoch
其中,为双差算子,S表示GPS系统,卫星n和卫星m为该GPS系统内的卫星,表示观测站r和基准站b在同一历元观测卫星n和卫星m时卫星中心到接收机相位中心间的几何距离,λ为载波波长;为观测站r和基准站b在同一历元观测卫星n和卫星m的载波相位未知整周模糊度,表示基准站b和观测站r在同一历元观测卫星n和卫星m的伪距多路径误差,表示观测站r和基准站b在同一历元观测卫星n和卫星m的载波相位多路径误差,表示观测站r和基准站b在同一历元观测卫星n和卫星m的伪码,表示观测站r和基准站b在同一历元观测卫星n和卫星m的载波相位测量噪声。in, is the double difference operator, S represents the GPS system, satellite n and satellite m are the satellites in the GPS system, Represents the geometric distance between the satellite center and the receiver phase center when the observation station r and the base station b observe satellite n and satellite m in the same epoch, and λ is the carrier wavelength; is the unknown integer ambiguity of the carrier phase of satellite n and satellite m observed at the same epoch by observation station r and base station b, represents the pseudorange multipath error of satellite n and satellite m observed by base station b and observation station r in the same epoch, represents the carrier phase multipath error of satellite n and satellite m observed by observation station r and reference station b in the same epoch, Pseudocode indicating that observation station r and base station b observe satellite n and satellite m in the same epoch, It represents the carrier phase measurement noise of satellite n and satellite m observed by observation station r and reference station b in the same epoch.
S102,根据基准站的位置坐标、各天的双差伪距和双差载波相位进行位置解算,得到观测站在各天的差分坐标序列。S102 , perform position calculation according to the position coordinates of the reference station, the double-difference pseudorange and the double-difference carrier phase of each day, and obtain the differential coordinate sequence of the observation station in each day.
本发明实施例中,可以在RTKLIB(定位结算软件)中进行差分处理,消除或削弱所得差分坐标序列中的电离层延迟、对流层延迟、卫星轨道误差、接收机和卫星钟差等误差。In the embodiment of the present invention, differential processing can be performed in RTKLIB (positioning settlement software) to eliminate or weaken errors such as ionospheric delay, tropospheric delay, satellite orbit error, receiver and satellite clock errors in the obtained differential coordinate sequence.
S103,通过二维移动加权平均法对观测站在各天的差分坐标序列进行处理,得到观测站在各天的平滑坐标序列,二维移动加权平均法基于时域空间和值域空间进行移动加权平均。S103: Process the differential coordinate sequence of the observation station on each day by a two-dimensional moving weighted average method to obtain a smooth coordinate sequence of the observation station on each day, and the two-dimensional moving weighted average method performs moving weighting based on the time domain space and the value domain space average.
具体的,观测站对基准站(周围环境不发生变化)进行静态观测时,其多路径误差在相邻历元间的变化趋于平稳,而由于观测噪声的存在,使得相邻历元间的坐标误差波动加剧。为了削弱差分定位所得坐标序列中的观测噪声,得到精确的多路径误差提取模型,可以通过二维移动加权平均法进行滤波预处理,削弱观测噪声。二维移动加权平均法指的是时域空间和值域空间两个维度的移动加权平均法,下文将对该方法进行详细介绍。Specifically, when the observation station performs static observation on the reference station (the surrounding environment does not change), the change of its multipath error between adjacent epochs tends to be stable, and due to the existence of observation noise, the difference between adjacent epochs is reduced. Coordinate error fluctuations intensify. In order to weaken the observation noise in the coordinate sequence obtained by differential positioning and obtain an accurate multi-path error extraction model, the two-dimensional moving weighted average method can be used for filtering preprocessing to weaken the observation noise. The two-dimensional moving weighted average method refers to the moving weighted average method of the two dimensions of the time domain space and the value domain space, and the method will be described in detail below.
S104,针对任一参考天,计算观测站该参考天与其他邻近天的平滑坐标序列的互相关系数,选取该参考天以及与该参考天互相关系数最大的邻近天的平滑坐标序列。S104 , for any reference day, calculate the cross-correlation coefficient of the smoothed coordinate sequence of the reference day and other adjacent days of the observation station, and select the smoothed coordinate sequence of the reference day and the adjacent day with the largest cross-correlation coefficient with the reference day.
本发明实施例中,多路径在临近天同一时间段内具有强相关性,经过二维移动加权平均法平滑处理后的观测点平滑坐标序列仍然含有观测噪声,但多路径误差占主导地位,故而也具有很强的相关性。互相关系数ρab的大小反映了两个随机变量间的相关程度,相关系数|ρab|越接近1,则A和B越接近线性相关,即相关程度越大。相关系数的计算公式可表示为:In the embodiment of the present invention, the multipath has strong correlation in the same time period of the adjacent days, and the smoothed coordinate sequence of observation points after smoothing by the two-dimensional moving weighted average method still contains observation noise, but the multipath error is dominant, so also has a strong correlation. The size of the cross-correlation coefficient ρ ab reflects the degree of correlation between two random variables. The closer the correlation coefficient |ρ ab | is to 1, the closer A and B are to a linear correlation, that is, the greater the degree of correlation. The formula for calculating the correlation coefficient can be expressed as:
其中,A和B为随机变量,ρab为变量A和B的互相关系数, 为变量A和B的协方差,为变量A的方差,为变量B的方差,为变A的平均值,为变量B的平均值。Among them, A and B are random variables, ρ ab is the cross-correlation coefficient of variables A and B, is the covariance of variables A and B, is the variance of variable A, is the variance of variable B, is the average value of variable A, is the mean of variable B.
本发明实施例中,A和B分别为参考天的平滑坐标序列和其他邻近天的平滑坐标序列。但由于受到残余观测噪声、卫星可视条件以及观测站天线与反射面之间的距离变化等因素的影响,使得临近天之间的平滑坐标序列的互相关系数变小。通过对观测点各天之间平滑坐标序列进行相关性分析,选取相关性最好的临近天的平滑坐标序列进行多路径误差的提取。In the embodiment of the present invention, A and B are respectively the smooth coordinate sequence of the reference day and the smooth coordinate sequence of other adjacent days. However, due to the influence of residual observation noise, satellite visibility conditions, and changes in the distance between the observation station antenna and the reflector, the cross-correlation coefficient of the smooth coordinate sequence between adjacent days becomes smaller. Through the correlation analysis of the smooth coordinate sequence between the observation points, the smooth coordinate sequence of the adjacent days with the best correlation is selected to extract the multi-path error.
S105,通过小波包算法对该参考天以及与该参考天互相关系数最大的邻近天的平滑坐标序列进行M层分解,计算各分解层中低频部分的互相关系数,M为预先设置的分解层数。S105, perform M-level decomposition on the smooth coordinate sequence of the reference day and the adjacent day with the largest cross-correlation coefficient with the reference day by using the wavelet packet algorithm, and calculate the cross-correlation coefficient of the low-frequency part in each decomposition layer, where M is a preset decomposition layer number.
其中,观测站的平滑坐标序列是离散序列,现有技术中多采用Mallat塔式算法进行离散二进小波变换来进行多路径误差的提取。但由于其尺度是按二进制变化,在高频频段其频率分辨率较差,而在低频段其时间分辨率较差,故而不能很好的分析观测站平滑坐标序列高频部分的多路径误差。本发明实施例通过小波包算法对观测站平滑坐标序列中的多路径误差在二维移动加权平均法的基础上进一步进行提取。该方法能将频带进行多层次划分,对二进小波变换中没有细分的高频部分进一步降半划分,并能够根据被分析信号的特征,使之与信号频谱相匹配。基于此特性,观测站平滑坐标序列中的高频部分的多路径误差将得到有效提取,将进一步削弱多路径效应在使用高精度定位进行变形观测中的影响。Among them, the smooth coordinate sequence of the observation station is a discrete sequence. In the prior art, the Mallat tower algorithm is often used to perform discrete binary wavelet transform to extract multi-path errors. However, because its scale changes in binary, its frequency resolution is poor in the high frequency band, and its time resolution is poor in the low frequency band, so the multipath error of the high frequency part of the smooth coordinate sequence of the observation station cannot be well analyzed. The embodiment of the present invention further extracts the multi-path error in the smooth coordinate sequence of the observation station on the basis of the two-dimensional moving weighted average method by using the wavelet packet algorithm. The method can divide the frequency band at multiple levels, further divide the high frequency part that is not subdivided in the binary wavelet transform into half, and can match the signal spectrum according to the characteristics of the analyzed signal. Based on this feature, the multi-path errors of the high-frequency part of the smooth coordinate sequence of the observation station will be effectively extracted, which will further weaken the influence of the multi-path effect in the deformation observation using high-precision positioning.
预先设置的分解层数M关系着多路径误差提取的有效性,分解层数过多,将会丢失平滑坐标序列中的重要信息,分解层数过少,则不能将平滑坐标序列中的多路径误差完全提取。分解层数M与平滑坐标序列的长度I之间满足关系:I=2M+l,l为余数项。The preset number of decomposition layers M is related to the effectiveness of multi-path error extraction. If the number of decomposition layers is too large, important information in the smooth coordinate sequence will be lost. Errors are fully extracted. The relationship between the number of decomposition levels M and the length I of the smooth coordinate sequence is satisfied: I=2 M +l, where l is the remainder item.
参见图2,图2为本发明实施例中3层分解后的小波包分解树结构图,其中,A表示低频部分,D表示高频部分,末尾的序号表示小波分解的层数。分解关系为S=AAA3+DAA3+ADA3+DDA3+AAD3+DAD3+ADD3+DDD3。Referring to FIG. 2, FIG. 2 is a structure diagram of a wavelet packet decomposition tree after 3-layer decomposition in an embodiment of the present invention, wherein A represents the low-frequency part, D represents the high-frequency part, and the sequence number at the end represents the number of layers of wavelet decomposition. The decomposition relationship is S=AAA3+DAA3+ADA3+DDA3+AAD3+DAD3+ADD3+DDD3.
S106,判断目标层J是否小于M,目标层J为各分解层中低频部分的互相关系数中最大值所在的层数。S106, determine whether the target layer J is less than M, and the target layer J is the layer number where the maximum value of the cross-correlation coefficients of the low-frequency parts in each decomposition layer is located.
本发明可以根据所选临近天对应低频部分的互相关系数的大小来确定最优分解层,即目标层,当确定对应的分解层低频部分的互相关系数达到最大时,该分解层为提取多路径误差的目标层。互相关系数越大表明多路径误差的占比在该分解层越大,但是低频部分的互相关系数随着分解层数的增加会表现出一条凸抛物线的变化趋势(分解层数一开始增加会更多地滤除噪声,但继续增大,原本是多路径误差的部分也被滤除了),也就是说,找到最高点即找到了最大互相关系数。假设预先设置的分解层数为M,若在M-1层找到了最大互相关系数,则可确定该点为最最大点,但如果最大相关系数出现在第M层,那么不能确定是不是最大点,只有再分解一层进行比较,如此循环反复就能确定最大互相关系数及需要分解的层数。因此,判断目标层J是否小于M,如果是,执行S107,否则,执行S108。The present invention can determine the optimal decomposition layer, that is, the target layer, according to the size of the cross-correlation coefficient of the low-frequency part corresponding to the selected adjacent days. Path error of the target layer. The larger the cross-correlation coefficient is, the larger the proportion of multi-path errors is in the decomposition layer, but the cross-correlation coefficient of the low-frequency part will show a convex parabolic change trend with the increase of the decomposition layer number (the decomposition layer will increase at the beginning. The noise is filtered out more, but it continues to increase, and the original part of the multi-path error is also filtered out), that is to say, the maximum cross-correlation coefficient is found when the highest point is found. Assuming that the preset number of decomposition layers is M, if the maximum cross-correlation coefficient is found in the M-1 layer, it can be determined that this point is the largest point, but if the largest correlation coefficient appears in the Mth layer, it is uncertain whether it is the largest point or not. point, only one more layer is decomposed for comparison, so the maximum cross-correlation coefficient and the number of layers to be decomposed can be determined by repeating the cycle. Therefore, it is judged whether the target layer J is smaller than M, and if so, S107 is executed, otherwise, S108 is executed.
S107,提取目标层J中的低频部分进行小波重构,得到多路径误差。S107, extracting the low-frequency part in the target layer J for wavelet reconstruction to obtain a multi-path error.
本步骤中,在确定目标层J后,将目标层J中包含随机噪声等影响的高频部分的小波包系数置零,对目标层J中的低频部分进行小波重构,即得到多路径误差。In this step, after the target layer J is determined, the wavelet packet coefficients of the high-frequency parts in the target layer J that contain random noise and other effects are set to zero, and the low-frequency part in the target layer J is reconstructed by wavelet, that is, the multi-path error is obtained. .
S108,将M的值加1,返回S105,直至目标层J小于M。S108, add 1 to the value of M, and return to S105 until the target layer J is less than M.
本发明实施例提供的误差提取方法,分别获取基准站和观测站在各天预设时间段内的接收信号,针对基准站和观测站每天的接收信号,对基准站和观测站在预设时间段内各历元的接收信号的伪距及载波相位进行差分处理,并根据得到的各天的双差伪距、双差载波相位以及基准站的位置坐标进行位置解算,得到观测站在各天的差分坐标序列;通过二维移动加权平均法对差分坐标序列进行处理,得到观测站在各天的平滑坐标序列,二维移动加权平均法基于时域空间和值域空间进行移动加权平均;通过小波包算法对互相关系数最大的参考天以及邻近天的平滑坐标序列进行M层分解,计算各分解层中低频部分的互相关系数,判断目标层J是否小于M,目标层J为各分解层中低频部分的互相关系数中最大值所在的层数;如果是,提取目标层J中的低频部分进行小波重构,得到多路径误差;如果否,将M的值加1,返回通过小波包算法对互相关系数最大的参考天以及邻近天的平滑坐标序列进行M层分解的步骤,直至目标层J小于M;提取目标层J中的低频部分进行小波重构,得到多路径误差。本发明实施例中,通过二维移动加权平均法可以削弱低频部分观测噪声对信号的影响,并解决信号平滑处理中的端部效应,达到保边去噪的目的。小波包算法可以将信号频带进行多层次划分,对高频部分进一步分解,消除高频部分中的观测噪声,从而提高多路径误差提取的准确性。In the error extraction method provided by the embodiment of the present invention, the received signals of the reference station and the observation station within the preset time period of each day are respectively obtained, and for the daily received signals of the reference station and the observation station, the preset time of the reference station and the observation station is obtained. The pseudorange and carrier phase of the received signal of each epoch in the segment are subjected to differential processing, and the position calculation is carried out according to the obtained double-difference pseudorange, double-difference carrier phase and the position coordinates of the reference station, and the observation station is obtained. The difference coordinate sequence of the day; the difference coordinate sequence is processed by the two-dimensional moving weighted average method to obtain the smooth coordinate sequence of the observation station on each day, and the two-dimensional moving weighted average method is based on the time domain space and the value domain space to perform moving weighted average; The wavelet packet algorithm is used to decompose the smooth coordinate sequence of the reference day with the largest cross-correlation coefficient and the adjacent days in M layers, calculate the cross-correlation coefficient of the low-frequency part in each decomposition layer, and judge whether the target layer J is less than M, and the target layer J is the decomposition of each layer. The number of layers where the maximum value of the cross-correlation coefficient of the low-frequency part in the layer is located; if yes, extract the low-frequency part in the target layer J for wavelet reconstruction to obtain the multi-path error; if not, add 1 to the value of M and return to the The packet algorithm performs M-level decomposition on the smooth coordinate sequence of the reference day with the largest cross-correlation coefficient and the adjacent days until the target layer J is less than M; extracts the low-frequency part of the target layer J for wavelet reconstruction, and obtains the multipath error. In the embodiment of the present invention, the influence of low-frequency part observation noise on the signal can be weakened by the two-dimensional moving weighted average method, and the end effect in the signal smoothing process can be solved, so as to achieve the purpose of edge preservation and denoising. The wavelet packet algorithm can divide the signal frequency band into multiple layers, further decompose the high-frequency part, and eliminate the observation noise in the high-frequency part, thereby improving the accuracy of multi-path error extraction.
本发明的一种实现方式中,图1实施例S103中,通过二维移动加权平均法对观测站在各天的差分坐标序列进行处理,得到观测站在各天的平滑坐标序列,包括以下步骤:In an implementation manner of the present invention, in the embodiment S103 of FIG. 1 , the differential coordinate sequence of the observation station on each day is processed by the two-dimensional moving weighted average method to obtain the smooth coordinate sequence of the observation station on each day, including the following steps :
第一步,若样本历元坐标的个数为I,确定时域空间的移动平均周期T为不大于I的奇数。In the first step, if the number of sample epoch coordinates is I, determine that the moving average period T of the time domain space is an odd number not greater than I.
其中,时域空间的移动平均周期T与多路径效应的影响周期一致,因此取为不大于I的奇数。Among them, the moving average period T of the time domain space is consistent with the influence period of the multipath effect, so it is taken as an odd number not greater than I.
第二步,若第i个样本历元坐标为xi, 表示值域空间,i为1~I的整数,根据公式:计算各样本历元坐标的平均值 In the second step, if the coordinate of the ith sample epoch is x i , Represents the value range space, i is an integer from 1 to I, according to the formula: Calculate the mean of the coordinates of each sample epoch
根据公式:计算各样本历元坐标的标准差σI;According to the formula: Calculate the standard deviation σ I of the coordinates of each sample epoch;
将标准差σI作为选取值域空间历元坐标的阈值θI,则θI=σI,Taking the standard deviation σ I as the threshold θ I for selecting the spatial epoch coordinates of the range, then θ I =σ I ,
第三步,若对第r个样本历元坐标进行平滑处理,r为1~I的整数;In the third step, if smoothing is performed on the coordinates of the rth sample epoch, r is an integer from 1 to I;
根据公式:计算值域空间内各样本历元坐标与当前样本历元坐标xr的标准差σr,According to the formula: Calculate the standard deviation σ r of each sample epoch coordinate and the current sample epoch coordinate x r in the range space,
将标准差σr不大于阈值θI的样本历元坐标作为分权历元坐标xk,分权历元坐标xk用于确定在值域空间第r个样本历元坐标平滑后的值,分权历元坐标xk的个数为第r个样本历元坐标在值域空间的移动平均周期Nr,Nr≤I。The sample epoch coordinates whose standard deviation σ r is not greater than the threshold θ I are taken as the decentralized epoch coordinates x k , and the decentralized epoch coordinates x k are used to determine the smoothed value of the rth sample epoch coordinates in the value domain space, The number of decentralized epoch coordinates x k is the moving average period N r of the r-th sample epoch coordinates in the value domain space, where N r ≤I.
本发明实施例中,值域空间的移动平均周期按照样本历元坐标的离散分布情况进行取值。第r个样本历元坐标在值域空间的移动平均周期Nr,也就是说,不同的样本历元坐标,对应不同的值域空间的移动平均周期。In the embodiment of the present invention, the moving average period of the range space is valued according to the discrete distribution of the coordinates of the sample epoch. The moving average period N r of the rth sample epoch coordinates in the value range space, that is, different sample epoch coordinates correspond to the moving average periods of different value range spaces.
第四步,对于第r个样本历元坐标,若αj表示时域空间第j个延迟历元坐标的加权系数,j为的整数,[]表示取整数, The fourth step, for the r-th sample epoch coordinate, if α j represents the weighting coefficient of the j-th delay epoch coordinate in the time domain space, j is Integer of , [] means to take an integer,
βk表示值域空间第k个分权历元坐标的加权系数,k为1~Nr的整数, β k represents the weighting coefficient of the k-th decentralized epoch coordinate in the range space, k is an integer from 1 to N r ,
根据公式:计算时域空间的移动平均周期T内各样本历元坐标相对于第r个样本历元坐标的平均延迟历元坐标为Δτj为第j个样本历元坐标相对于第r个样本历元坐标xr的延迟历元,According to the formula: Calculate the average delay epoch coordinates of each sample epoch coordinate relative to the rth sample epoch coordinate within the moving average period T of the time domain space as Δτ j is the delay epoch of the jth sample epoch coordinate relative to the rth sample epoch coordinate x r ,
根据公式:计算时域空间各延迟历元坐标的异化系数为 根据公式:计算时域空间第j个延迟历元坐标的加权系数αj。According to the formula: The dissimilation coefficient of each delay epoch coordinate in the time domain space is calculated as According to the formula: Calculate the weighting coefficient α j of the coordinate of the jth delay epoch in the time domain space.
在时域空间,计算样本历元坐标与第r个样本历元坐标的历元延迟,历元延迟越小,相应历元的值的权重越大。In the time domain space, the epoch delay between the sample epoch coordinates and the r-th sample epoch coordinates is calculated. The smaller the epoch delay, the greater the weight of the value of the corresponding epoch.
第五步,根据公式:计算值域空间内分权历元坐标xk的平均值Δμk为第k个分权历元坐标相对于第r个样本历元坐标的差值;The fifth step, according to the formula: Calculate the mean of the weighted epoch coordinates x k in the range space Δμ k is the difference between the coordinate of the k-th decentralized epoch relative to the coordinate of the r-th sample epoch;
根据公式:计算值域空间第k个分权历元坐标的异化系数 According to the formula: Calculate the dissimilation coefficient of the k-th weighted epoch coordinate in the range space
根据公式:计算值域空间第k个分权历元坐标的加权系数βk。According to the formula: Calculate the weighting coefficient β k of the coordinate of the k-th decentralized epoch in the range space.
在值域空间,计算所有分权历元坐标与第r个样本历元坐标值的差值,按照差值的大小对相应各个分权历元坐标进行权重分配,差值越小,说明相似度越高,相应分权历元坐标的权重也越大。In the value domain space, calculate the difference between the coordinates of all decentralized epochs and the r-th sample epoch coordinates, and assign weights to the corresponding decentralized epoch coordinates according to the size of the difference. The smaller the difference, the greater the similarity. The higher the weight, the greater the weight of the corresponding decentralized epoch coordinates.
第六步,根据公式:The sixth step, according to the formula:
将第r个样本历元坐标xr分别在时域空间和值域空间以xr为中心进行移动加权平均处理,得到第r个样本历元坐标的时域空间移动加权平均值以及值域空间移动平均加权值 Perform the moving weighted average processing of the rth sample epoch coordinate x r in the time domain space and the value domain space with x r as the center, and obtain the time domain space moving weighted average of the rth sample epoch coordinates and the range space moving average weighted value
第七步,根据公式:计算xr在二维空间下的平均值 The seventh step, according to the formula: Calculate the mean of x r in 2D space
根据公式:计算时域空间的异化系数 According to the formula: Calculate the dissimilation coefficient in the time-domain space
根据公式:计算值域空间的异化系数 According to the formula: Calculate the dissimilation coefficient of the range space
根据公式:计算时域空间的加权系数γT,根据公式:计算值域空间的加权系数 According to the formula: Calculate the weighting coefficient γ T of the time domain space, according to the formula: Calculate the weighting factor of the range space
根据公式:将和进行加权处理,得到经过平滑处理后的坐标 According to the formula: Will and Perform weighting processing to obtain smoothed coordinates
通过上述二维移动加权平均法可以解决平滑处理中的端部效应,达到保边去噪的目的。Through the above-mentioned two-dimensional moving weighted average method, the end effect in the smoothing process can be solved, and the purpose of edge-preserving and denoising can be achieved.
本发明的一种实现方式中,若该参考天以及与该参考天互相关系数最大的邻近天的平滑坐标序列为xi,根据公式:In an implementation manner of the present invention, if the smooth coordinate sequence of the reference day and the adjacent day with the largest cross-correlation coefficient with the reference day is x i , according to the formula:
对xi进行小波包分解, Perform wavelet packet decomposition on x i ,
i为1~I的整数,ω为1~2M的整数,q∈I,q为参考历元;i is an integer from 1 to I, ω is an integer from 1 to 2M, q∈I, q is a reference epoch;
表示未分解时的小波包,表示第M层上的第2ω-1个小波包系数,表示第M-1层上的第ω个小波包系数,表示第M层上的第2ω个小波包系数,G为尺度函数分解滤波器,H为小波函数分解滤波器。 represents the wavelet packet when it is not decomposed, represents the 2ω-1 wavelet packet coefficient on the Mth layer, represents the ωth wavelet packet coefficient on the M-1th layer, Represents the second ω wavelet packet coefficient on the Mth layer, G is the scale function decomposition filter, and H is the wavelet function decomposition filter.
本发明实施例中,在进行小波包分解时,由于子空间有不同的分解方式,即子空间有不同的正交基。在理论分析中,通常可以定义一个序列的代价函数,从小波库的所有小波基中寻找使代价函数最小的基,对一个给定向量来说,代价最小就是最有效的表示,该小波基即为“最优基”。在实际应用中,一般根据信号分析的目的和经验来选择小波基。在GPS数据处理方面在理论上提出利用Symlets小波,但是在实际应用中大多采用Daubechies小波。In the embodiment of the present invention, when the wavelet packet is decomposed, because the subspaces have different decomposition modes, that is, the subspaces have different orthogonal bases. In theoretical analysis, a cost function of a sequence can usually be defined, and the basis that minimizes the cost function can be found from all wavelet bases in the wavelet library. For a given vector, the least cost is the most effective representation, and the wavelet basis is is the "optimal basis". In practical applications, the wavelet basis is generally selected according to the purpose and experience of signal analysis. In the aspect of GPS data processing, Symlets wavelet is proposed in theory, but Daubechies wavelet is mostly used in practical applications.
提取目标层J中的低频部分进行小波重构,得到多路径误差,包括:Extract the low-frequency part in the target layer J for wavelet reconstruction, and obtain the multi-path error, including:
根据公式:对xi进行小波包重构;According to the formula: Perform wavelet packet reconstruction on xi ;
表示第J+1层上的第2ω-1个小波包系数,表示第J+1层上的第2ω个小波包系数,g为尺度函数重构滤波器,h为小波函数重构滤波器。 represents the 2ω-1 wavelet packet coefficient on the J+1th layer, represents the second ω wavelet packet coefficient on the J+1th layer, g is the scale function reconstruction filter, and h is the wavelet function reconstruction filter.
具体的,将包含随机噪声等影响的高频部分的小波包系数置零,即得到坐标残差序列中的多路径误差。Specifically, the wavelet packet coefficients of the high-frequency part including random noise and the like are set to zero, that is, the multi-path error in the coordinate residual sequence is obtained.
实施例一Example 1
若模拟信号模型为: If the analog signal model is:
s由一个周期为1200s的余弦信号和两个周期分别为900s和300s的正弦信号以及一个高斯白噪声序列e组成。取e的噪声水平为N(0,1.52),数据采样率设为1s,采样个数为3000。s consists of a cosine signal with a period of 1200s, two sinusoidal signals with periods of 900s and 300s, and a Gaussian white noise sequence e. Take the noise level of e as N(0,1.52), the data sampling rate as 1s, and the number of samples as 3000.
将该模拟信号分别通过小波阈值去噪算法、二维移动加权平均算法以及基于二维移动加权平均处理的小波包算法进行滤波去噪,小波函数均选用db8,分解层数为5层,去噪信号如图3所示,残差信号如图4所示。The simulated signal is filtered and denoised by wavelet threshold denoising algorithm, two-dimensional moving weighted average algorithm and wavelet packet algorithm based on two-dimensional moving weighted average processing. The wavelet functions are all selected from db8, and the number of decomposition layers is 5 layers. The signal is shown in Figure 3, and the residual signal is shown in Figure 4.
由图3和图4可以看出,本发明实施例的基于二维移动加权平均处理的小波包算法去噪后的信号能更好地还原原始信号,而经二维移动加权平均算法处理过的信号还原性稍好于小波阈值去噪算法。It can be seen from FIG. 3 and FIG. 4 that the signal denoised by the wavelet packet algorithm based on the two-dimensional moving weighted average processing according to the embodiment of the present invention can better restore the original signal, while the signal processed by the two-dimensional moving weighted average algorithm can better restore the original signal. The signal restoration is slightly better than the wavelet threshold denoising algorithm.
将以上三种算法去噪后还原的信号和产生的残差信号进行定量分析,分析指标包括去噪后信号中的SRMS以及去噪后的信号与原始信号的相关系数ρ。为了科学评价三种算法的性能,分别在4种不同的高斯白噪声模拟环境下进行了仿真分析,结果如表1、表2和表3所示。表1为不同高斯白噪声模拟环境下小波阈值去噪的性能指标,表2为不同高斯白噪声模拟环境下二维移动加权平均法去噪的性能指标,表3为不同高斯白噪声模拟环境下基于二维移动加权平均处理的小波包去噪的性能指标。Quantitative analysis is carried out on the restored signal and the generated residual signal after denoising by the above three algorithms. The analysis indicators include S RMS in the denoised signal and the correlation coefficient ρ between the denoised signal and the original signal. In order to scientifically evaluate the performance of the three algorithms, simulation analysis was carried out under four different Gaussian white noise simulation environments, and the results are shown in Table 1, Table 2 and Table 3. Table 1 shows the performance indexes of wavelet threshold denoising under different Gaussian white noise simulation environments, Table 2 shows the performance indexes of two-dimensional moving weighted average denoising under different Gaussian white noise simulation environments, and Table 3 shows the performance indexes under different Gaussian white noise simulation environments Performance metrics for wavelet packet denoising based on two-dimensional moving weighted average processing.
表1Table 1
表2Table 2
表3table 3
由表1、表2和表3可知,上述三种算法的相关系数在不同噪声水平下指标大部分都超过了0.95,说明它们都能较好的还原原始信号的信息。但是随着噪声水平的不断增大,小波阈值去噪算法和二维移动加权平均处理算法的SRMS明显增加,但小波阈值去噪在高噪声水平情况下增加的较为缓慢,而二维移动加权平均处理算法相对较快,但其SRMS在高噪声水平下任然优于小波阈值去噪算法。It can be seen from Table 1, Table 2 and Table 3 that the correlation coefficients of the above three algorithms mostly exceed 0.95 under different noise levels, indicating that they can restore the original signal information well. However, with the increasing noise level, the S RMS of the wavelet threshold denoising algorithm and the two-dimensional moving weighted average processing algorithm increase significantly, but the wavelet threshold denoising increases slowly in the case of high noise level, while the two-dimensional moving weighted The averaging processing algorithm is relatively fast, but its S RMS is still better than the wavelet threshold denoising algorithm at high noise levels.
基于二维移动加权平均处理的小波包算法的SRMS在不同噪声水平下都相对最优,且其相关系数都大于0.98。The S RMS of the wavelet packet algorithm based on two-dimensional moving weighted average processing is relatively optimal under different noise levels, and its correlation coefficients are all greater than 0.98.
综上所述,基于二维移动加权平均处理的小波包算法可以更好地对加噪后的信号进行滤波去噪,表明本发明可以精确地提取多路径误差。To sum up, the wavelet packet algorithm based on the two-dimensional moving weighted average processing can better filter and denoise the noised signal, indicating that the present invention can accurately extract multipath errors.
相应于上述方法实施例,本发明实施例提供了一种误差提取装置,参见图5,图5为本发明实施例的误差提取装置的结构图,包括:Corresponding to the above method embodiments, an embodiment of the present invention provides an error extraction apparatus. Referring to FIG. 5, FIG. 5 is a structural diagram of the error extraction apparatus according to the embodiment of the present invention, including:
双差计算模块501,用于分别获取基准站和观测站在各天预设时间段内的接收信号,针对基准站和观测站每天的接收信号,对基准站和观测站在预设时间段内各历元的接收信号的伪距及载波相位进行差分处理,得到双差伪距和双差载波相位;The double-
差分坐标序列计算模块502,用于根据基准站的位置坐标、各天的双差伪距和双差载波相位进行位置解算,得到观测站在各天的差分坐标序列;The differential coordinate
平滑处理模块503,用于通过二维移动加权平均法对观测站在各天的差分坐标序列进行处理,得到观测站在各天的平滑坐标序列,二维移动加权平均法基于时域空间和值域空间进行移动加权平均;The smoothing
平滑坐标序列选取模块504,用于针对任一参考天,计算观测站该参考天与其他邻近天的平滑坐标序列的互相关系数,选取该参考天以及与该参考天互相关系数最大的邻近天的平滑坐标序列;The smooth coordinate
分解模块505,用于通过小波包算法对该参考天以及与该参考天互相关系数最大的邻近天的平滑坐标序列进行M层分解,计算各分解层中低频部分的互相关系数,M为预先设置的分解层数;The
判断模块506,用于判断目标层J是否小于M,目标层J为各分解层中低频部分的互相关系数中最大值所在的层数;The
重构模块507,用于在判断模块的判断结果为是时,提取目标层J中的低频部分进行小波重构,得到多路径误差;The
循环模块508,用于在判断模块的判断结果为否时,将M的值加1,返回通过小波包算法对该参考天以及与该参考天互相关系数最大的邻近天的平滑坐标序列进行M层分解的步骤,直至目标层J小于M。The
本发明实施例提供的误差提取装置,分别获取基准站和观测站在各天预设时间段内的接收信号,针对基准站和观测站每天的接收信号,对基准站和观测站在预设时间段内各历元的接收信号的伪距及载波相位进行差分处理,并根据得到的各天的双差伪距、双差载波相位以及基准站的位置坐标进行位置解算,得到观测站在各天的差分坐标序列;通过二维移动加权平均法对差分坐标序列进行处理,得到观测站在各天的平滑坐标序列,二维移动加权平均法基于时域空间和值域空间进行移动加权平均;通过小波包算法对互相关系数最大的参考天以及邻近天的平滑坐标序列进行M层分解,计算各分解层中低频部分的互相关系数,判断目标层J是否小于M,目标层J为各分解层中低频部分的互相关系数中最大值所在的层数;如果是,提取目标层J中的低频部分进行小波重构,得到多路径误差;如果否,将M的值加1,返回通过小波包算法对互相关系数最大的参考天以及邻近天的平滑坐标序列进行M层分解的步骤,直至目标层J小于M;提取目标层J中的低频部分进行小波重构,得到多路径误差。本发明实施例中,通过二维移动加权平均法可以削弱低频部分观测噪声对信号的影响,并解决信号平滑处理中的端部效应,达到保边去噪的目的。小波包算法可以将信号频带进行多层次划分,对高频部分进一步分解,消除高频部分中的观测噪声,从而提高多路径误差提取的准确性。The error extraction device provided by the embodiment of the present invention separately obtains the received signals of the reference station and the observation station within the preset time period of each day, and for the daily received signals of the reference station and the observation station, the preset time of the reference station and the observation station is obtained. The pseudorange and carrier phase of the received signal of each epoch in the segment are subjected to differential processing, and the position calculation is carried out according to the obtained double-difference pseudorange, double-difference carrier phase and the position coordinates of the reference station, and the observation station is obtained. The difference coordinate sequence of the day; the difference coordinate sequence is processed by the two-dimensional moving weighted average method to obtain the smooth coordinate sequence of the observation station on each day, and the two-dimensional moving weighted average method is based on the time domain space and the value domain space to perform moving weighted average; The wavelet packet algorithm is used to decompose the smooth coordinate sequence of the reference day with the largest cross-correlation coefficient and the adjacent days in M layers, calculate the cross-correlation coefficient of the low-frequency part in each decomposition layer, and judge whether the target layer J is less than M, and the target layer J is the decomposition of each layer. The number of layers where the maximum value of the cross-correlation coefficient of the low-frequency part in the layer is located; if yes, extract the low-frequency part in the target layer J for wavelet reconstruction to obtain the multi-path error; if not, add 1 to the value of M and return to the The packet algorithm performs M-level decomposition on the smooth coordinate sequence of the reference day with the largest cross-correlation coefficient and the adjacent days until the target layer J is less than M; extracts the low-frequency part of the target layer J for wavelet reconstruction, and obtains the multipath error. In the embodiment of the present invention, the influence of low-frequency part observation noise on the signal can be weakened by the two-dimensional moving weighted average method, and the end effect in the signal smoothing process can be solved, so as to achieve the purpose of edge preservation and denoising. The wavelet packet algorithm can divide the signal frequency band into multiple layers, further decompose the high-frequency part, and eliminate the observation noise in the high-frequency part, thereby improving the accuracy of multi-path error extraction.
可选的,平滑处理模块具体用于,若样本历元坐标的个数为I,确定时域空间的移动平均周期T为不大于I的奇数;Optionally, the smoothing processing module is specifically used to, if the number of sample epoch coordinates is 1, determine that the moving average period T of the time domain space is an odd number not greater than 1;
若第i个样本历元坐标为xi, 表示值域空间,i为1~I的整数,根据公式:计算各样本历元坐标的平均值 If the coordinate of the ith sample epoch is x i , Represents the value range space, i is an integer from 1 to I, according to the formula: Calculate the mean of the coordinates of each sample epoch
根据公式:计算各样本历元坐标的标准差σI;According to the formula: Calculate the standard deviation σ I of the coordinates of each sample epoch;
将标准差σI作为选取值域空间历元坐标的阈值θI,则θI=σI,Taking the standard deviation σ I as the threshold θ I for selecting the spatial epoch coordinates of the range, then θ I =σ I ,
若对第r个样本历元坐标进行平滑处理,r为1~I的整数;If smoothing is performed on the coordinates of the rth sample epoch, r is an integer from 1 to I;
根据公式:计算值域空间内各样本历元坐标与当前样本历元坐标xr的标准差σr,According to the formula: Calculate the standard deviation σ r of each sample epoch coordinate and the current sample epoch coordinate x r in the range space,
将标准差σr不大于阈值θI的样本历元坐标作为分权历元坐标xk,分权历元坐标xk用于确定在值域空间第r个样本历元坐标平滑后的值,分权历元坐标xk的个数为第r个样本历元坐标在值域空间的移动平均周期Nr,Nr≤I;对于第r个样本历元坐标,若αj表示时域空间第j个延迟历元坐标的加权系数,j为的整数,[]表示取整数, The sample epoch coordinates whose standard deviation σ r is not greater than the threshold θ I are taken as the decentralized epoch coordinates x k , and the decentralized epoch coordinates x k are used to determine the smoothed value of the rth sample epoch coordinates in the value domain space, The number of decentralized epoch coordinates x k is the moving average period N r of the r-th sample epoch coordinates in the value domain space, N r ≤I; for the r-th sample epoch coordinates, if α j represents the time domain space Weighting coefficient of the jth delay epoch coordinate, j is Integer of , [] means to take an integer,
βk表示值域空间第k个分权历元坐标的加权系数,k为1~Nr的整数, β k represents the weighting coefficient of the k-th decentralized epoch coordinate in the range space, k is an integer from 1 to N r ,
根据公式:计算时域空间的移动平均周期T内各样本历元坐标相对于第r个样本历元坐标的平均延迟历元坐标为Δτj为第j个样本历元坐标相对于第r个样本历元坐标xr的延迟历元,According to the formula: Calculate the average delay epoch coordinates of each sample epoch coordinate relative to the rth sample epoch coordinate within the moving average period T of the time domain space as Δτ j is the delay epoch of the jth sample epoch coordinate relative to the rth sample epoch coordinate x r ,
根据公式:计算时域空间各延迟历元坐标的异化系数为 根据公式:计算时域空间第j个延迟历元坐标的加权系数αj;According to the formula: The dissimilation coefficient of each delay epoch coordinate in the time domain space is calculated as According to the formula: Calculate the weighting coefficient α j of the coordinate of the jth delay epoch in the time domain space;
根据公式:计算值域空间内分权历元坐标xk的平均值Δμk为第k个分权历元坐标相对于第r个样本历元坐标的差值;According to the formula: Calculate the mean of the weighted epoch coordinates x k in the range space Δμ k is the difference between the coordinate of the k-th decentralized epoch relative to the coordinate of the r-th sample epoch;
根据公式:计算值域空间第k个分权历元坐标的异化系数 According to the formula: Calculate the dissimilation coefficient of the k-th weighted epoch coordinate in the range space
根据公式:计算值域空间第k个分权历元坐标的加权系数βk;According to the formula: Calculate the weighting coefficient β k of the coordinate of the k-th decentralized epoch in the value domain space;
根据公式:According to the formula:
将第r个样本历元坐标xr分别在时域空间和值域空间以xr为中心进行移动加权平均处理,得到第r个样本历元坐标的时域空间移动加权平均值以及值域空间移动平均加权值 Perform the moving weighted average processing of the rth sample epoch coordinate x r in the time domain space and the value domain space with x r as the center, and obtain the time domain space moving weighted average of the rth sample epoch coordinates and the range space moving average weighted value
根据公式:计算xr在二维空间下的平均值 According to the formula: Calculate the mean of x r in 2D space
根据公式:计算时域空间的异化系数 According to the formula: Calculate the dissimilation coefficient in the time-domain space
根据公式:计算值域空间的异化系数 According to the formula: Calculate the dissimilation coefficient of the range space
根据公式:计算时域空间的加权系数γT,根据公式:计算值域空间的加权系数 According to the formula: Calculate the weighting coefficient γ T of the time domain space, according to the formula: Calculate the weighting factor of the range space
根据公式:将和进行加权处理,得到经过平滑处理后的坐标 According to the formula: Will and Perform weighting processing to obtain smoothed coordinates
可选的,双差计算模块,具体用于根据公式:Optional, double difference calculation module, specifically used according to the formula:
计算观测站r和基准站b在同一历元观测卫星n和卫星m的双差伪距以及观测站r和基准站b在同一历元观测卫星n和卫星m的双差载波相位 Calculate the double-difference pseudoranges of satellite n and satellite m observed by observation station r and base station b in the same epoch and the double-difference carrier phase of satellite n and satellite m observed by observation station r and reference station b in the same epoch
其中,为双差算子,S表示全球定位系统GPS系统,卫星n和卫星m为该GPS系统内的卫星,表示观测站r和基准站b在同一历元观测卫星n和卫星m时卫星中心到接收机相位中心间的几何距离,λ为载波波长,为观测站r和基准站b在同一历元观测卫星n和卫星m的载波相位未知整周模糊度,表示基准站b和观测站r在同一历元观测卫星n和卫星m的伪距多路径误差,表示观测站r和基准站b在同一历元观测卫星n和卫星m的载波相位多路径误差,表示观测站r和基准站b在同一历元观测卫星n和卫星m的伪码,表示观测站r和基准站b在同一历元观测卫星n和卫星m的载波相位测量噪声。in, is a double difference operator, S represents the global positioning system GPS system, satellite n and satellite m are satellites in the GPS system, represents the geometric distance between the satellite center and the receiver phase center when the observation station r and the base station b observe satellite n and satellite m in the same epoch, λ is the carrier wavelength, is the unknown integer ambiguity of the carrier phase of satellite n and satellite m observed at the same epoch by observation station r and base station b, represents the pseudorange multipath error of satellite n and satellite m observed by base station b and observation station r in the same epoch, represents the carrier phase multipath error of satellite n and satellite m observed by observation station r and reference station b in the same epoch, Pseudocode indicating that observation station r and base station b observe satellite n and satellite m in the same epoch, It represents the carrier phase measurement noise of satellite n and satellite m observed by observation station r and reference station b in the same epoch.
可选的,分解模块,具体用于若该参考天以及与该参考天互相关系数最大的邻近天的平滑坐标序列为xi,根据公式:Optionally, the decomposition module is specifically used if the smooth coordinate sequence of the reference day and the adjacent day with the largest cross-correlation coefficient with the reference day is x i , according to the formula:
对xi进行小波包分解, Perform wavelet packet decomposition on x i ,
i为1~I的整数,ω为1~2M的整数,q∈I,q为参考历元;i is an integer from 1 to I, ω is an integer from 1 to 2M, q∈I, q is a reference epoch;
表示未分解时的小波包,表示第M层上的第2ω-1个小波包系数,表示第M-1层上的第ω个小波包系数,表示第M层上的第2ω个小波包系数,G为尺度函数分解滤波器,H为小波函数分解滤波器, represents the wavelet packet when it is not decomposed, represents the 2ω-1 wavelet packet coefficient on the Mth layer, represents the ωth wavelet packet coefficient on the M-1th layer, represents the second ω wavelet packet coefficient on the Mth layer, G is the scale function decomposition filter, H is the wavelet function decomposition filter,
重构模块,具体用于根据公式:Refactoring the module, specifically for use according to the formula:
对xi进行小波包重构;表示第J+1层上的第2ω-1个小波包系数,表示第J+1层上的第2ω个小波包系数,g为尺度函数重构滤波器,h为小波函数重构滤波器。 Perform wavelet packet reconstruction on xi ; represents the 2ω-1 wavelet packet coefficient on the J+1th layer, represents the second ω wavelet packet coefficient on the J+1th layer, g is the scale function reconstruction filter, and h is the wavelet function reconstruction filter.
本发明实施例还提供了一种电子设备,参见图6,图6为本发明实施例的电子设备的结构图,包括:处理器601、通信接口602、存储器603和通信总线604,其中,处理器601、通信接口602、存储器603通过通信总线604完成相互间的通信;An embodiment of the present invention also provides an electronic device. Referring to FIG. 6, FIG. 6 is a structural diagram of the electronic device according to the embodiment of the present invention, including: a
存储器603,用于存放计算机程序;a
处理器601,用于执行存储器603上所存放的程序时,实现上述任一误差提取方法的步骤。The
需要说明的是,上述电子设备提到的通信总线604可以是PCI(PeripheralComponent Interconnect,外设部件互连标准)总线或EISA(Extended Industry StandardArchitecture,扩展工业标准结构)总线等。该通信总线604可以分为地址总线、数据总线、控制总线等。为便于表示,图6中仅用一条粗线表示,但并不表示仅有一根总线或一种类型的总线。It should be noted that the
通信接口602用于上述电子设备与其他设备之间的通信。The
存储器603可以包括RAM(Random Access Memory,随机存取存储器),也可以包括非易失性存储器(non-volatile memory),例如至少一个磁盘存储器。可选的,存储器还可以是至少一个位于远离前述处理器的存储装置。The
上述的处理器601可以是通用处理器,包括:CPU(Central Processing Unit,中央处理器)、NP(Network Processor,网络处理器)等;还可以是DSP(Digital SignalProcessing,数字信号处理器)、ASIC(Application Specific Integrated Circuit,专用集成电路)、FPGA(Field-Programmable Gate Array,现场可编程门阵列)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件。The above-mentioned
本发明实施例还提供了一种计算机可读存储介质,计算机可读存储介质内存储有计算机程序,计算机程序被处理器执行时,实现上述任一误差提取方法的步骤。Embodiments of the present invention further provide a computer-readable storage medium, where a computer program is stored in the computer-readable storage medium, and when the computer program is executed by a processor, the steps of any of the above error extraction methods are implemented.
需要说明的是,在本文中,诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来,而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且,术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含,从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素,而且还包括没有明确列出的其他要素,或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下,由语句“包括一个……”限定的要素,并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。It should be noted that, in this document, relational terms such as first and second are only used to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any relationship between these entities or operations. any such actual relationship or sequence exists. Moreover, the terms "comprising", "comprising" or any other variation thereof are intended to encompass a non-exclusive inclusion such that a process, method, article or device that includes a list of elements includes not only those elements, but also includes not explicitly listed or other elements inherent to such a process, method, article or apparatus. Without further limitation, an element qualified by the phrase "comprising a..." does not preclude the presence of additional identical elements in a process, method, article or apparatus that includes the element.
本说明书中的各个实施例均采用相关的方式描述,各个实施例之间相同相似的部分互相参见即可,每个实施例重点说明的都是与其他实施例的不同之处。尤其,对于装置、电子设备及可读存储介质实施例而言,由于其基本相似于方法实施例,所以描述的比较简单,相关之处参见方法实施例的部分说明即可。Each embodiment in this specification is described in a related manner, and the same and similar parts between the various embodiments may be referred to each other, and each embodiment focuses on the differences from other embodiments. Especially, for the embodiments of the apparatus, electronic device and readable storage medium, since they are basically similar to the method embodiments, the description is relatively simple, and reference may be made to some descriptions of the method embodiments for related parts.
以上所述仅为本发明的较佳实施例而已,并非用于限定本发明的保护范围。凡在本发明的精神和原则之内所作的任何修改、等同替换、改进等,均包含在本发明的保护范围内。The above descriptions are only preferred embodiments of the present invention, and are not intended to limit the protection scope of the present invention. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention are included in the protection scope of the present invention.
Claims (10)
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910096259.3A CN109725290B (en) | 2019-01-31 | 2019-01-31 | Error extraction method and device, electronic equipment and readable storage medium |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910096259.3A CN109725290B (en) | 2019-01-31 | 2019-01-31 | Error extraction method and device, electronic equipment and readable storage medium |
Publications (2)
Publication Number | Publication Date |
---|---|
CN109725290A CN109725290A (en) | 2019-05-07 |
CN109725290B true CN109725290B (en) | 2020-10-13 |
Family
ID=66300438
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201910096259.3A Active CN109725290B (en) | 2019-01-31 | 2019-01-31 | Error extraction method and device, electronic equipment and readable storage medium |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN109725290B (en) |
Families Citing this family (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112180408B (en) * | 2020-09-29 | 2023-06-23 | 中山大学 | An intelligent terminal-based multipath error extraction method and related device |
CN113075706A (en) * | 2021-03-25 | 2021-07-06 | 上海海洋大学 | GNSS-R based snow depth inversion method and application thereof |
CN113189624B (en) * | 2021-04-30 | 2023-10-03 | 中山大学 | Self-adaptive classification multipath error extraction method and device |
CN115060160A (en) * | 2022-05-23 | 2022-09-16 | 首钢京唐钢铁联合有限责任公司 | Hot blast stove installation space precision positioning method, device and electronic equipment |
CN115752399A (en) * | 2022-11-15 | 2023-03-07 | 中国电建集团成都勘测设计研究院有限公司 | Annual adjustment reservoir vertical line displacement result correction method |
CN115808703B (en) * | 2022-12-14 | 2024-02-23 | 北京六分科技有限公司 | Multipath influence degree detection method, device, storage medium and program product |
Family Cites Families (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102594472A (en) * | 2012-03-22 | 2012-07-18 | 北京邮电大学 | Method and system for wireless channel measurement based on wavelet decomposition threshold de-nosing |
US9864064B2 (en) * | 2012-06-27 | 2018-01-09 | Mitsubishi Electric Corporation | Positioning device |
JP2017143459A (en) * | 2016-02-12 | 2017-08-17 | 国立研究開発法人産業技術総合研究所 | Method and device for measuring propagation delay characteristics |
CN106646538B (en) * | 2016-10-31 | 2019-06-04 | 东南大学 | A Multipath Correction Method for Deformation Monitoring GNSS Signals Based on Monodifference Filtering |
-
2019
- 2019-01-31 CN CN201910096259.3A patent/CN109725290B/en active Active
Also Published As
Publication number | Publication date |
---|---|
CN109725290A (en) | 2019-05-07 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN109725290B (en) | Error extraction method and device, electronic equipment and readable storage medium | |
Su et al. | A new multipath mitigation method based on adaptive thresholding wavelet denoising and double reference shift strategy | |
Lau | Wavelet packets based denoising method for measurement domain repeat-time multipath filtering in GPS static high-precision positioning | |
Jacobs et al. | The Murchison widefield array 21 cm power spectrum analysis methodology | |
Prasad et al. | FLAGCAL: a flagging and calibration package for radio interferometric data | |
van der Veen et al. | Signal processing tools for radio astronomy | |
Offringa et al. | Precision requirements for interferometric gridding in the analysis of a 21 cm power spectrum | |
Karaçaylı et al. | Optimal 1D Ly α forest power spectrum estimation–III. DESI early data | |
Shen et al. | Site-specific real-time GPS multipath mitigation based on coordinate time series window matching | |
Mao et al. | An Improved DOA Estimation Algorithm Based on Wavelet Operator. | |
Sun et al. | Improving the resolution of underwater acoustic image measurement by deconvolution | |
Zhang et al. | Parameterized CLEAN deconvolution in radio synthesis imaging | |
Ansari et al. | Wavelet and power spectrum analysis of global navigation satellite system multipath error modeling and mitigation | |
Cheng et al. | Comprehensive analysis of multipath estimation algorithms in the framework of information theoretic learning | |
CN112926190A (en) | Multi-path weakening method and device based on VMD algorithm | |
Carroll et al. | Detecting and quantifying stellar magnetic fields-Sparse Stokes profile approximation using orthogonal matching pursuit | |
CN110657742A (en) | Aquifer deformation signal separation method, device, device and readable storage medium | |
CN115795272A (en) | Spectral line denoising method based on fractional order iterative discrete wavelet transform | |
Bahr et al. | Wavespace-based coherent deconvolution | |
Leyang et al. | Sequential bias-corrected weighted least squares solution of mixed additive and multiplicative error models | |
CN119001660B (en) | A method, device and medium for deambiguating uniform circular array spatial spectrum direction finding | |
CN113658606B (en) | A Beamforming Method Based on Adaptive Compressive Sensing under Low SNR Conditions | |
CN111694025A (en) | Unambiguous multi-path suppression method suitable for MBOC navigation signal | |
CN118465684B (en) | A distributed two-dimensional nested array and its algorithm for obtaining target position | |
CN110554367B (en) | Target scattering characteristic measurement interference removal method based on compressed sensing |
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 |