2. 中国人民解放军 95169部队, 南宁 530031
2. Unit 95169 of the People's Liberation Army of China, Nanning 530031, China
线性调频(Linear Frequency Modulation，LFM)是一种具有低截获性能的大时宽带宽信号，有效地解决了雷达作用距离和距离分辨率的矛盾，在现代雷达中得到了广泛的应用[1-2]。频谱弥散(Smeared Spectrum，SMSP)是由Sparrow和Cakilo[3]于2006年提出来一种有效的抗LFM的相参干扰样式。这一干扰样式利用数字射频存储器(Digital Radio Frequency Memory, DRFM)产生，并且由多个子脉冲组成，每个子脉冲具有与雷达发射信号相同的带宽，因此可以部分获得雷达脉冲压缩处理的增益。经过脉冲压缩后，该干扰样式可以生成大量近似等幅的假目标，能够有效遮盖真实目标回波，并且可以通过控制子脉冲的数目和调频斜率实现对假目标覆盖范围的控制。Sparrow和Cakilo[3]对SMSP干扰的性能进行了研究，指出SMSP可以有效对抗最大搜索目标数小于500的雷达。这一干扰样式对雷达造成了严重威胁，目前雷达系统中常用的抗干扰措施无法有效对抗，因此迫切需要研究新的抗干扰技术。
当前研究针对SMSP干扰的对抗方法可分为两大类：干扰识别和干扰抑制。SMSP干扰信号的调频斜率与雷达信号存在显著的差异，因此可以考虑用时频分析的方法来识别SMSP干扰。但是常规的线性时频分析方法(如短时傅里叶及其改进方法)时频分辨率有限[4-5]，难以准确估计干扰信号的参数，而一些改进的时频分析方法(如分数阶傅里叶变换和小波分析等)则存在计算量较大的问题[6-7]。非线性时频分析方法(如魏格纳变换、模糊函数等)具有较好的时频聚集性，但是会受到交叉项的干扰，导致参数估计误差[8-9]。文献[10]提出了一种基于匹配信号变换(Matched Signal Transform, MST)的SMSP干扰识别方法，该方法首先通过应用匹配信号变换估计干扰子脉冲的调频斜率，之后根据子脉冲调频斜率设计了一个特殊的匹配滤波器来进行SMSP干扰的识别。文献[11]根据干扰信号的模糊函数与雷达信号的差异，提取峰值个数和峰值间距作为特征来进行干扰识别。文献[12]提出了一种基于匹配傅里叶变换(Matched Fourier Transform, MFT)的SMSP干扰识别方法，该方法首先利用MFT估计子脉冲的调频斜率，之后采用两级步进搜索算法，估计出干扰信号的脉冲起始时刻，并进行干扰识别。在干扰抑制方面，文献[13]提出了一种基于复值盲源分离算法的SMSP干扰抑制方法，该方法首先通过引入复数混合矩阵将阵元间延时量转化为信号混合系数，避免了实数域盲源分离需要精确补偿阵元延时的弊端，之后应用改进的复值非高斯最大算法对接收信号进行分离，具有较好的分离性能。文献[14]同样将盲源分离算法应用于SMSP干扰抑制中，提出了一种新的基于峭度(kurtosis)的盲源分离代价函数，并应用自然梯度法求解优化问题，实现了干扰信号和回波信号的分离。但是以上2种基于盲源分离的抗干扰方法，需要应用多通道雷达，并且对于接收信号信噪比有较高的要求，这限制了该方法的应用范围。文献[15]提出了基于分数阶傅里叶变换(Fractional Fourier Transform，FRFT)滤波、原子分解(Atomic Decomposition, AD)和频域限幅的SMSP干扰对抗方法，其中FRFT滤波方法通过在变换域构建带通滤波器来实现干扰滤波，AD方法则通过选择最优的原子实现信号分解和目标回波重构，频域限幅方法则通过对接收信号的频谱进行限幅来抑制干扰信号。文献[16]提出了一种基于分段解线调的SMSP干扰抑制方法，该方法首先利用离散MFT估计干扰信号的子脉冲调频斜率，根据估计的结果构造带通滤波器对分段解线调的结果进行滤波来抑制干扰信号。文献[17]则首先根据干扰信号和回波信号经FRFT后的差异，在变换域滤除干扰信号，之后利用剩余信号与雷达回波信号频谱之间的线性关系，应用压缩感知方法精确重构回波信号，相对于文献[15]中的FRFT滤波方法，具有更好的干扰抑制效果。文献[18]根据SMSP干扰为全脉冲干扰的特性，通过在慢时间域调制信息，利用干扰信号滞后于调制信息的特性，实现对干扰的抑制。
本文根据SMSP干扰子脉冲之间的周期特性，首先利用自相关运算确定子脉冲的周期，重构各个子脉冲; 然后通过对重构子脉冲做2次赋相处理，实现了对干扰信号初相位的粗略估计; 最后根据峭度最大化原理，利用重构的干扰信号对接收信号进行对消，实现干扰抑制。
1 SMSP干扰原理及其信号模型
根据干扰产生原理，该干扰的产生可分为3步：第1步，接收天线接收雷达发射信号，经下变频转换为中频信号；第2步，在模数转换器中以采样频率fs将中频模拟信号转变为数字信号，并在移位寄存器中重复存储n次；第3步，在数模转换器中以频率nfs将移位寄存器中的数字信号转换为模拟信号，并经过上变频后形成干扰信号发射出去。
雷达发射幅度归一化LFM基带信号可表示为
|
(1) |
式中：rect(t/T)为脉冲宽度为T的矩形脉冲；K为雷达信号的调频斜率。信号的带宽为B=KT。
设目标所在距离为R，则雷达接收回波信号可表示为
|
(2) |
式中：Δt1=2R/c为目标回波的双程传播延迟时间，c为光速；φ0为回波信号的相位。
设SMSP干扰由n个子脉冲组成，第i个子脉冲可表示为
|
(3) |
式中：Δt2为SMSP干扰信号的脉冲前沿时刻；φj为干扰信号的初始相位。式(3)表明每一个干扰子脉冲具有与雷达信号相同的带宽。
因此，幅度归一化的SMSP干扰信号可表示为
|
(4) |
雷达对回波信号做去斜处理，对应的参考信号可表示为
|
(5) |
式中：Tref为接收窗的窗口宽度；tref为接收窗的起始时刻；φref为参考信号的初始相位。
经过与参考信号混频后，回波信号将转变为一个单频信号，其频率与目标距离成线性关系，而第i个干扰子脉冲则转变为
|
(6) |
式中：φji=πK{n[(i-1)T/n+Δt2]2-tref2}+φj-φref为第i个子脉冲经过去斜处理后的剩余相位。
从式(6)可以发现，经过混频后，干扰子脉冲信号变为一个LFM信号。根据去斜处理的原理，对式(6)做傅里叶变换，可得第i个干扰子脉冲所形成的假目标分布的时间范围约为
首先，根据对去斜处理结果的分析可知，由于调频斜率失配的影响，SMSP干扰必须以较大的干扰功率，才能保证对目标回波的压制效果，并且SMSP干扰波形为恒包络调制波形，因此可以根据接收信号幅度突变确定SMSP干扰的脉冲起始时刻。其次，根据SMSP干扰信号子脉冲的周期性，确定每段信号的调频斜率，即可实现干扰信号的重构。最后，引入峭度作为评价指标，通过寻找最佳的干扰抑制因子实现干扰抑制。
2.1 子脉冲周期估计根据SMSP干扰产生的原理，子脉冲的数目和子脉冲调频斜率存在对应关系。因此，只需要估计出子脉冲的周期，即可根据对应关系计算出子脉冲的调频斜率。
周期估计的方法主要有自相关法和平均幅度差函数法等[19]，这里根据接收信号的性质，采用自相关法来估计子脉冲的周期。对于任意信号x(t), 基于自相关的周期估计方法可表示为
|
(7) |
式中：Tmin和Tmax分别为子脉冲周期的下限和上限; τ为延迟时间。
雷达接收信号可表示为
|
(8) |
式中：a为干信比对应的电压增益；j(t)为干扰信号；n(t)为接收机内噪声，服从高斯分布。
根据自相关周期估计方法的定义，可将雷达接收信号的自相关函数表示为
|
(9) |
通常情况下，由于信号失配的影响，SMSP干扰所形成的假目标要想完全覆盖目标回波，必须保持较高的干信比，即a≫1。根据式(9)，对接收信号做自相关处理后，自相关函数的峰值即对应干扰信号自相关函数的峰值。因此，对接收信号不需要做额外的处理，根据自相关函数即可得到干扰子脉冲的周期信息。
根据式(7), 在对干扰进行处理之前，需要确定Tmin和Tmax的取值。由于雷达接收机无法获得关于干扰子脉冲周期的先验信息，此处根据SMSP干扰的基本原理，来确定Tmin和Tmax的取值。根据文献[3]的分析可知，子脉冲的个数至少应当为2，即周期的上限应当为T/2，而子脉冲数目为5~7时，可以获得较好的干扰效果，即子脉冲周期的下限可取T/7，为了避免出现周期估计错误，可将下限取为T/10。但是这样设定周期的范围会导致周期估计模糊。例如，当周期为T/7时，设定周期上限为T/2，会导致自相关的结果出现2个相关峰，导致周期的估计出现问题，针对这一问题，可以根据第一个峰值出现的时刻来确定周期。
假设估计出来的周期为T1，则SMSP的子脉冲个数为
|
(10) |
式中：round(·)表示朝最近方向取整数。
根据估计得到的子脉冲数目和SMSP干扰的原理，可知子脉冲的调频斜率K′=nK，由此可得到第i个重构子脉冲为
|
(11) |
式中：φjri为第i个重构子脉冲的相位，是一个未知量。
此时重构信号可表示为
|
(12) |
根据SMSP干扰的产生原理，各个子脉冲均由第一个子脉冲延迟得到，因此各个子脉冲的初始相位相同。通过自相关估计子脉冲周期，实现了干扰信号脉内调制特征的重构，但是由于重构相位未知，会导致干扰抑制的结果出现偏差。本文通过给不同子脉冲设定不同的初始相位，根据重构子脉冲与接收信号的相位差异来确定重构信号的相位。
忽略通道噪声的影响，将接收信号与重构信号共轭相乘，可得
|
(13) |
从式(13)可以发现，共轭相乘后的信号由2部分组成，第1部分为回波信号与重构信号的共轭，第2部分为一个分段函数。每一段的函数值对应着重构子脉冲与干扰子脉冲的相位差，s(t)*jr*(t)仍然为恒包络信号，因此对式(13)取实部，可得
|
(14) |
根据式(14)可知，重构子脉冲与干扰子脉冲相位差越小，则对应的实部取值越大, 对式(14)按子脉冲周期分段取均值，可得
|
(15) |
式中：E(·)为均值函数。
根据式(15), 分段取均值的结果只与重构子脉冲和干扰子脉冲的相位差有关，相位差越小，分段均值越大。根据这一关系，可以通过2次设定重构子脉冲相位来确定干扰子脉冲的相位。
步骤1 设定第i个重构子脉冲的相位为
步骤2 以步骤1所得到的相位为中心，以更小的相位间隔重新设定n个重构子脉冲的相位，并重复步骤1，即可得到更接近的相位值。
由于通过多次搜索的方法增加了相位估计的复杂性，并且噪声的存在，会对相位的估计造成影响，因此这里通过2步搜索，确定一个较为接近真实相位的重构相位。将此相位代入式(11)和式(12)，即可得到与真实干扰信号较为接近的重构信号。
2.3 基于峭度最大化的干扰抑制方法峭度是一个衡量信号非高斯性的统计量，又称为四阶累积量[20]。对于复信号，峭度的定义可表示为
|
(16) |
当kurt(y)=0时，y服从高斯分布；当kurt(y)＞0时，y服从超高斯分布；当kurt(y)＜0时，y服从亚高斯分布。并且对于2个相互独立的零均值随机变量x和y，峭度满足以下关系：
|
其中：p为任意常数。
干扰抑制的过程可表示为用雷达接收信号减去重构的干扰信号：
|
(17) |
式中：w为一个实数，表示干扰抑制因子。
由于目标回波信号和SMSP干扰信号分别由不同的物理过程产生，因此二者可视作2个相互独立的随机变量。干扰抑制前后的峭度可分别表示为
|
(18) |
|
(19) |
通常情况下，SMSP干扰所形成的假目标要想完全覆盖目标回波，必须保持较高的干信比，即a≫1。在此条件下，根据式(18)可知，雷达接收信号的峭度将主要由干扰信号决定。根据式(19)可以发现，通过选择合适的干扰抑制因子，可以使混合信号的峭度发生较大的变化。
目标回波信号和SMSP干扰信号均为亚高斯信号。根据峭度的定义，当a-w的取值与回波信号的幅度在同一数量级时，混合信号接近于高斯分布；峭度将取得最大值，最大值所对应的w即为干扰抑制因子。
整个干扰抑制方法的步骤如下：
步骤1 根据接收信号的幅度突变确定干扰信号的起始时刻。
步骤2 对接收信号做自相关处理，根据自相关的峰值时刻确定干扰子脉冲的周期。
步骤3 对重构的子脉冲设定不同的相位，根据共轭相乘结果确定干扰信号的初相。
步骤4 根据峭度最大化的原理，对干扰抑制因子的幅度进行搜索，确定最佳干扰抑制因子w。
步骤5 将干扰抑制因子w代入式(17), 实现干扰抑制。
3 仿真分析雷达信号的参数为：脉冲宽度为T=100 μs，调频斜率K=4×104 MHz/s，对应带宽为B=4 MHz；SMSP干扰的参数为：子脉冲个数为n=5, 子脉冲周期为T/n=20 μs。雷达信号和SMSP干扰信号的波形如图 2所示。
|
| 图 2 雷达信号和SMSP干扰信号波形 Fig. 2 Waveform of radar signal and SMSP jamming signal |
令干信比JSR=20 dB，信噪比SNR=0 dB时，目标回波对应的延迟时间Δt1=100 μs，干扰信号相对目标延迟Δt2-Δt1=5 μs，则雷达接收信号的时频分析结果和脉冲压缩输出如图 3所示。
|
| 图 3 接收信号时频分析和脉冲压缩结果 Fig. 3 Time-frequency analysis of received signal and pulse compression results |
从时频分析结果来看，当JSR=20 dB，即a=10时，在时频平面内，目标回波信号已经被完全覆盖。而脉冲压缩的结果显示，SMSP干扰能够在真实目标两侧形成密集的假目标，有效地遮盖了真实目标，使雷达无法实现正常的目标检测。
设延迟时间的范围为τ∈[T/10, T/2]，则接收信号的自相关函数如图 4所示。
|
| 图 4 接收信号的自相关函数 Fig. 4 Autocorrelation function of received signal |
根据自相关函数可以发现，第一个峰值对应的延迟时间为20 μs，刚好等于干扰信号的周期，表明通过自相关处理可以准确估计出干扰型号的周期，与理论分析结果是一致的。
根据估计得到的子脉冲周期，即可重构出干扰信号。不失一般性，设接收到的SMSP干扰信号初始相位φj= 
|
| 图 5 重构信号与接收信号共轭相乘结果 Fig. 5 Conjugate multiplication results between reconstructed signal and received signal |
从图 5(b)可以发现，第4段对应的均值更大，即可以认为干扰信号的相位为
设定w∈[0, 30]，针对不同的w取值，分别计算对应的峭度，则信号的峭度随w的变化规律如图 6所示。
|
| 图 6 峭度变化规律 Fig. 6 Variation curve of kurtosis |
从图 6可以发现，峭度取得最大值时的w值与a=10十分接近。实际计算得到的干扰抑制因子取值为w=10.4，将此数值代入式(17)，可得干扰抑制后的结果如图 7所示。
|
| 图 7 干扰抑制结果 Fig. 7 Jamming suppression results |
从图 7(a)可以发现，经过干扰抑制后，在时频分布中，虽然仍存在干扰能量，相对于图 3(a)来说，可以明显看到回波信号的能量已经强于干扰的能量。而由图 7(b)可以看出，对干扰抑制结果做脉冲压缩处理后，干扰所产生的假目标受到了很好的抑制，与干扰抑制前的脉冲压缩结果(见图 3(b))相比，经过干扰抑制后，干扰所产生的假目标幅度普遍下降了20 dB以上，体现了较好的干扰抑制性能。
根据LFM信号的参数，以2/B作为检测单元的宽度，待检测单元两侧保护单元数为1，参考单元数为32(即左右各16个)，采用单元平均恒虚警检测，虚警概率为Pfa=10-6，则干扰抑制前后的单元平均恒虚警检测器的检测门限如图 8所示。
|
| 图 8 干扰抑制前后的检测门限(SNR=0 dB) Fig. 8 Detection threshold before and after jamming suppression (SNR=0 dB) |
根据2.2节的分析，干扰抑制效果主要受到信噪比和相位估计误差的影响。为了定量分析该方法的性能，定义信干噪比增益G=SINR2-SINR1，其中SINR1为抑制前的信干比，SINR2为抑制后的信干噪比，均以dB为单位。图 9给出了不同信噪比和相位误差条件下，信干噪比增益的变化曲线。
|
| 图 9 信干噪比增益变化曲线 Fig. 9 Variation curves of signal-to-interference and noise ratio gain |
由图 9(a)可以看出，信干噪比增益随着信噪比的增大而增大，这是因为信噪比越大，噪声对峭度计算的影响越小，计算出来的干扰抑制因子越接近于真实的干扰幅度，干扰抑制效果就越好。从图 9(b)可以看出，相位误差越小，信干噪比增益就越大，即干扰抑制的效果越好。从图中可以发现，当相位误差Δφ∈[-0.2, 0.2]时，信干噪比增益较大，实际相位估计时，通过2次分段设置相位，将相位误差范围控制在
本文提出了一种基于信号重构和峭度最大化的SMSP干扰抑制方法，主要得到的如下结论：
1) 本文方法可以有效实现干扰的抑制，在干信比为20 dB时，经过干扰抑制和脉冲压缩，干扰所产生的假目标幅度下降了20 dB以上，并且该方法对噪声和估计的相位误差不敏感。
2) 利用延迟自相关方法可以有效估计SMSP干扰信号的子脉冲周期，以此可以准确重构SMSP干扰的脉内调制特征；对重构信号分段赋相并与接收信号做共轭相乘，对共轭相乘的结果分段取平均，即可根据分段均值较为准确地确定干扰信号的初始相位，从而实现干扰重构。
3) 根据峭度最大化的原理，本文方法可以较为精确地计算出干扰抑制因子，并且该方法仅进行一维峰值搜索，具有较快的运算速度。
本文方法并不能完全正确地计算出干扰信号初始相位和干扰抑制因子，但是利用SMSP干扰和雷达信号的失配特性，经过脉冲压缩后，仍可有效实现干扰抑制。但是本文方法未考虑多目标场景和目标多普勒频率等情况的影响，在下一步仍需要做进一步的研究。
| [1] | SCHLEHER D C. Electronic warfare in the information age[M]. Boston: Artech House, 2000: 198-199. |
| [2] | CAMP W, MAYHAN J T, DONNELL R M. Wideband radar for ballistic missile defense and range-Doppler imaging of satellites[J]. Lincoln Laboratory Journal, 2000, 12 (2): 267–280. |
| [3] | SPARROW J, CAKILO J. ECM techniques to counter pulse compression radar: 7081846[P]. 2006-07-25. |
| [4] | DURAK L, ARIKAN O. Short-time fourier transform:Two fundamental:Properties and an optimal implementation[J]. IEEE Transactions on Signal Processing, 2003, 51 (5): 1231–1242. DOI:10.1109/TSP.2003.810293 |
| [5] | ZHAO Z, SHI X Q. FM interference suppression for PRC-CW radar based on adaptive STFT and time-varying filtering[J]. Journal of Systems Engineering and Electronics, 2010, 21 (2): 219–223. DOI:10.3969/j.issn.1004-4132.2010.02.008 |
| [6] | CHEN R, WANG Y M. Universal FRFT-based algorithm for parameter estimation of chirp signals[J]. Journal of Systems Engineering and Electronics, 2012, 23 (4): 495–501. DOI:10.1109/JSEE.2012.00063 |
| [7] | MARTIS J, ACHARYA U R, MIN L C. ECG beat classification using PCA.LDA.ICA and discrete wavelet transform[J]. Biomedical Signal Processing and Control, 2013, 5 (8): 437–438. |
| [8] |
粟嘉, 陶海红, 宋大伟, 等. 基于WD和时频滑窗掩膜的SAR干扰抑制算法[J].
电子学报, 2015, 43 (12): 2345–2351.
SU J, TAO H H, SONG D W, et al. Interference suppression algorithm for SAR based on WD and sliding window masking technique in time-frequency domain[J]. Acta Electronica Sinica, 2015, 43 (12): 2345–2351. DOI:10.3969/j.issn.0372-2112.2015.12.001 (in Chinese) |
| [9] | WANG Y, ZHOU G. On the use of high-order ambiguity function for multi-component polynomial phase signals[J]. Signal Processing, 1998, 65 (2): 283–296. DOI:10.1016/S0165-1684(97)00224-7 |
| [10] | LI Y P, YING X, TANG B. SMSP jamming identification based on matched signal transform[C]//2011 International Conference on Computational Problem-Solving. Piscataway, NJ: IEEE Press, 2011: 182-185. |
| [11] |
李永平, 卢刚, 田晓, 等. 基于模糊函数的SMSP和C & I干扰识别算法[J].
航空兵器, 2011 (4): 51–54.
LI Y P, LU G, TIAN X, et al. Jamming identification algorithms of SMSP and C & I based on ambiguity function[J]. Aero Weaponry, 2011 (4): 51–54. (in Chinese) |
| [12] |
卜天翔, 董扬, 伍光新, 等. 一种污染谱干扰的识别技术研究[J].
现代雷达, 2015, 37 (3): 42–45.
BU T X, DONG Y, WU G X, et al. A study on technology of recognizing SMSP jamming[J]. Modern Radar, 2015, 37 (3): 42–45. (in Chinese) |
| [13] |
尹洪伟, 李国林, 路翠华. 一种基于复值盲分离的欺骗干扰抑制算法[J].
上海交通大学学报, 2015, 49 (10): 1564–1569.
YIN H W, LI G L, LU C H. An algorithm of deception jamming suppression based on complex-value blind source separation[J]. Journal of Shanghai Jiaotong University, 2015, 49 (10): 1564–1569. (in Chinese) |
| [14] |
李飞, 李国林, 粘朋雷. 基于盲源分离的雷达信号欺骗干扰抑制[J].
海军航空工程学院学报, 2015, 30 (5): 424–428.
LI F, LI G L, NIAN P L. Radar signal deception jamming suppressing based on blind source separation[J]. Journal of Naval Aeronautical and Astronautical University, 2015, 30 (5): 424–428. DOI:10.7682/j.issn.1673-1522.2015.05.005 (in Chinese) |
| [15] | SUN M H, TANG B. Suppression of smeared spectrum ECM signal[J]. Journal of the Chinese Institute of Engineers, 2009, 32 (3): 407–413. DOI:10.1080/02533839.2009.9671521 |
| [16] |
李永平, 唐斌. 基于分段解线条的SMSP干扰抑制方法[J].
电子信息对抗技术, 2012, 27 (1): 41–45.
LI Y P, TANG B. Suppression of SMSP jamming based on subsection de-chirp[J]. Electronic Information Warfare Technology, 2012, 27 (1): 41–45. (in Chinese) |
| [17] |
卢云龙, 李明, 曹润清, 等. 联合时频分布和压缩感知对抗频谱弥散干扰[J].
电子与信息学报, 2016, 38 (12): 3275–3281.
LU Y L, LI M, CAO R Q, et al. Jointing time-frequency distribution and compressed sensing for countering smeared spectrum jamming[J]. Journal of Electronics & Information Technology, 2016, 38 (12): 3275–3281. (in Chinese) |
| [18] | SOUMEKH M. SAR-ECCM using phase-perturbed LFM chirp signals and DRFM repeat jammer penalization[J]. IEEE Transactions on Aerospace and Electronic Systems, 2006, 42 (1): 191–205. DOI:10.1109/TAES.2006.1603414 |
| [19] |
贺思三, 赵会宁, 张永顺. 基于时频域联合滤波的中段群目标信号分离[J].
雷达学报, 2015, 4 (5): 545–551.
HE S S, ZHAO H N, ZHANG Y S. Signal separation for target group in midcourse based on time-frequency filtering[J]. Journal of Radars, 2015, 4 (5): 545–551. (in Chinese) |
| [20] |
季策, 王艳茹, 王晓宇. 一种基于标准峭度的新型复数盲分离算法[J].
东北大学学报(自然科学版), 2015, 36 (5): 614–617.
JI C, WANG Y R, WANG X Y. A new complex blind source separation algorithm based on standard kurtosis[J]. Journal of Northeastern University(Natural Science), 2015, 36 (5): 614–617. DOI:10.12068/j.issn.1005-3026.2015.05.002 (in Chinese) |