基于模拟退火算法的频域受限序列搜索算法

Spectrally-Constrained Sequence Search Based on Simulated Annealing Algorithm

摘要: 低自相关和低峰值平均功率比的序列是通信或信号处理应用中所需要的。本文研究了在频谱约束下具有低自相关和低峰均比的序列的合成。频谱约束限制了每个子载波上的最大允许功率，以避免特定保留频带上的干扰。为了获得性能更好的频域受限序列，本文提出采用模拟退火算法来对序列进行数值优化。基本原理是为序列搜索建立适当的目标函数，调试出适当的退火和停止规则。数值实验结果表明，该算法在满足频谱约束的条件下，具有更好的抑制PAPR和非周期自相关的能力，且该算法收敛性良好。

Abstract: Sequences with low Autocorrelation (AC) and low Peak to Average Power Ratio (PAPR) are desired in many communications or signal processing applications. This work investigates the synthesis of sequence with the desired low AC and low PAPR under spectral constraints. The spectral constraints limit the maximum allowable power on each subcarrier to avoid interference on particular reserved bands. In order to obtain a Spectrally-Constrained Sequence (SCS) with better performance, Simulated Annealing (SA) algorithm is proposed to optimize the sequence numerically. The basic principle is to establish the appropriate objective function for sequence search and debug the appropriate annealing and stopping rules. Numerical results indicate that this algorithm features better suppression capabilities for both PAPR and aperiodic auto-correlations under the condition of spectrum constraint, and its convergence is good.

1. 引言

几十年来，具有良好相关性的序列设计一直是经典的研究课题 [1]。传统的相关序列的研究通常是由诸如码分多址、信道估计、安全性等应用领域设计问题引起的。近年来，该主题一直在与动态频谱相关的研究中不断发展。接入系统，除了要求低相关性外，还引入了新的频谱约束。这种系统的示例包括认知无线电 [2] [3] 和认知雷达 [4]，它们是由于频谱日益拥挤而出现的。

通常，在传统的序列设计中假定了整个频谱带的可用性，这意味着将序列能量分配给所有连续的载波。对于当今日益拥挤和碎片化的频谱，现代通信和雷达感测中已引入了许多依赖有限频谱资源的应用。认知无线电/雷达被认为是解决可用频谱稀缺的一种有前途的范例，该频谱能够感知和搜索周围环境中的频谱机会，并将其提供给次要用户进行传输。频谱机会被称为频谱限制，它是被许可用户未使用的频谱波段。在某些频谱约束中具有低频谱功率的序列称为频谱受限序列(SCS)。传统的序列通常不适用于频谱受限的系统，如 [5] 中所示。面对设计具有良好相关性和低PAPR的SCS的挑战，从数值优化的观点来看，已经有许多关于SCS设计的工作 [6] [7] [8] [9]。Liu等人 [10] 和Tsai等人 [11] 分别给出了SCSs和SNCSs和相关下界，可作为SCSs和SNCSs的最优性的度量。Hu等人 [12] 介绍了一种在任意频谱零约束下的准零相关区域(ZCZ)序列的解析构造，该序列显示零互相关和近零自相关区域性质。同时，提出了一种对准ZCZ序列的峰均比和非周期自相关进行联合优化的算法，其优点是带外功耗为零，生成速度快，计算负荷低。

作为 [12] 的后续，本文从数值优化的角度出发，提出采用模拟退火算法 [13] 对 [12] 中数值优化所得结果进行进一步优化，使所获序列在满足频谱约束的同时，还具有更低的PAPR和更低的异相非周期自相关。

2. 模拟退火算法

该算法是一种全局优化算法，其原理与序列最优相位搜索有一定的内在联系。如果我们建立起SA与序列最优相位搜索之间的联系，就可以绘制SA算法中相似的目标函数(或代价函数)。一旦建立了合适的目标函数，就可以搜索目标函数的最小值，即目标函数的最优解。

基本原理

模拟退火之所以如此命名，是因为它类似于固体的物理退火过程，在这一过程中，结晶固体被加热，然后缓慢冷却，直到它达到最规则的可能晶格构型(即，其最小晶格能态)，从而没有晶体缺陷。如果冷却进度足够慢，最终的配置会产生具有如此优异结构完整性的固体。模拟退火建立了这种类型的热力学行为和对离散优化问题的全局最小值的搜索之间的联系。此外，它提供了利用这种联系的算法手段。

在应用于离散优化问题的模拟退火算法的每次迭代中，目标函数为两个解(当前解和新选择的解)生成比较值。改进的解决方案总是被接受，而一小部分非改进的(劣质的)解决方案也被接受，希望在寻找全局最优解时避开局部最优解。接受非改进解的概率取决于温度参数，该参数通常不随算法的每次迭代而增加。这样的一个过程叫做Metropolis过程。

简单描述一下这个过程： $\left(s,f\right)$ 表示一个优化问题。i和j是这个问题的两个解， $f\left(i\right)$ 和 $f\left(j\right)$ 是成本因子。从i到j的接收标准由接收概率设计：

$P=\{\begin{array}{l}1f\left(j\right)\le f\left(i\right)\\ \exp \left(\frac{f\left(i\right)-f\left(j\right)}{T}\right)f\left(j\right)>f\left(i\right)\end{array}$ (1)

其中T是控制参数，在控制参数逐渐减小的条件下，SA算法是Metropolis过程的迭代。每个还原过程将产生L个可能的解，以形成长度为L的马尔科夫链，因此每个控制参数(温度)都有L个Metropolis过程。随着控制参数的减小，搜索范围逐渐缩小，算法将停止直到满足收敛条件。

模拟退火的关键算法特征是它提供了一种通过允许爬山移动(即恶化目标函数值的移动)来逃避局部最优的方法。随着温度参数降低到零，爬山运动发生得越来越少，与模拟算法行为的非均匀马尔可夫链相关联的解分布收敛到一种形式，其中所有概率都集中在全局最优解集上。

3. 基于模拟退火算法的SCS搜索

3.1. SC序列

假设整个频带分为N个频率槽。用 $M={\left[m_0,m_1,\cdots ,m_{N-1}\right]}^{\text{T}}$ 表示“频率槽标记向量”，给出所有N个频率槽的状态。具体而言，如果第k个频率槽有效，即可供使用，则 $m_k\left(0\le k\le N-1\right)$ 的值设置为1；否则， $m_k=0$，表示该频率槽被其他应用占用或保留。用 $\Omega$ 表示所有空频槽的位置，即 $\Omega =\left\{k|m_k=0\right\}$。在本文中， $\Omega$ 也被称为“频谱空穴约束”。如果 $\Omega$ 是非空的，即 $\left|\Omega \right|>0$，就称频谱空穴约束是非平凡的。

在本文中，用矩阵 $F_N={\left[f_{i,j}\right]}_{i,j=0}^{N-1}$ 表示N阶傅里叶变换矩阵，即：

$f_{i,j}=\frac{1}{\sqrt{N}}{\omega }_N^{-ij},0\le i,j\le N-1$ (2)

其中 ${\omega }_N=\exp \left(\sqrt{-12}\pi /N\right)$。值得注意的是 $F_N$ 是酉矩阵，即 $F_N{F}_N^H=I_N$，因此N阶离散傅里叶逆变换矩阵是 $F_N^H$

设

$\Omega =\left\{i_0,i_1,\cdots ,i_{k-1}\right\}\subset \left\{0,1,2,\cdots ,N-1\right\}$ (3)

其中 $0\le i_0<i_1<\cdots <i_{k-1}\le N-1$，为了便于表达，令

$\bar{\Omega }=\left\{0,1,2,\cdots ,N-1\right\}\backslash \Omega =\left\{j_0,j_1,\cdots ,j_{N-k-1}\right\}$ (4)

其中 $0\le j_0<j_1<\cdots <j_{N-k-1}\le N-1$。又设 $B={\left[B_0,B_1,\cdots ,B_{\mu },\cdots ,B_{N-1}\right]}^{\text{T}}$ 是长度为N的频域序列，如果 $\mu \in \Omega$，则 $B_{\mu }=0$。基于B，我们有相应的时域序列b如下：

$b={\left[b_0,b_1,\cdots ,b_{\mu },\cdots ,b_{N-1}\right]}^{\text{T}}=F_N^{H}B$ (5)

如果b的空频率槽特征为 $\Omega$，则b被称为SCS。在这个意义上， $\Omega$ 被称为b的频谱空穴约束。

3.2. 目标函数

时域序列 $\left\{b_n,n=0,2,\cdots ,N-1\right\}$ 的非周期( $C_b\left(\tau \right)$ )自相关函数(Aperiodic Auto-correlation Function, AACF)定义如下：

$C_b\left(\tau \right)={\displaystyle \underset{n=\tau +1}{\overset{N}{\sum }}{b}_n{b}_{n-\tau }^{\ast }=r{\left(-\tau \right)}^{\ast },0\le \tau \le \left(N-1\right)}$ (6)

设 ${\delta }_b$ 为时域序列b的峰值相关旁瓣(Peak Sidelobe Level, PSL)，其定义如下：

${\delta }_b=\max \left\{\left|C_b\left(\tau \right)\right||1\le \tau \le N-1\right\}$ (7)

时域序列b的峰均比(PAPR)定义如下：

$\text{PAPR}\left(b\right)=10{\log }_{10}\left(\frac{\underset{0\le n\le N-1}{\max }{\left|b_n\right|}^2}{\frac{1}{N}{\displaystyle \underset{n=0}{\overset{N-1}{\sum }}{\left|b_n\right|}^2}}\right)$ (8)

在本篇文章中，我们通过引入权重因子 $\lambda \in \left[0,1\right]$ 来控制时域序列b的PSL和PAPR的相对权重，因此我们优化问题的目标函数就可以表述为：

$E=\lambda {\delta }_b+\left(1-\lambda \right)\text{PAPR}\left(b\right)$ (9)

3.3. 接收规则和退火规则

式(1)中表示的Metropolis规则为接收规则，接收概率计算函数为：

$P=\exp \left(\frac{-\Delta E}{T_k}\right)$ (10)

退火规则为：

$T_{k+1}=\alpha \cdot T_k$ (11)

在以上设置中， $\Delta E$ 为目标函数的变化量； $T_k$ 是控制参数； $\alpha$ 为常数，在本次设计中， $\alpha$ 取0.99；k为Metropolis过程中的循环指数，为模拟退火的第k个阶段。

3.4. 借助模拟退火优化SCS

本文是在 [12] 的数值优化结果的基础上，对SCS的PAPR值和异相非周期自相关进行进一步优化。上文，我们已经建立了此次优化问题的目标函数。在对序列性能进行优化的过程中，通过序列相位值“突变”进行随机搜索。对于每次相位“突变”，计算相位变化前后的目标函数，采用(10)计算的概率接受相位变化，得到新的序列。为了成功实施SA，需要确定退火过程的几个关键参数或标准。这些是起始温度、温度的降低速率，即冷却时间表、在每个温度下的平衡条件的确定以及退火停止标准。

在温度 $T_k\left(k>0\right)$ 下，序列不断“突变”并按(10)的概率接收，直到目标函数分布达到平衡状态。然后根据(11)将温度降至 $T_{k+1}$，重复“突变”过程，直到在更新后的温度下达到新的平衡状态。在本次设计中，我们选择将每个温度 $T_k$ 下的最小目标函数值对应的序列作为下一次温度 $T_{k+1}$ 下的初始序列，这样能从一定程度上加快算法的收敛速度。

如果在连续三次温度降低中都没有“突变”的相位被接受或 $T<\epsilon$ (其中 $\epsilon$ 为最小停止温度)，则退火过程停止。选择的 $\epsilon$ 值是0.001。

算法具体流程如下所示：

步骤1：加载出 [12] 中数值优化后的结果，得到一条满足频谱空穴的频域受限序列B，其对应的时域序列具有低PAPR和低AACF；

步骤2：设置初始温度 $T_0$ 、马尔科夫链长度L (固定)、外循环终止迭代条件 $\epsilon =0.001$ 、降温速率 $\alpha$ ；

步骤3：初始化随机序列：初始序列X为一条单模频域序列，将 [12] 中通过数值优化后得到的最优频域序列B的相位作为X的初始相位；

步骤4：序列“突变”：

a) 确定“突变”位置：在序列X中随机选择两个“突变”位置；

b) 确定“突变”相位值：在 $\left[0,2\pi \right]$ 之间随机选择两个数去替代a中已选位置上的相位值，由此可得到新的单模频域序列 $X^{\prime }$，也得到了满足频域受限的新频域序列 $B^{\prime }=\left|B\right|\cdot X^{\prime }$，其中 $\left|B\right|$ 是 [12] 中通过数值优化后得到的最优频域序列B的频域幅值；

步骤5：计算“突变”前后目标函数的变化量 $\Delta E$，以Metropolis规则决定是否接收 $X^{\prime }$ ；

步骤6：重复步骤4~5；当每个温度下的迭代次数达到阈值，则根据降温速率进行降温，并选取每个温度下使得目标函数E最小的单模频域序列作为下一个温度下的初始序列；

步骤7：重复步骤4~6；达到外循环终止迭代条件时，循环结束。

4. 序列搜索结果

应用我们提出的算法并使用参考文献 [12] 实例3中相同的频谱空穴约束，即设序列长度 $N=64$，“频率槽标记向量”为 $M={\left[1_{14}，0_6，1_{20}，0_8，1_{16}\right]}^{\text{T}}$，其中 $1_m$ 和 $0_n$ 分别表示长度为m，元素全为1的行向量以及长度为n，元素全为0的行向量。因此，我们有：

$\Omega =\left\{14,15,\cdots ,19\right\}\cup \left\{40,41,\cdots ,47\right\}$ (12)

图1(a)和图1(b)分别显示了 [12] 中数值优化后和本文算法优化后的序列时域和频域幅值及其AACF。可以看出，应用SA算法， $\lambda =0$ 的SC序列显示出更平坦的时域幅度，产生0.5954分贝的PAPR和0.1349的最大异相AACF幅值(归一化)。和 [12] 中的数值优化结果相比，我们的PAPR降低了0.5946分贝，最大异相AACF幅值(归一化)降低了0.0028。

Figure 1. Comparison of two different SC sequences obtained after optimization by the algorithm in this paper and the algorithm in [12]

图1. 由本文算法和 [12] 中算法优化后得到的两条不同SC序列的比较

5. 结论

为了进一步优化 [12] 中的SC序列，使其不仅满足频谱空穴约束，而且同时具有更低的PAPR和更低的AACF，我们采用模拟退火算法在 [12] 的数值优化结果上进行进一步优化。通过本文所述算法，所得序列性能相较于 [12] 中的序列性能有较好的改善，且该序列在新兴的频谱受限系统中起着关键作用，例如认知无线电和认知雷达。除此以外，统计优化算法的全局收敛性和鲁棒性都在序列设计中得到了体现。在设计中没有挖掘算法的全部潜力，通过更慎重地选择退火控制参数可能会得到更好的结果。

致谢

本工作受到教育部春晖计划项目(No.Z2017065)的支持。从论文选题到完稿，刘子龙教授给我提出了宝贵的意见。