CN109639612A - 一种基于非线性最小二乘法的zpw-2000信号解调方法 - Google Patents

Abstract

本发明公开了一种基于非线性最小二乘法的ZPW‑2000信号解调方法。该方法基于非线性最小二乘法法则，包含数据预处理、粗估计、网格搜索及精确搜索五个步骤。本发明具有较高的解调准确度。高信噪比环境下采样持续时间大于等于0.11s的ZPW‑2000移频信号均可被本发明正确解调。同时，本发明具备实时解调能力，较好的抗白噪声能力以及较好的抗单频干扰能力。

Description

一种基于非线性最小二乘法的ZPW-2000信号解调方法

技术领域

本发明属于铁路信号轨道电路领域，涉及一种基于非线性最小二乘法的 ZPW-2000信号解调方法，应用于普速及配备CTCS-2级列控系统的列车车载设备中ZPW-2000轨道电路移频信号的解调。

背景技术

ZPW-2000轨道电路作为一种地-车不间断通信设备，被广泛应用于我国既有线及客运专线，在ZPW-2000轨道电路区段，车载TCR(track circuit reader)天线不断感应以钢轨为传输媒介的ZPW-2000移频信号，车载设备利用感应电压采样值解调得到低频码序，从而获知前方区段占用情况，指导列车安全运行。

由于ZPW-2000移频信号(以下简称ZPW-2000信号)与列车牵引回流共用钢轨作为传输媒介，不平衡牵引电流干扰对ZPW-2000信号影响较为严重。同时，相关规范要求车载设备接收信息的应变时间不能超过2s。因此，ZPW-2000信号解调方法应兼备抗干扰能力和实时解调能力。

ZPW-2000信号为相位连续二进制频移键控(CPBFSK,continuous phasemodulated binary frequency shift keying)信号，其二进制符号”0”、”1”始终处于交替变化状态，使得其频谱形态基本恒定。考虑采样过程，ZPW-2000信号s(n)的傅立叶级数展开形式为：

式中，A₀、f₀、f₁分别为s(n)的包络幅值、载频频率和低频频率；m＝f₁/f_Δ为调制指数；f_Δ为频偏，恒为11Hz；T_s为采样间隔。式(1)描述了ZPW-2000信号理想频谱的形态，依据该理想频谱形态，发展出一系列基于信号频域的解调方法。

不平衡牵引电流干扰主要包括牵引谐波干扰和暂态冲击干扰。牵引谐波干扰谐波频率间隔为50Hz左右，因此在ZPW-2000信号有效频带[f₀-40，f₀+40](Hz) 内只可能存在一个谐波频率分量。因此本发明将牵引谐波干扰描述为单频干扰。暂态冲击干扰频率成分复杂，不失一般地，本发明将其假定为高斯白噪声。因此，在本发明实施例的仿真试验中，将抗白噪声能力和抗单频干扰能力作为解调方法抗干扰能力的评判指标。

现有机车信号设备通常采用基于FFT(fast Fourier transform)的周期图法(亦称频域分析法)来解调ZPW-2000信号，主要有：文献“移频轨道电路测试系统的设计”(魏学业,汪希时,丁正庭.铁道学报,18(05):67-72,1996)提出的结合欠采样、 FFT滤波及补零FFT的解调方法，文献“铁路移频信号处理方法研究”(胡幸江,黄文君,何伟挺,谭平.仪器仪表学报,33(08):1729-1734,2012)提出的结合欠采样和补零FFT的解调方法(以下简称“欠采样补零法”)，文献“铁路移频信号检测系统设计与实现”(杨帆,戴胜华,刘泽.电子测量与仪器学报,24(05):500-505,2010)提出的 Zoom-FFT结合重心法插值的解调方法(以下简称“Zoom-FFT插值法”)。由于 ZPW-2000信号各频率成分分布较集中，而周期图法存在频谱泄漏现象，能量较大的载频分量的泄漏值容易使附近能量较小的边频分量谱峰位置发生偏移，当采样持续时间较小时，这种偏移现象会更加严重，此时周期图法的解调准确度会大幅降低，解调正确性难以保障。当列车高速通过较短的道岔区段以及信号受到暂态冲击干扰影响时，有效信号的采样持续时间长度被缩短，周期图法的解调正确性会受到挑战。因此，通过采用新的解调方法，提升解调准确度，进而降低解调方法对采样持续时间长度的限制，对于改善ZPW-2000轨道电路这一地-车通信方式的可靠性具有重要意义。

发明内容

为解决现有技术存在的上述问题，本发明提出基于非线性最小二乘法的 ZPW-2000信号解调方法，其解调准确度优于目前常用的周期图法而且具备实时解调能力。

一种基于非线性最小二乘法的ZPW-2000信号解调方法，该方法的实现目的为寻求以下价值函数的最大值：

式中，分别为载频频率f₀、低频频f₁估计值；x＝[x(0)x(1)x(2)…x(N-1)]^T为接收信号采样序列；N为采样序列长度，序列下标从0开始；[·]^T为矩阵转置符号；Z为下文定义的ZPW-2000信号近似模型矩阵形式中的参数矩阵；

依据ZPW-2000信号的傅立叶级数展开形式，定义ZPW-2000信号s(n)的近似模型如下：

式中，L为近似模型中频率分量个数，本发明令L＝5；ω_l，l＝1,2,...,L为近似模型中各频率分量的频率值，其定义为ω_l＝2π(f₀+(l-3)f₁)T_s，T_s为接收信号采样间隔； A_l，l＝1,2,...,L和l＝1,2,...,L分别为近似模型中各频率分量的幅值和相位；a_l， l＝1,2,...,L和b_l，l＝1,2,...,L共同组成近似模型中各频率分量的系数，它们与A_l及的关系为： ZPW-2000信号近似模型的矩阵形式定义如下：

式中，为近似模型序列；Z为参数矩阵，定义如下：

Z＝[c(ω₁) c(ω₂) … c(ω_L) s(ω₁) s(ω₂) … s(ω_L)] (5)

式中，c(ω_l)，l＝1,2,...,L和s(ω_l)，l＝1,2,…,L分别为近似模型中各频率分量的余弦、正弦序列，对于频率分量ω_l，c(ω_l)和s(ω_l)的定义分别如下：

c(ω_l)＝[1cos(ω_l) cos(2ω_l) … cos(ω_l(N-1))]^T (6)

s(ω_l)＝[0sin(ω_l) sin(2ω_l) … sin(ω_l(N-1))]^T (7)

α为近似模型中各频率分量的系数向量，定义如下：

α＝[a₁ a₂ … a_L -b₁ -b₂ … -b_L]^T (8)

为求解价值函数的最大值，基于非线性最小二乘法的ZPW-2000信号解调方法的步骤为：

A.计算解调过程中与接收信号不相关的参数，以全局变量的形式存入内存；

B.依次提取内存中针对粗估计的144种载频、低频配置方案对应的参数，采用高斯消元法求解信号在不同载频、低频配置方案下各频率分量幅值，并计算相应的价值函数值。依据各频率分量幅值特征，排除不符合ZPW-2000信号频谱特征的载频、低频配置方案，在剩余方案中，选取价值函数值最大者对应的载频、低频值作为粗估计结果；

C.提取内存中针对网格搜索的相关参数中与载频、低频粗估计结果相对应的参数；采用FFT快速构造网格搜索所需的中间变量，采用高斯消元法求解信号在不同载频、低频配置方案下各频率分量幅值，并计算相应的价值函数值。选取价值函数值最大者对应的载频、低频值作为网格搜索结果；

D.在网格搜索结果附近，采用二维二分搜索算法搜寻价值函数峰值位置，获得载频、低频精确估计值及对应的各频率分量幅值；

E.汇总估计结果，通过条件判断，判断接收信号是否为ZPW-2000信号以及信号是否受到干扰，输出载频、低频及提示信息。

本发明具有以下优点：本发明的解调准确度优于目前常用的周期图法；在高信噪比条件下，本发明可正确解调采样持续时间大于等于0.11s的ZPW-2000移频信号的低频、载频频率，降低了解调方法对采样持续时间的限制；本发明具备实时解调能力，当信号长度为8000时，本发明在MATLAB环境下的平均运行时间为0.32s；本发明的抗白噪声性能优于目前常用的周期图法；本发明给出了具有容错特性的ZPW-2000信号频谱特征检测方法，当信号频谱形态被干扰信号破坏时，可发现译码错误，当信号受信干比大于等于1:3.6的单频干扰影响时，可正确解调。

本发明基于式(2)所示目标函数提出，详细给出了一种实现方法，使得解调方法能够满足实时解调需求。凡是以式(2)为目标函数的ZPW-2000移频信号解调方法都属于本发明的范围。

附图说明

图1为本发明流程图；

图2为本发明实施例1的仿真结果，即在采样持续时间为0.08s，以及以0.01s 为间隔从0.10s枚举至0.19s，以及以0.1s为间隔从0.2s枚举至1s共20种情形下，采用本发明方法，欠采样补零法及Zoom-FFT插值法解调ZPW-2000仿真信号所得的载频估计平均绝对误差结果；其中，每种采样持续时间下共有4320个高斯白噪声信噪比为30dB，采样频率为8000Hz的ZPW-2000仿真信号，涵盖了144种载频、低频配置方案，且对每种方案进行了30次MonteCarlo模拟；因此各解调方法的平均绝对误差为4320个绝对误差的平均值；

图3为本发明实施例1的仿真结果，与图2相对应，是实施例1仿真过程得出的低频估计平均绝对误差；

图4为本发明实施例1的仿真结果，与图2相对应，是实施例1仿真过程得出的错误译码次数统计数据；

图5为本发明实施例2的仿真结果，即当采样持续时间以0.1s为间隔从0.1s 枚举至1s时，采用本发明方法，欠采样补零法及Zoom-FFT插值法解调1000次载频、低频、采样频率分别为2001.4Hz，29Hz，8000Hz的ZPW-2000移频信号所花的平均运行时间。因此任意采样持续时间下各解调方法的平均运行时间为 1000个运行时间的平均值；

图6为本发明实施例3的仿真结果，即当高斯白噪声信噪比以2dB为间隔从-20dB枚举至50dB时，采用本发明方法，欠采样补零法及Zoom-FFT插值法解调ZPW-2000仿真信号所得的载频估计平均绝对误差结果。其中，每种高斯白噪声信噪比下共有4320个采样持续时间为0.3s，采样频率为8000Hz的ZPW-2000 仿真信号，涵盖了144种载频、低频配置方案，且对每种方案进行了30次Monte Carlo模拟。因此任意高斯白噪声信噪比下各解调方法的平均绝对误差为4320个绝对误差的平均值；

图7为本发明实施例3的仿真结果，与图6相对应，是实施例3仿真过程得出的低频估计平均绝对误差；

图8为本发明实施例3的仿真结果，与图6相对应，是实施例3仿真过程得出的错误译码次数统计数据；

图9为本发明实施例4的仿真结果，即当单频干扰的频率f_sfi以1Hz为间隔从(2001.4-40)Hz枚举至(2001.4+40)Hz，幅值A_sfi以0.2为间隔从0.2枚举至7时，采用本发明方法，欠采样补零法及Zoom-FFT插值法解调ZPW-2000信号与单频干扰叠加而成的仿真信号所得的译码错误次数统计结果。其中，ZPW-2000信号的采样频率、载频、低频、采样持续时间、高斯白噪声信噪比分别为8000Hz、 2001.4Hz、29Hz、1s、5dB，且在每种单频干扰情况下进行了30次Monte Carlo 模拟。因此任意频率、幅值的单频干扰下各解调方法的译码错误次数为30次仿真结果的统计值。

具体实施方式

下面对本发明的技术方案进行详细描述。

参考图1，一种基于非线性最小二乘法的ZPW-2000信号解调方法，该方法的实现目的为寻求以下价值函数的最大值：

式中，分别为载频频率f₀、低频频f₁估计值；x＝[x(0) x(1) x(2) … x(N-1)]^T为接收信号采样序列；N为采样序列长度，序列下标从0开始；[·]^T为矩阵转置符号；Z为下文定义的ZPW-2000信号近似模型矩阵形式中的参数矩阵。

由于ZPW-2000信号能量主要集中于载频及一、二次边频分量，依据ZPW-2000信号的傅立叶级数展开形式，定义如下ZPW-2000信号s(n)的近似模型如下：

式中，L为近似模型中频率分量个数，本发明令L＝5；ω_l，l＝1,2,...,L为近似模型中各频率分量的频率值，其定义为ω_l＝2π(f₀+(l-3)f₁)T_s，T_s为接收信号采样间隔； A_l，l＝1,2,...,L和l＝1,2,...,L分别为近似模型中各频率分量的幅值和相位；a_l， l＝1,2,...,L和b_l，l＝1,2,...,L共同组成近似模型中各频率分量的系数，它们与A_l及的关系为： ZPW-2000信号近似模型的矩阵形式定义如下：

式中，为近似模型序列；Z为参数矩阵，定义如下：

Z＝[c(ω₁) c(ω₂) … c(ω_L) s(ω₁) s(ω₂) … s(ω_L)] (5)

式中，c(ω_l)，l＝1,2,...,L和s(ω_l)，l＝1,2,…,L分别为近似模型中各频率分量的余弦、正弦序列，对于频率分量ω_l，c(ω_l)和s(ω_l)的定义分别如下：

c(ω_l)＝[1cos(ω_l) cos(2ω_l) … cos(ω_l(N-1))]^T (6)

s(ω_l)＝[0sin(ω_l) sin(2ω_l) … sin(ω_l(N-1))]^T (7)

α为近似模型中各频率分量的系数向量，定义如下：

α＝[a₁ a₂ … a_L -b₁ -b₂ … -b_L]^T (8)

为求解价值函数的最大值，基于非线性最小二乘法的ZPW-2000信号解调方法的步骤为：

A.计算解调过程中与接收信号不相关的参数，以全局变量的形式存入内存。

B.依次提取内存中针对粗估计的144种载频、低频配置方案对应的参数，采用高斯消元法求解信号在不同载频、低频配置方案下各频率分量幅值，并计算相应的价值函数值。依据各频率分量幅值特征，排除不符合ZPW-2000信号频谱特征的载频、低频配置方案，在剩余方案中，选取价值函数值最大者对应的载频、低频值作为粗估计结果。

C.提取内存中针对网格搜索的相关参数中与载频、低频粗估计结果相对应的参数。采用FFT快速构造网格搜索所需的中间变量，采用高斯消元法求解信号在不同载频、低频配置方案下各频率分量幅值，并计算相应的价值函数值。选取价值函数值最大者对应的载频、低频值作为网格搜索结果。

D.在网格搜索结果附近，采用二维二分搜索算法搜寻价值函数峰值位置，获得载频、低频精确估计值及对应的各频率分量幅值。

E.汇总估计结果，通过条件判断，判断接收信号是否为ZPW-2000信号以及信号是否受到干扰，输出载频、低频及提示信息。

以下是基于非线性最小二乘法的ZPW-2000信号解调方法的详细实施步骤：

A.计算与接收信号不相关参数存入内存

(A.1)首先确认接收信号序列的长度N和采样频率f_s。定义载频序列以及低频序列依次枚举载频和低频，形成144组编号为(p,q)的载频、低频配置方案分别构造对应的矩阵，·^(p,q)表示编号为(p,q)的载频、低频配置方案下的数据。矩阵的参数配置方法参考式(3)和式(5)，其中的第一行元素采用CORDIC算法生成，的剩余元素由三角函数和角公式求得。由可得将144组与以全局变量的形式存入内存，供后续步骤(B.1)使用。

(A.2)定义载频序列2601.4,2598.7)(Hz)以及低频序列依次枚举载频和低频，形成144组编号为(p,q)的载频、低频配置方案对于编号(p,q) 的配置方案，载频和低频的网格搜索间隔均为Δω_AII，定义如下：

式中：

式中，为向下取整符号，意为四舍五入取整；以单位为Hz的实际频率值来描述，网格搜索中载频搜索频率值定义在(Hz)区间内，低频搜索频率值定义在(Hz)区间内；转化为数字频率，载频搜索频率值的定义如下：

式中，为：

当载频的网格搜索频率值确定后，即可以j为系数设定低频搜索值。j以及低频搜索值确定后各频率分量频率值l＝1,2,…,L，L＝5的定义如下所示：

式中，·^(p,q,i,j)表示编号(p,q)的载频、低频配置方案下编号(i,j)的网格搜索值对应的数据；和的定义如下：

将式(16)至式(20)代入式(5)，可得矩阵计算出因此，144组载频、低频配置方案中，每组方案对应一组 M^(p,q)，同时还包括组按次序以全局变量的形式将这些参数存入内存，供后续步骤(C.1)使用。

B.对ZPW-2000信号的载频、低频进行粗估计

(B.1)令J_BImax＝0，f_0c＝-3000(Hz)，f_1c＝-3000(Hz)。确认长度为N的接收信号序列x＝[x(0)x(1)…x(N-1)]^T。依次提取内存中(A.1)步骤存入的144组载频、低频配置方案对应的数据。其中，对于载频、低频配置方案其包含和采用高斯消元法求解方程：

将解得的α_BI代入下式：

式中，A_l，l＝2,3,4为各频率分量幅值，参考式(3)和式(8)由α_BI求得；A_threshold为ZPW-2000信号包络幅值最小值，可实测获得；A_e为高信噪比下各频率分量幅值估计的最大误差，可实测获得，依据仿真结果，令A_e＝0.1；m为调制指数，其定义为r₀₁为载频分量幅值与一次边频分量幅值比值的绝对值，定义如下：

若式(24)结果为真，则将α_BI代入式(26)，得到该载频、低频配置方案的价值函数值：

否则令J_BI＝0；此时，若J_BI>J_BImax，则更新粗估计价值函数及载频、低频粗估计值，即当完成144组载频、低频配置方案的枚举后，若f_0c＝-3000(Hz)或f_1c＝-3000(Hz)，则令并转至步骤E，否则转至步骤C。

C.采用网格搜索获得更精确的低频、载频估计值

(C.1)提取内存中(A.2)存入的数据中满足条件的载频、低频配置方案对应的计算：

ω_shift＝ω_0c-k₀Δω_CI (27)

式中k₀及其定义式中ω_0c和Δω_CI定义分别如下：

定义向量x_shift：

将序列x·x_shift补零至M点，执行M点FFT得X_M＝[X_M(0)X_M(1)X_M(2)… X_M(M-1)]。其中，·为向量点积符号。令J_CImax＝0。依次提取(A.2)存入数据中载频、低频配置方案组以(i,j)为编号的网格搜索数据。对于其中某一网格搜索编号(i,j)，计算η＝[η₁ η₂ η₃ ... η_2L]^T，其中：

式中，Re(·),Im(·)分别为取实部符号和取虚部符号；ξ_l，l＝1,2,…,L，L＝5的定义如下：

ξ₃＝X_M(k₀+i) (35)

采用高斯消元法计算方程：

将α_CI代入下式：

J_CI＝η^Tα_CI (39)

此时，若J_CI>J_CImax，则J_CImax＝J_CI，ω_0g＝ω_0c+iΔω_CI，ω_1g＝(k₁ ^(p,q)+j)Δω_CI，Δω_g＝Δω_CI。当完成(2L_h+1)²组网格搜索值的枚举，转至步骤D。

D.采用精确搜索获得ZPW-2000信号载频、低频精确估计值

(D.1)精确搜索由二维二分搜索算法实现。定义载频频率二分搜索初始区间的左边界ω_0l、右边界ω_0r及区间半长Δω₀定义为：ω_0l＝ω_0g-Δω_g，ω_0r＝ω_0g+Δω_g，Δω₀＝Δω_g。搜索终止条件为2Δω₀≤ω_rp，ω_rp为解调方法的精度要求。载频频率二分搜索步骤如下：

Step 1：以载频频率ω_0l、ω_0r为参数进行低频频率二分搜索，分别得到对应的价值函数值J_0l、J_0r。令ω_0m＝(ω_0l+ω_0r)/2。

Step 2：以载频频率ω_0m为参数进行低频频率二分搜索，得到对应的低频估计值ω_{0m_1}以及对应的频率分量系数α_0m和价值函数值J_0m。令Δω₀＝Δω₀/2，ω_{0m_last}＝ω_0m。若J_0l<J_0r，则J_0l＝J_0m，ω_0m＝ω_0m+Δω₀，否则J_0r＝J_0m，ω_0m＝ω_0m-Δω₀。

Step 3：若2Δω₀≤ω_rp，则结束二分搜索，进而载频、低频及各频率分量系数的精确估计值分别为否则回到 step 2。

低频频率二分搜索步骤与载频频率二分搜索基本一致。相对应地，低频频率二分搜索初始区间的左边界ω_1l、右边界ω_1r及区间半长Δω₁定义为：ω_1l＝ω_1g-Δω_g，ω_1r＝ω_1g+Δω_g，Δω₁＝Δω_g。搜索终止条件为2Δω₁≤ω_rp。令代入的载频频率参数值为ω_0p。低频频率二分搜索步骤如下：

Step 1：以ω_0p为载频频率参数，ω_1l、ω_1r为低频频率参数，代入式(3)和式(5)，得到矩阵Z_1l、Z_1r，参考式(23)和式(26)，结合高斯消元法，得到对应的价值函数值J_1l、J_1r。令ω_1m＝(ω_1l+ω_1r)/2。

Step 2：同样地，参考式(3)、式(5)、式(23)和式(26)求得载频频率ω_0p，低频频率ω_1m对应的价值函数值J_1m以及频率分量系数α_1m。令Δω₁＝Δω₁/2，ω_{1m_last}＝ω_1m。若J_1l<J_1r，则J_1l＝J_1m，ω_1m＝ω_1m+Δω₁，否则J_1r＝J_1m，ω_1m＝ω_1m-Δω₁。

Step 3：若2Δω₁≤ω_rp，则结束低频频率二分搜索，低频频率二分搜索输出的低频估计值，价值函数值以及频率分量系数分别为ω_{1m_last}，J_1m，α_1m，否则回到 step 2。

E.汇总估计结果，通过条件判断，决定输出信息

(E.1)若或则提示信号接收错误，结束解调过程，否则将代入下式：

式中，l＝2,3,4为各频率分量幅值，参考式(3)和式(8)由求得；为载频分量幅值与一次边频分量幅值的比值，参考式(24)和式(25)求得。若式(40)结果为真，则正常显示载频、低频信息结束解调过程，否则显示载频、低频信息的同时，提示信号受到干扰，结束解调过程。

以下从具体仿真实验进一步说明本发明及其效果。

实施例1，为测试本发明的解调准确度，本实施例中，仿真生成了不同采样持续时间，144种载频、低频配置方案下，高斯白噪声信噪比、采样频率分别为 30dB、8000Hz的ZPW-2000仿真信号。其中，高斯白噪声信噪比的定义为σ²为高斯白噪声的方差。采用本发明方法，欠采样补零法及 Zoom-FFT插值法解调仿真信号，得到载频、低频平均绝对误差分别如图2、图 3所示。需指出，实施例1以及以下各实施例中，本发明的ω_rp＝2π·(0.0005/8000)(rad/sample)，同时，通过补零，欠采样补零法及Zoom-FFT插值法中FFT的频域采样间隔为0.0005Hz。可看出，本发明方法的载频、低频估计平均绝对误差均小于欠采样补零法及Zoom-FFT插值法。当采样持续时间为0.3s 时，本发明的载频、低频平均绝对误差分别可低至0.0039Hz及0.0035Hz。因此，本发明的解调准确度优于目前常用的周期图法。

上述仿真过程中，若解调方法输出的载频或低频频率估计值绝对误差超过0.4Hz，则判定为译码错误。通过对每次仿真结果进行统计，在不同采样持续时间下，本发明方法，欠采样补零法及Zoom-FFT插值法译码错误次数统计结果如图4所示。当采样持续时间大于等于0.11s时，本发明译码错误次数可降为0次。相比之下，欠采样补零法及Zoom-FFT插值法都需要0.19s的采样持续时间才能保证信号正确解调。因此，本发明方法显著降低了解调方法对采样持续时间的限制。

实施例2，为测试本发明的实时性能，本实施例中，仿真生成了不同采样持续时间下，载频、低频、采样频率分别为2001.4Hz，29Hz，8000Hz的ZPW-2000 移频信号。采用本发明方法，欠采样补零法及Zoom-FFT插值法分别对每种仿真信号解调1000次，得到平均运行时间如图5所示。需指出，由于步骤A可在接收信号前完成，仿真中记录的运行时间不包括步骤A所花时间。可看出，虽然本发明的平均运行时间高于欠采样补零法及Zoom-FFT插值法，但未超过0.32s。因此本发明能使解调过程的应变时间不大于2s，基本可满足实时解调需求。由于本发明的粗估计、网格搜索过程中计算过程相互独立，利用相关硬件的并行计算功能，可进一步减少运行时间。此外，为了在满足实时解调需求的同时在时间上避开暂态冲击干扰的影响，可首先利用粗估计过程搜寻出粗估计结果能够保持一致的不重叠区间，之后再采用网格搜索及精确搜索获得精确的载频、低频估计值，进而可提高解调方法的可靠性。

实施例3，为测试本发明的抗白噪声能力，本实施例中，仿真生成了不同高斯白噪声信噪比，144种载频、低频配置方案下，采样持续时间、采样频率分别为0.3s、8000Hz的ZPW-2000移频信号。采用本发明，欠采样补零法及Zoom-FFT 插值法解调仿真信号，得到载频、低频平均绝对误差分别如图6、图7所示。可看出，相比于Zoom-FFT插值法，本发明需在信噪比大于等于-2dB后，载频、低频估计平均绝对误差才趋于稳定。这是因为本发明中，式(25)的判决条件可以甄别出部分译码错误情况，输出-3000Hz的载频、低频估计结果，因此部分译码错误情况不会被求平均过程掩盖，因此在信噪比小于-2dB时，本发明的载频、低频平均绝对误差高于Zoom-FFT插值法。

上述仿真过程中，若解调方法输出的载频或低频频率估计值绝对误差超过0.4Hz，则判定为译码错误。通过对每次仿真结果进行统计，在不同高斯白噪声信噪比下，本发明，欠采样补零法及Zoom-FFT插值法译码错误次数统计结果如图8所示。可看出，本发明在低信噪比下错误译码次数要小于前两种方法，在信噪比大于等于2dB时，本发明错误译码次数降为0。相比之下，欠采样补零法及 Zoom-FFT插值法分别需在信噪比大于等于16dB、4dB时方可保证正确解调。因此，本发明的抗白噪声能力优于欠采样补零法及Zoom-FFT插值法。同时，由图 6、图7与图8相比较可知，本发明中式(25)的判决条件尚不能甄别出所有的译码错误情况，因此当白噪声干扰较严重时，本发明的解调结果不再可信。

实施例4，由于29Hz的ZPW-2000信号边频幅值最小，最易受单频干扰影响，为测试本发明的抗单频干扰能力，本实施例仿真生成了低频频率为29Hz的 ZPW-2000信号以及幅值、频率分别为A_sfi、f_sfi的单频正弦干扰。两者相叠加，采用本发明，欠采样补零法及Zoom-FFT插值法解调该信号。当解调方法输出的载频或低频频率估计值绝对误差超过0.4Hz，则判定为译码错误。由此得到在不同幅值、频率单频干扰下各方法的译码错误次数统计情况如图9所示。图中，每个像素点颜色越深代表译码错误次数越多。可看出，当A_sfi小于等于3.6时，本发明的译码错误次数均为0，可保证信号正确解调。因此，当信号受信干比大于等于1:3.6的单频干扰影响时，本发明可正确解调信号。

本发明是以文献"Fast fundamental frequency estimation:Making astatistically efficient estimator computationally efficient"(Nielsen J K,Jensen T L,Jensen J R,et al.Signal Processing,135:188-197,2017)所提的非线性最小二乘法为基础，结合 ZPW-2000移频信号特点提出一系列改进措施，而最终完成。此外，本发明参考了文献“A survey of CORDIC algorithms for FPGA based computers”(Andraka R. Proceedings of the 1998ACM/SIGDA sixth international symposium onField programmable gate arrays:191-200,1998)提出的CORDIC算法，用于计算任意角度的sin函数和cos函数值，以及文献“Linear algebra and its application 4th edition”(G. Strang.Cengage Learning:1-64,2006)介绍的高斯消元法，用于求解线性方程组。

Claims (7)

1.一种基于非线性最小二乘法的ZPW-2000信号解调方法，其特征在于，该方法的实现目的为寻求以下价值函数的最大值：

式中，分别为载频频率f0、低频频f1估计值；x＝[x(0) x(1) x(2) … x(N-1)]^T为接收信号采样序列；N为采样序列长度，序列下标从0开始；[·]^T为矩阵转置符号；Z为下文定义的ZPW-2000信号近似模型矩阵形式中的参数矩阵；

依据ZPW-2000信号的傅立叶级数展开形式，定义ZPW-2000信号s(n)的近似模型如下：

式中，L为近似模型中频率分量个数，本发明令L＝5；ω_l，l＝1,2,...,L为近似模型中各频率分量的频率值，其定义为ω_l＝2π(f₀+(l-3)f₁)T_s，T_s为接收信号采样间隔；A_l，l＝1,2,...,L和分别为近似模型中各频率分量的幅值和相位；a_l，l＝1,2,...,L和b_l，l＝1,2,...,L共同组成近似模型中各频率分量的系数，它们与A_l及的关系为：ZPW-2000信号近似模型的矩阵形式定义如下：

式中，为近似模型序列；Z为参数矩阵，定义如下：

Z＝[c(ω₁) c(ω₂) … c(ω_L) s(ω₁) s(ω₂) … s(ω_L)] (5)

式中，c(ω_l)，l＝1,2,...,L和s(ω_l)，l＝1,2,…,L分别为近似模型中各频率分量的余弦、正弦序列，对于频率分量ω_l，c(ω_l)和s(ω_l)的定义分别如下：

c(ω_l)＝[1 cos(ω_l) cos(2ω_l) … cos(ω_l(N-1))]^T (6)

s(ω_l)＝[0 sin(ω_l) sin(2ω_l) … sin(ω_l(N-1))]^T (7)

α为近似模型中各频率分量的系数向量，定义如下：

α＝[a₁ a₂ … a_L -b₁ -b₂ … -b_L]^T (8) 。

2.根据权利要求1所述的一种基于非线性最小二乘法的ZPW-2000信号解调方法，其特征在于，为求解价值函数的最大值，基于非线性最小二乘法的ZPW-2000信号解调方法的步骤为：

A.计算解调过程中与接收信号不相关的参数，以全局变量的形式存入内存；

B.依次提取内存中针对粗估计的144种载频、低频配置方案对应的参数，采用高斯消元法求解信号在不同载频、低频配置方案下各频率分量幅值，并计算相应的价值函数值。依据各频率分量幅值特征，排除不符合ZPW-2000信号频谱特征的载频、低频配置方案，在剩余方案中，选取价值函数值最大者对应的载频、低频值作为粗估计结果；

C.提取内存中针对网格搜索的相关参数中与载频、低频粗估计结果相对应的参数；采用FFT快速构造网格搜索所需的中间变量，采用高斯消元法求解信号在不同载频、低频配置方案下各频率分量幅值，并计算相应的价值函数值。选取价值函数值最大者对应的载频、低频值作为网格搜索结果；

D.在网格搜索结果附近，采用二维二分搜索算法搜寻价值函数峰值位置，获得载频、低频精确估计值及对应的各频率分量幅值；

E.汇总估计结果，通过条件判断，判断接收信号是否为ZPW-2000信号以及信号是否受到干扰，输出载频、低频及提示信息。

3.根据权利要求2所述的一种基于非线性最小二乘法的ZPW-2000信号解调方法，其特征在于，所述步骤A的过程为：

A.1首先确认接收信号序列的长度N和采样频率f_s；定义载频序列以及低频序列依次枚举载频和低频，形成144组编号为(p,q)的载频、低频配置方案分别构造对应的矩阵，·(p,q)表示编号为(p,q)的载频、低频配置方案下的数据；矩阵的参数配置方法参考式(3)和式(5)，其中的第一行元素采用CORDIC算法生成，的剩余元素由三角函数和角公式求得；由可得将144组与以全局变量的形式存入内存，供后续步骤(B.1)使用；

A.2定义载频序列 以及低频序列依次枚举载频和低频，形成144组编号为(p,q)的载频、低频配置方案对于编号(p,q)的配置方案，载频和低频的网格搜索间隔均为Δω_AII，定义如下：

式中，M^(p,q)及其定义式中的N_Δ和Δω_d定义分别如下：

式中，为向下取整符号，意为四舍五入取整；以单位为Hz的实际频率值来描述，网格搜索中载频搜索频率值定义在区间内，低频搜索频率值定义在区间内；转化为数字频率，载频搜索频率值的定义如下：

其中，定义为：

当载频的网格搜索频率值确定后，即可以j为系数设定低频搜索值；j以及低频搜索值确定后各频率分量频率值的定义如下所示：

式中，·(p,q,i,j)表示编号(p,q)的载频、低频配置方案下编号(i,j)的网格搜索值对应的数据；和的定义如下：

将式(16)至式(20)代入式(5)，可得矩阵计算出因此，144组载频、低频配置方案中，每组方案对应一组M^(p,q)，同时还包括组按次序以全局变量的形式将这些参数存入内存，供后续步骤(C.1)使用。

4.根据权利要求2所述的一种基于非线性最小二乘法的ZPW-2000信号解调方法，其特征在于，所述步骤B的过程为：

(B.1)令J_BImax＝0，f_0c＝-3000Hz，f_1c＝-3000Hz；确认长度为N的接收信号序列x＝[x(0)x(1)…x(N-1)]^T；依次提取内存中(A.1)存入的144组载频、低频配置方案对应的数据；其中，对于载频、低频配置方案其包含和采用高斯消元法求解方程：

将解得的α_BI代入下式：

式中，A_l，l＝2,3,4为各频率分量幅值，参考式(3)和式(8)由α_BI求得；A_threshold为ZPW-2000信号包络幅值最小值，可实测获得；A_e为高信噪比下各频率分量幅值估计的最大误差，可实测获得，依据仿真结果，令A_e＝0.1；m为调制指数，其定义为r₀₁为载频分量幅值与一次边频分量幅值比值的绝对值，定义如下：

若式(24)结果为真，则将α_BI代入式(26)，得到该载频、低频配置方案的价值函数值：

否则令J_BI＝0；此时，若J_BI>J_BImax，则更新粗估计价值函数及载频、低频粗估计值，即J_BImax＝J_BI，当完成144组载频、低频配置方案的枚举后，若f_0c＝-3000(Hz)或f_1c＝-3000Hz，则令并转至步骤E，否则转至步骤C。

5.根据权利要求2所述的一种基于非线性最小二乘法的ZPW-2000信号解调方法，其特征在于，所述步骤C的过程为：

C.1提取内存中A.2步骤存入的数据中满足条件的载频、低频配置方案对应的M^(p,q)，计算：

ω_shift＝ω_0c-k₀Δω_CI (27)

式中k₀及其定义式中ω_0c和Δω_CI定义分别如下：

定义向量x_shift：

将序列x·x_shift补零至M点，执行M点FFT得X_M＝[X_M(0) X_M(1) X_M(2) … X_M(M-1)]；其中，·为向量点积符号；令J_CImax＝0；依次提取A.2步骤存入数据中载频、低频配置方案组以(i,j)为编号的网格搜索数据；对于其中某一网格搜索编号(i,j)，计算η＝[η₁ η₂ η₃ ... η_2L]^T，其中：

式中，Re(·),Im(·)分别为取实部符号和取虚部符号；ξ_l，l＝1,2,…,L，L＝5的定义如下：

ξ₃＝X_M(k₀+i) (35)

采用高斯消元法计算方程：

将α_CI代入下式：

J_CI＝η^Tα_CI (39)

此时，若J_CI>J_CImax，则J_CImax＝J_CI，ω_0g＝ω_0c+iΔω_CI，Δω_g＝Δω_CI；当完成(2L_h+1)²组网格搜索值的枚举，转至步骤D。

6.根据权利要求2所述的一种基于非线性最小二乘法的ZPW-2000信号解调方法，其特征在于，所述步骤D的过程为：

D.1精确搜索由二维二分搜索算法实现。定义载频频率二分搜索初始区间的左边界ω_0l、右边界ω_0r及区间半长Δω₀定义为：ω_0l＝ω_0g-Δω_g，ω_0r＝ω_0g+Δω_g，Δω₀＝Δω_g；搜索终止条件为2Δω₀≤ω_rp，ω_rp为解调方法的精度要求；载频频率二分搜索步骤如下：

Step 1：以载频频率ω_0l、ω_0r为参数进行低频频率二分搜索，分别得到对应的价值函数值J_0l、J_0r。令ω_0m＝(ω_0l+ω_0r)/2；

Step 2：以载频频率ω_0m为参数进行低频频率二分搜索，得到对应的低频估计值ω_{0m_1}以及对应的频率分量系数α_0m和价值函数值J_0m；令Δω₀＝Δω₀/2，ω_{0m_last}＝ω_0m；若J_0l<J_0r，则J_0l＝J_0m，ω_0m＝ω_0m+Δω₀，否则J_0r＝J_0m，ω_0m＝ω_0m-Δω₀；

Step 3：若2Δω₀≤ω_rp，则结束二分搜索，进而载频、低频及各频率分量系数的精确估计值分别为否则回到step 2；

低频频率二分搜索步骤与载频频率二分搜索基本一致；相对应地，低频频率二分搜索初始区间的左边界ω_1l、右边界ω_1r及区间半长Δω₁定义为：ω_1l＝ω_1g-Δω_g，ω_1r＝ω_1g+Δω_g，Δω₁＝Δω_g；搜索终止条件为2Δω₁≤ω_rp；令代入的载频频率参数值为ω_0p。低频频率二分搜索步骤如下：

Step 1：以ω_0p为载频频率参数，ω_1l、ω_1r为低频频率参数，代入式(3)和式(5)，得到矩阵Z_1l、Z_1r，参考式(23)和式(26)，结合高斯消元法，得到对应的价值函数值J_1l、J_1r。令ω_1m＝(ω_1l+ω_1r)/2；

Step 2：同样地，参考式(3)、式(5)、式(23)和式(26)求得载频频率ω_0p，低频频率ω_1m对应的价值函数值J_1m以及频率分量系数α_1m。令Δω₁＝Δω₁/2，ω_{1m_last}＝ω_1m；若J_1l<J_1r，则J_1l＝J_1m，ω_1m＝ω_1m+Δω₁，否则J_1r＝J_1m，ω_1m＝ω_1m-Δω₁；

Step 3：若2Δω₁≤ω_rp，则结束低频频率二分搜索，低频频率二分搜索输出的低频估计值，价值函数值以及频率分量系数分别为ω_{1m_last}，J_1m，α_1m，否则回到step 2。

7.根据权利要求2所述的一种基于非线性最小二乘法的ZPW-2000信号解调方法，其特征在于，所述步骤E的过程为：

E.1若或则提示信号接收错误，结束解调过程，否则将代入下式：

式中，为各频率分量幅值，参考式(3)和式(8)由求得；为载频分量幅值与一次边频分量幅值的比值，参考式(24)和式(25)求得；若式(40)结果为真，则正常显示载频、低频信息结束解调过程，否则显示载频、低频信息的同时，提示信号受到干扰，结束解调过程。