0 引 言
在天波超视距雷达信号处理领域中，多径时延估计技术是核心技术之一。短波信号在传播过程中到达电离层时，由于电离层不同高度的电子密度分布存在差异，信号会在不同层高发生反射，进而形成多条传播路径，最终经这些路径到达接收端时因传播距离不同而产生时间延迟，引发多径效应。雷达受多径效应的影响在实际探测中产生假目标，影响对真实目标的识别与定位，准确估计多径时延信息是有效消除假目标、提升雷达探测精度，为后续信号分离与干扰抑制奠定基础的关键所在。
目前，主流时延估计方法主要有互相关法、子空间法及稀疏优化方法。互相关法原理简单，易于工程实现，Knapp等提出了广义互相关（Generalized Cross Correlation,GCC）算法，GCC算法通过对接收到的信号进行加权处理，使得参考信号与接收信号的相关峰更加尖锐，提高时延的估计精度，该方法在工程领域得到了广泛的应用［1］。子空间法及稀疏优化方法估计精度高，但计算复杂度高，Oziewicz将多重信号分类（Multiple Signal Classification,MUSIC）算法应用于正交频分复用（Orthogonal Frequency Division Multiplexing, OFDM）系统的时延估计问题［2］。研究者主要集中于提升算法的估计精度，并随之设计了一系列复杂算法。但对于计算资源有限的天波超视距雷达系统而言，难以满足系统对实时性的需求［3-5］。在对快速多径时延估计算法的研究中，研究者发现多径信道由多个复指数分量构成，精确估计其复指数分量的频率信息［6］，能够实现对多径时延的有效提取。Liu等提出了一种基于快速傅里叶变换（Fast Fourier Transform,FFT）及离散傅里叶变换（Discrete Fourier Transform,DFT）插值的正弦信号频率估计算法，相较于传统的优化算法，在计算效率上实现了显著提升［7］。陈李等基于信号的谱稀疏特性，提出了一种基于稀疏傅里叶变换的快速频谱分析方法，该方法利用延时采样的相位旋转效应，在低采样率下实现了对宽带信号的快速频谱感知［8］。杨君将一维的稀疏傅里叶变换扩展到二维非均匀稀疏分数傅里叶变换，利用信号的稀疏性构建了非均匀采样下的二维离散分数傅里叶变换快速算法［9］。离散余弦变换（Discrete Cosine Transform, DCT）相较于傅里叶变换在处理一些具有特定特征的信号时，能够将信号能量高度集中在少数变换稀疏上，其快速算法实现更加简单且计算量更小［10-11］。Serbes设计了一种快速的正弦频率估计算法，通过使用移位技术和DFT插值器对DFT稀疏进行插值，在两次迭代中完成收敛［12］并进一步对多个正弦曲线的移位DFT稀疏进行插值，在选取适当参数时，所提算法的估计误差达到了克拉美罗下界，且总计算成本为O（KNlogN）［13］。
尽管现有诸多算法按估计方法在不同方面展现出优势，但在计算效率和估计精度的平衡，尤其针对天波超视距雷达这类计算资源有限且对实时性要求较高的场景，仍存在不足［14］。针对在计算资源有限的情况下，实现快速准确的多径时延估计这一问题，我们提出了基于稀疏离散余弦变换（Sparse Discrete Cosine Transform,SDCT）的多径时延估计算法。基于多径信号的路径稀疏特性，对应其多径信道向量由有限谐波分量组成，实现了对多径时延的快速估计，同时有效提升了算法的时延分辨率，实现了在更小时延间隔下对多径时延的准确估计。
1 信号模型
天波超视距雷达进行探测时，信号由发射端发出后，并非沿单一路径直射至接收端，而是依托电离层这一特殊传播介质，与电离层不同高度分层、不均匀体等产生复杂交互。信号会经电离层的E层、F层等不同区域反射，还可能出现多次反射，形成数量有限的传播路径［15］。天波超视距雷达接收端天线接收到的信号，是经多条路径传输的信号在时域、频域上的线性叠加，其中既包含携带目标信息的多径有用分量，又会混入电离层噪声、杂波干扰等，给后续目标检测、多径时延参数估计及抗干扰处理带来诸多挑战。具体天波超视距雷达信号多径传播模型如图1所示。
图1 天波超视距雷达探测视图
在实际的天波超视距雷达探测场景中，信号回波在电离层不同高度发生反射，沿不同路径到达接收端，由于电离层的层状特性以及电子密度分布的非均匀特性，信号的发射位置主要集中在E层与F层，信号差传播路径呈现出显著的稀疏特性［15-17］，即信号能量主要集中在有限的传播路径上。这种能量集中的现象使得电离层多径信道表现出显著的稀疏特性［18-20］，即能量较弱的路径可以忽略不计，只有少数几条主导路径对接收信号具有实质性贡献，接收端所获取的时域信号可以表示为发射信号经过不同延迟和衰减之后的线性叠加。进一步考虑到电离层的电子密度随日照、季节和地理位置动态变化，信号传播路径参数实际上具有时变特性，因此，可以将接收信号的时域模型表示为
式中，y(t)表示接收信号，s(t)表示发射信号，n(t)表示噪声，αl(t)表示第l条路径的衰减系数,τl(t)表示第l条路径的真实时延，L(t)表示显著多径的数量，αl(t)、τl(t)、L(t)均为时间相关函数，随电离层波动实时变化。通过实测数据表明［18-19］，在典型的天波传播环境下，显著路径的数量通常不超过4～6条，满足L(t) ≪N，其中N表示信号的采样点数，满足稀疏性条件。同时，在一个相干处理间隔内，一般认为在几十毫秒到几秒的跨度内，电离层参数变化缓慢，可以近似认为该系统是一个时不变系统，即信号传播路径参数在一个相干处理间隔内保持不变［20-21］。因此，我们将相干时间内接收信号的时域模型简化为
利用傅里叶变换的线性性质和时移性质，对式（2）进行傅里叶变换，得到接收信号的频域模型：
式中，Y( f)、S( f)、N( f)分别表示接收信号、发射信号与噪声的频域形式，
表示第l条路径的时延带来的相位偏移信息。
信号经过电离层反射后，信号能量集中在少数路径上，基于该特性，可以将多径信道的信道频率响应（Channel Frequency Response, CFR）表示为L个复指数分量的叠加，利用接收信号和发射信号的频域特性，提取多径CFR：
式中，H( f)表示多径CFR，W( f) = Y( f)/X( f)。
进一步，对连续多径CFR进行离散采样，得到离散频率点fk = kΔf,k = 1,2,…,N，Δf表示频率采样间隔，得到多径CFR的离散形式［21］：
进一步，将多径CFR表示为一个稀疏向量H ∈ℂN，其中，H中非零元素的位置对应时延τl，非零元素的值对应衰减系数αl。从稀疏信号处理理论出发，信道向量H的稀疏特性可以定义为［22］
式中，K表示信号显著路径的最大数量，在一个相干处理间隔内，多径稀疏度保持稳定［20-21］，满足‖ H ‖0 = L,L ≤K。
2 离散余弦变换
离散余弦变换是一种众所周知的技术，广泛应用于信号与图像处理。与DFT相比，DCT可以通过消除基于DFT方法中不连续边缘的影响来减少变换域中的高频分量。因此，与傅里叶变换域方法相比，DCT在低频区域具有更好的功率集中性，对原始信道冲激的频率响应有更好的频率逼近，并且混叠误差更低。在应用离散余弦变换进行多径时延估计时，需推导DCT系数与多径时延的关系，当多径信道为无噪信道，提取多径CFR的实部序列：
式中，Re{·} 表示对·取实部，X[k]表示H[k]的实部，αl为复增益，包含幅度和载波相位信息。令αl = al + blj，则第l条路径CFR的实部为
因此，得到多径信道频率响应H[k]实部信息：
采用DCT-II变换，对于长度为N的离散序列X［k］，k=0,1,…,N-1进行DCT变换，其DCT-II系数定义如下：
式中，C[m]表示X[k]的DCT变换结果，α(m)表示归一化因子，用于保持DCT与IDCT的正交性和单位能量，定义为
进一步，将多径CFR的实部代入式（10）中，此时C[m]的表达式为
为了简化计算与分析，交换求和顺序：
定义求和项S1,l(m)和S2,l(m)：
令
，利用三角函数积化和差公式得
式（15）中的求和项是正余弦函数的几何级数，利用几何级数的求和性质得到：
由正余弦函数几何级数求和的性质可知，当wl + φm = 2πn(n ∈ℤ)时，S1,l(m)和S2,l(m)达到最大，因此其最值对应出现在：
对于m = 2NΔfl + 2Nn(n ∈ℤ)而言，一个时延为τl的多径分量，会在DCT域中周期性地产生峰值。但是对于一个天波超视距雷达系统，其最大不模糊时延τmax = 1/Δf，对于任意时延超过τmax的多径信号，都会在后续信号处理过程中产生模糊，无法分辨其真实时延，因此对于m = 2NΔfl +2Nn(n ∈ℤ)，有且只有n = 0这一个解，对于所有|n |＞0(n ∈ℤ)的解所对应的时延，必然已经超出该系统最大不模糊时延范围，属于物理上不可实现的模糊解，基于该特性，可以得到多径CFR经过DCT变换后的峰值位置：
3 基于稀疏余弦变换的快速多径时延估计算法
如式（5）所示，多径信道的CFR可表示为具有稀疏特性的冲激响应之和。本文基于该特性，提出一种稀疏离散余弦变换算法实现对多径时延信息的快速、准确估计。具体流程如图2所示：①信号先经过随机置换，打乱原始谱结构特征；②加窗滤波，抑制噪声并提取信号在DCT域中的峰值信息；③利用子采样和余弦变换提取稀疏峰值并分配到对应区间中；④通过哈希映射返回原始峰值信息并通过迭代统计分析，选择概率最高的L个峰值并输出L稀疏谱。根据L稀疏谱中的位置，结合式（18）中时延和DCT系数的关系即可定位多径时延。
图2 稀疏余弦变换示意图
3.1 随机置换
根据文献［22］，伸缩变换因子σ 与变换长度K互质，且τ ∈{1,2,…,K-1}，定义置换k →σk+τmodK，多径CFR上的相关置换表示为
该操作将余弦域中的稀疏系数随机性映射到新位置，并将重排后的频谱均匀划分为B个区间，每个区间包含K B个元素，对应的置换方式为
式中，xp(i)表示对Xp[k]进行DCT变换的结果，x(i)表示对X[k]进行DCT变换的结果。
3.2 加窗滤波
对随机置换后的信号进行加窗滤波，抑制噪声并提取信号在DCT域中的峰值信息。利用辛格函数对信号能量的聚焦特性及高斯函数的旁瓣衰减特性，将辛格函数与高斯函数相乘设计窗函数G［k］：
窗函数的响应曲线更平滑，旁瓣相对较低，能够更好地识别和区分不同频率分量，从而减少干扰引起的误差。将窗函数与变换后的信号相乘，利用窗函数聚焦信号的非零分量：
式中，Z[k]表示对Xp[k]进行加窗滤波后的结果。
3.3 子采样与离散余弦变换
将Z[k]划分为B个区间，在每个区间内进行子采样提取关键信号特征：
式中，M = K/B，表示每个区间的宽度，Q[i]表示按间隔M对Z[k]采样的求和结果。对子采样信号进行傅里叶逆变换得到其时域形式：
式中，q[i]表示Q[i]的DCT变换结果，z[Mi]表示对Z[k]的傅里叶逆变换结果按间隔M采样的结果。
3.4 哈希映射
提取B个区间内的峰值索引集Γ，并将区间内对应元素作为潜在峰值集Λ：
式中，findpeak(·)表示搜索局部峰值。对集合Λ 中的元素进行哈希映射，得到真实峰值位置集P：
循环执行式（19）至式（26）的操作，得到综合集合：
式中，PT表示第T次迭代得到的峰值索引集，Π 表示第T次迭代的统计结果，其索引对应时域信道冲激响应中出现次数最高的前L个峰值位置，利用式（18）中峰值与DCT系数的关系，得到多径时延估计结果
：
4 仿真结果与分析
本节将SDCT算法与经典DCT算法、FFT算法及SFFT算法进行对比，通过评估4种算法的多径时延估计结果、不同信噪比下的多径时延估计的平均绝对误差（Mean Absolute Error, MAE）以及不同数据长度下的运行时间，验证所提SDCT算法的性能，定义MAE如下：
式中，τl表示真实时延，
表示估计时延，T表示蒙特卡洛模拟次数，一般设置为1 000。
为了评估不同算法的时延分辨率，设置信号长度N为4 096采样点，信号带宽为20 kHz，采样率为120 kHz。多径数量L分别设置为2，信噪比设为10 dB。设置多径时延参数依次为1.038 8和1.068 8 ms，幅度参数依次为0.6+0.7j、0.8+0.3j。算法的多径时延估计结果对比如图3所示，相较于FFT、SFFT、DCT算法，SDCT算法的归一化幅度曲线更为尖锐且集中，其估计结果更接近真实时延。在小时延间隔场景下，SDCT算法能精准捕捉多径时延特征，有效区分邻近时延分量，展现出更优的多径时延分辨率，可实现小时延间隔下精确的多径时延估计，在多径时延估计任务中表现出显著优势。
图3 多径时延估计结果（Δτ = 0.03 ms）
为了分析在不同多径数下的算法性能，我们将SNR设置为10 dB，多径数L依次设置为2、3、4、5，对应多径时延τ依次设置为0.52、3.52、1.99、3.03和2.79 ms，对应幅度系数α 依次设置为0.6+0.7j、0.8+0.3j、0.6+0.5j、0.6+0.4j、0.1+0.5j，对 比 所 提SDCT算法与DCT、SFFT、FFT算法的估计性能。图4（a）、（b）、（c）、（d）分别展示了多径数L为2、3、4、5四类算法的估计结果。从图中可以看出，所提SDCT算法在所有多径数下均能准确估计出所有多径时延，且估计值与真实时延值高度重合，表现出优异的估计精度，DCT、SFFT、FFT算法也实现了对不同多径设置下的精确参数估计，实验结果表明，四类算法的性能随多径数量的变化具有鲁棒性，为进一步对比算法时延估计精度，将对不同多径数下四类算法的MAE随SNR的变化关系进行仿真，通过对比四类算法的MAE评估算法的估计精度。
图4 不同稀疏度下多径时延估计结果（Δτ = 0.03 ms）
我们进一步进行蒙特卡洛实验，对比在不同多径数下四类算法MAE随SNR的变化情况，评估所提SDCT算法的估计性能，将SNR设置为-20 dB至30 dB间隔5 dB，多径数、多径时延及多径幅度设置不变。图5（a）、（b）、（c）、（d）分别展示了多径数L为2、3、4、5四类算法的MAE-SNR曲线图，通过对比在不同稀疏度下的四类算法的MAE随SNR的变化情况可知，在多径数改变时，SDCT算法保持了稳定的性能，且所提的SDCT算法在所有信噪比区间内的MAE均小于FFT、SFFT、DCT算法，其优势在信噪比较低时更为突出，这表明SDCT算法在具有更高的估计精度的同时，也具有较强的抗噪声能力。
图5 不同多径数下各算法的MAE随SNR的变化曲线
最后对比在不同信号规模下的算法运行时间，不改变其他参数，以2为倍数，设置信号长度为211至216。采用算法运行时间作为效率评估指标进行蒙特卡洛仿真，并在同一平台和环境中进行比较。从图6结果可知，在所有信号长度内，SDCT算法的运行时间始终低于DCT算法，且具有一个量级的提升，这说明SDCT算法在运行效率上具有一定优势。但是相较于FFT及SFFT算法，其运行时间更高，且随着数据规模的增大，运行时间的差异进一步提升。
图6 不同信号长度下各个算法运行时间
5 结束语
在本文的研究中，我们推导了多径CFR的DCT变换与多径时延的关系，基于多径信号的路径稀疏特性，提出了一种稀疏离散余弦变换算法，实现了对多径时延的快速准确估计。该算法针对基于多径信道在离散余弦变换域的稀疏信息进行设计，通过随机置换、加窗滤波、子采样离散余弦变换和哈希映射的方式估计多径时延，有效提升了算法的时延分辨率，同时提升了算法的计算效率，为复杂电离层环境下天波超视距雷达多径干扰抑制提供了有力的技术支撑。
［1］KNAPP C, CARTER G. The Generalized Correlation Method for Estimation of Time Delay［J］. IEEE Trans on Acoustics, Speech, and Signal Processing, 1976, 24（4）:320-327.
［2］OZIEWICZ M. On Application of MUSIC Algorithm to Time Delay Estimation in OFDM Channels［J］.IEEE Trans on Broadcasting,2005,51（2）:249-255.
［3］YE Hao,LI G Y,JUANG B H.Power of Deep Learning for Channel Estimation and Signal Detection in OFDM Systems［J］. IEEE Wireless Communications Letters, 2018, 7（1）:114-117.
［4］李雪梅，陶然，王越.时延估计技术研究［J］.雷达科学与技术，2010，8（4）:362-371.
［5］秘超鹏，杨松岩，张春杰，等.时差相差联合估计的多站信号合成增强方法［J］.舰船电子工程，2025，45（4）：44-50.
［6］张宇，景鑫磊，蒋忠进.一种基于稀疏贝叶斯学习的离网DOA 估计算法［J］.雷达科学与技术，2024，22（1）:35-42.
［7］LIU Jinyu,FAN Lei,JIN Jiyu,et al.An Accurate and Efficient Frequency Estimation Algorithm by Using FFT and DTFT［C］//2020 39th Chinese Control Conference, Shenyang,China:IEEE,2020:2913-2917.
［8］陈李，叶芃，杨慧芝，等.基于稀疏傅里叶变换的快速频谱感知方法［J］.电子科技大学学报，2025，54（2）：274-279.
［9］杨君.二维非均匀稀疏分数傅里叶变换算法研究［D］.西安：西安电子科技大学，2023.
［10］龚岩.基于稀疏傅里叶变换的电子通信弱信号采集系统设计［J］.电子技术与软件工程，2022（19）：113-116.
［11］田琳，胡津健.稀疏短时傅里叶变换谱分解方法及应用［J］.地球物理学进展，2021，36（6）：2581-2587.
［12］SERBES A. Fast and Efficient Sinusoidal Frequency Estimation by Using the DFT Coefficients［J］.IEEE Trans on Communications,2019,67（3）:2333-2342.
［13］DAVEY S J,FABRIZIO G A, RUTTEN M G. Detection and Tracking of Multipath Targets in Over-the-Horizon Radar［J］.IEEE Trans on Aerospace and Electronic Systems,2019,55（5）:2277-2295.
［14］SU Hongtao, LIU Hongwei, SHUI Penglang, et al. Estimation of the Doppler Frequency and Direction of Arrival of the Ionospherically Propagated HF Signals［J］. Radio Science,2009,44（4）:1-10.
［15］XU Chen, CAI Hongtao, GAO Shunzu, et al. A Method for HF Skywave Source Geolocation in Unknown Ionosphere Environments and Experimental Results［J］.IEEE Antennas and Wireless Propagation Letters, 2023,22（5）:1059-1063.
［16］TANG Wenqi, WAN Qun, QU Wenzhong, et al. Multipath Measurements for HF Signals with Polarimetric Array and Band-Limited Gaussian Noise Waveform［J］. IEEE Sensors Journal,2025,25（3）:5139-5151.
［17］YANG Xuguang, WANG Mingjie, HUANG Weimin, et al. Experimental Observation and Analysis of Ionosphere Echoes in the Mid-Latitude Region of China Using High-Frequency Surface Wave Radar and Ionosonde［J］.IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing，2020,13:4599-4606.
［18］YANG Longquan, FAN Junmei, GUO Lixin, et al. Simulation Analysis and Experimental Study on the Echo Characteristics of High-Frequency Hybrid Sky-Surface Wave Propagation Mode［J］.IEEE Trans on Antennas and Propagation,2018,66（9）:4821-4831.
［19］WANG Fen, OU Ming, GUO Yaping, et al. Validation of IRI and IRTAM Models Based on Global Ionosphere Radio Observatory Data［C］//2024 14th International Symposium on Antennas, Propagation and EM Theory, Hefei,China:IEEE,2024:1-4.
［20］VENKATA RATNAM D,SARMA A D.Modeling of Low-Latitude Ionosphere Using GPS Data with SHF Model［J］. IEEE Trans on Geoscience and Remote Sensing,2012,50（3）:972-980.
［21］HASSANIEH H, INDYK P, KATABI D, et al. Simple and Practical Algorithm for Sparse Fourier Transform［C］//In Proceedings of the Twenty-third Annual ACMSIAM Symposium on Discrete Algorithms, Kyoto, Japan:Society for Industrial and Applied Mathematics, 2012:1183-1194.
［22］TAO Ran, LI Xuemei, LI Yanlei, et al. Time-Delay Estimation of Chirp Signals in the Fractional Fourier Domain［J］. IEEE Trans on Signal Processing, 2009, 57（7）:2852-2855.