2. 西安中车永电电气有限公司, 陕西 西安 710016
2. CRRC Xi'an Yonge Electric Co., Ltd., Xi'an 710016, China
0 引言
交通运输是国民经济发展的坚实基础.无论是城市内公路网的控制系统,还是城际、省际间的轨道网控制系统,稳定可靠的信号识别与控制对于提高运输效率、降低运维成本、保障运输安全等都至关重要[1-2].在轨道交通运行系统中,轨道电路作为应用广泛的信号设备,可用来监督线路占用情况及将列车运行与信号显示等联系起来.轨道电路由钢轨线路、钢轨绝缘、电源、限流设备、发送/接收设备组成[3-4].其工作原理是利用电信号来控制继电器的吸合,当闭塞区间内无列车占用时,电流经由轨道流经继电器,使其励磁吸起,该区段可正常行车;当有列车占用时,低阻的列车轮对使轨道电路短路,继电器失磁落下,其它车辆禁止通行;当钢轨发生断裂时,轨道电流被切断,继电器因供电不足而失磁.为防止颠覆事故发生,即使此时线路空闲,列车仍被禁止驶入[5-6].新型数字化高压脉冲(high-voltage impulse,HVI)轨道电路与传统的轨道电路有所区别,继电器不是直接被脉冲波形驱动,而是由接收器来控制[7]. HVI轨道电路具有传输距离远、分路灵敏度高和能防止迷流干扰等优点[8].但是,当接收信号中出现干扰时,可能造成系统判断错误,导致轨道继电器错误动作.因此,HVI信号的滤波处理对轨道电路的安全运行尤为重要.文[9]采用小波变换对HVI信号识别进行了仿真研究,但该方法计算量大,不太适合实际工程应用.常用计算量较小、适合工程应用的滤波算法有最小二乘法(least squares,LS)[10]、滑动平均(moving average,MA)滤波法[11-12]等. LS滤波平滑性高,但是随着数据的增长会出现“数据饱和”现象[13];MA滤波法存在“端部效应”,不利于提高滤波后波形的信噪比[14].因此,针对HVI波形的特点,为消除数据饱和现象及端部效应,本文提出一种综合了LS及MA滤波法各自优点的改进型滤波算法.
1 高压脉冲轨道电路系统及仿真模型数字化高压脉冲轨道电路系统结构如图 1所示.每一区段内有两台主从式发送器与两台热备份接收器.发送器与接收器通过CAN总线方式将设备自检结果和/或轨道占用状态信息发送到本区段监测中心,各区段监测中心再通过GPRS通信链路将汇总信息输送至远程监控中心,供调度部门掌控管辖区域内的行车情况.正常状况下,主发送器产生3 Hz的HVI信号,自检到故障时,主发送器通过RS232串行通信方式触发从发送器工作,并同时失效自身工作状态[15].两台接收器均可使用,任何时刻,一台处于工作状态,另一台处于监测状态.如果工作状态下的接收器检测到自身硬件故障,或软件处理信息错误,则该区段监测中心在失效该接收器的同时启用另一台开始工作.
为了模拟HVI信号通过实际轨道电路的情况,本文基于传输线理论设计了如图 2所示的HVI轨道电路Simulink仿真模型.发送器产生频率为3 Hz的HVI信号,该信号通过电缆连接到发送端轨道变压器后送入钢轨.目前,HVI信号在轨道上所能传输的距离范围是900 m至1 200 m,本模型用10个级联的100 m轨道模型将传输距离设定为1 km.通过改变级联的数量可方便灵活地减少或增加轨道距离,便于模拟不同距离对信号传输的影响.当信号到达接收端后,经过轨道变压器升压后,通过电缆传送到接收器.实际应用中,作为传输线的钢轨参数会受到外界环境的影响,文中所采用的各项参数均设定为在钢轨周围土壤含水量为0.1%时的干燥条件下[16].
连接发送器/接收器与轨道变压器的电缆模型可以用图 3(a)的四端网络来表示.在HVI轨道电路模型中,发送端与接收端的电缆均选取1 km,各项参数及其取值为:1 km电缆的电阻Rc=30.6 Ω/km,电缆的电感Lc=750 mH/km,电缆的电容Cc=34.5 nF/km.根据传输线理论,钢轨同样可由分布式参数构成的电气网络模型来表示.典型参数是串联的电阻和电感及并联的电容和电导.电容值和电导值易受道碴和环境条件的影响,对称的四端口模型在这里表示轨道区间[17-18],整个区间可看作相同子区间的级联. 图 2展示的轨道模型中,每100 m长度子区间的电气网络模型如图 3(b)所示.电阻Rt和电感Lt代表了钢轨的阻抗,电导Gt和电容Ct则表示钢轨并联导纳.每100 m的电阻、电感、电导和电容四个参数分别用Rt100、Lt100、Ct100和Gt100来表示,其具体参数值为:Rt=5 Ω/km,Lt=2 mH/km,Ct=1.5 μF/km,Gt=20 mS/km.
2 最小二乘法最小二乘法是一种曲线拟合方法,是把离散的数据序列视为对一个连续的、自变量x与因变量y有明确函数关系的函数y=f(x)抽样所得来的.由于在抽样过程中存在误差,数据序列不可能都被y=f(x)的曲线所经过,但是这些离散点可以在一定程度上反映出原函数y=f(x)曲线的大致形状.实际中,虽然不能预知样本点与原函数y=f(x)的误差程度,但是可以拟合出一条曲线使得该曲线距离所有样本点的总误差和最小,这条曲线就是要得到的原函数[19-21].
设HVI信号曲线的标准函数为y=f(x),拟合曲线函数为y=S(x),单周期理想HVI信号波形如图 4所示.对单周期HVI信号进行采样,采样个数为200,将采样到的离散序列记为yi=f(xi):
(1) |
则有误差:
(2) |
以误差δ平方和E作为拟合曲线y=S(x)与HVI标准曲线y=f(x)偏离程度的标准:
(3) |
式中,ω(xi)≥0,是第i点处的权重函数,用于表示不同点在整个波形中所占的比重.在HVI波形的研究中并不区别点与点之间的比重,则有:
(4) |
为了使波形的拟合达到最优效果,即S(x)尽可能逼近f(x),E应尽可能达到最小.从式(4)可以看出,决定E值的唯一因素是拟合函数S(x).对于一些特定的数据来源,可以由经验推断出原函数f(x);也有的离散序列f(xi)可以明显地反映出其原函数类型.对于这两种情况,可以直接列出S(x)的表达式,将离散序列代入表达式中求解系数.若HVI波形不符合这两种情况,则要将它归入到多项式函数集中,设:
(5) |
只要确定了多项式最高次,解系数矩阵a使得E最小,即可得函数S(x).
在Matlab仿真平台下,使用最小二乘法多阶拟合滤除HVI信号中的噪声.考察多项式最高次分别为3、9、15时的信号处理时间及滤波后的信噪比.需要注意的是,本设计中的信号采样率为10 kHz,对每个周期的波形进行200个点的采样,理论上采样时间为20 ms,由于采取了“先采再算”的数据处理方式,本文中提到的信号处理时间(即计算时间)均不包括采样时间,而是指在2 GB内存、1.5 GHz处理器主频、Windows XP操作系统的配置下,切断网络,关闭所有其它进程,仅运行Matlab R2012b软件对信号进行处理时,整体信号处理过程完成所花费的时间,该时间是使用CPU的主频进行计算的,其Matlab程序语言:
{
Tstart=cputime;
主程序;
Tcost=cputime-Tstart;
}
其中,Tstart为程序开始的时间,cputime为运行该条语句时的当前CPU时间,最终输出的计算时间Tcost等于程序运行结束时的CPU时间减去程序开始的时间.
信噪比用来估计待检测信号与噪声功率比率,是衡量信号优劣程度的重要参数.式(6)是以分贝(dB)为单位对接收到的信号进行信噪比估计,其中Ps为待检测信号的能量均值,Pw为噪声采样点的能量均值[22].本文中所有信噪比的计算均是依据式(6)进行,其中,s(n)为待检测信号在每一个采样点下的值;w(n)为每一个采样点下噪声的值;l=1,…,L,L为采样点数.
(6) |
图 5中的实线分别是3阶、9阶和15阶拟合后的滤波效果. 3阶拟合信噪比最低,约为13.3 dB,未能还原信号波形.相比9阶拟合,15阶拟合将信噪比从26.1 dB提高至了30.6 dB,但是其计算时间却有了明显增加,不利于该算法在主频较低的嵌入式处理器上运行,因此在应用中选择9阶拟合较为合理.
3 滑动平均滤波采用滑动平均滤波法对离散振幅序列进行局部递推,对原序列有平滑作用.其核心思想是认为采样到的波形存在随机干扰ei,但是信号在一个适当小的序列空间里被视为是平稳的,在该序列空间中取均值来代替序列中的一点,以此减小随机干扰ei所造成的影响,并依次递推完成滤波[23-24].将该方法运用在实时滤波处理上时,往往用均值点取代序列空间的最后一点.本设计对HVI信号“先采后算”,采样序列的元素数量及数值均已确定,选择用均值来代替序列中间点的值,即构造一个尺寸为奇数N的窗口,将新的数据放入该窗口最后一个位置并将窗口中第1个位置的数据去除(先入先出原则),对窗口内的数据做平均处理并将结果作为窗口中间点的新值.由于滑动平均滤波算法是以采样点半径为(N-1)/2的邻域区间做均值计算,所以采样波形最开始(N-1)/2个点和最后(N-1)/2个点无法进行平均运算,该情况被称为“端部效应”,对于这些点保留其原始值.滑动平均滤波法在窗口尺寸N=5时的实现框图如图 6,其中∑/5表示对每5个点的值求平均,最终可得输出序列:
(7) |
对含有噪声的HVI信号用滑动平均滤波法进行处理. 图 7给出了在窗口尺寸N为5、11和21三种条件下含噪声信号与滤波处理后的波形对比情况.可以看出,窗口尺寸越大,滤波后波形的平滑性就越好.但是,并不是窗口尺寸越大得到的信噪比就越高.在这3种窗口尺寸下,当N=11时,拥有最高18.6 dB的信噪比,而当N=21时的信噪比为17.9 dB,因此,文中选择窗口尺寸N=11对信号做滤波处理.
4 改进型滤波算法工程上常使用端点平滑来消除滑动平均滤波法存在的“端部效应”,即:设原始波形序列为xn,m为相邻数据的个数,p和q为小于m的任意一个正整数且p+q+1=m,前端点平滑公式如式(8)所示,式(9)则为后端点平滑法[25]:
(8) |
(9) |
该算法可以理解为把序列第1个点取原始采样值,第2个点取前两个采样点的均值,依次类推.用端点平滑来消除“端部效应”的方法称为端点平滑滑动平均滤波法.
本文结合最小二乘法与滑动平均滤波的优点,提出了一种改进型滤波算法.该方法首先对原始波形序列进行滑动平均处理,之后对两端进行最小二乘拟合运算,消除“端部效应”.在找到HVI波形的峰值点后,在峰值点附近小的邻域内判断该点两侧的斜率变化规律,从而决定是否需要进行二级运算.改进型滤波算法如图 8所示.
HVI轨道电路接收器得到的波形序列为{x1,x2,x3,…,xn},序列长度为n.文中以窗口尺寸N=5为例详细介绍算法的计算过程. 图 9为对接收到的信号序列x进行滑动平均处理,最终得到序列y的处理框图,m表示该值在序列中的位置. xm是序列{3,4,…,n-2}中的任意一点(m∈[3,n-2]),经过n-4次运算,得到n-4个滤波后的点{y3,…,yn-2},而两端无法参与运算的点x1、x2、xn-1、xn则直接赋值给滤波后的序列.在一级滑动平均滤波算法得到一级滤波序列{x1,x2,y3,…,yn-2,xn-1,xn}后,对信号继续进行拟合处理.两个端部序列和与其相邻的(N-1)/2个滤波后的元素构成新的序列{x1,x2,y3,y4}和{yn-3,yn-2,xn-1,xn},对这两个序列进行拟合,用拟合后的点来更新原始端部数据.对于HVI信号,其波尾端部曲线近似于线性方程,波头端部曲线近似于二次方程,可以利用最小二乘法进行线性与二次拟合,拟合计算过程如图 10所示.
在得到拟合处理序列{z1,z2,y3,y4,…,yn-2,zn-1,zn}后开始进行平滑度判断,即判断峰值邻域内的斜率变化.对{z1,z2,y3,y4,…,yn-2,zn-1,zn}序列求取峰值并记为yk,取yk左右各5个点的邻域范围,得到离散序列{yk-5,yk-4,…,yk+4,yk+5},记Δy(m,n)=ym-yn,Δy的计算方法如图 11所示.当所有Δy的值都大于0时,即峰值点左侧邻域内序列递增,右侧邻域内序列递减,则认为波形平滑,得到最终序列{z1,z2,y3,y4,…,yn-2,zn-1,zn};否则就将{z1,z2,y3,y4,…,yn-2,zn-1,zn}视为含噪信号再次进行滑动平均滤波、端部拟合、判断平滑度的计算,直到最终得到的序列满足条件为止.
5 仿真及实验验证在图 2所示轨道模型中加入3种不同噪声后,分别对接收端得到的信号进行最小二乘法拟合、滑动平均滤波、端点平滑的滑动平均滤波及改进型滤波算法处理.并从平滑性、信噪比及计算时间三个方面对不同噪声环境下这4种算法的性能进行了对比.
首先,向轨道电路模型中加入随机噪声序列,应用4种不同滤波算法对单周期含噪信号进行处理,结果如图 12所示,滤波处理后的信号用实线表示.在随机噪声影响下,最小二乘法9阶拟合的平滑度最高,但信噪比最低(为26.2 dB)且计算时间最长;改进型滤波算法的信噪比最高(为34.9 dB)且平滑性优良,计算时间适中.
接下来,对含有周期性脉冲干扰的单周期信号进行4种滤波处理,滤波效果如图 13所示.改进型滤波算法的平滑性较高且信噪比最高,高出最小二乘法9阶拟合约2.30 dB、滑动平均滤波法约1.92 dB、端点平滑的滑动平均滤波法约1.82 dB.同时,最小二乘法的信噪比最低(为20.1 dB)且计算时间最长.改进型滤波算法的运算量低于最小二乘法而高于其它两种算法.
图 14是加入高斯白噪声后的信号滤波算法对比.虽然改进型滤波算法的计算时间不是最优的,但其平滑度理想,且信噪比为20.55 dB,高出最小二乘法9阶拟合约1.65 dB、滑动平均滤波法约1.94 dB,端点平滑的滑动平均滤波法约1.32 dB.最小二乘法的信噪比最低且计算时间最长.
在仿真的基础上,采用如图 15所示的新型数字化HVI轨道电路实验平台对改进型滤波算法的效果进行了实际验证.该实验平台由发送器、接收器、轨道变压器、钢轨和上位机构成.发送器产生HVI信号经轨道变压器降压送至钢轨,信号经钢轨传输至接收端,接收器收到HVI信号数据并处理后,将处理结果用串口发送至上位机.
随机提取接收端得到的一个周期内的采样数据,使用Matlab软件绘制出来,如图 16所示,信号的波头、波尾及中间部分均存在不平滑的数据点.
对图 16所示的信号首先进行11点的滑动平均滤波,接下来结合滑动平均滤波后的前5个点及后5个点,对未处理到的波头及波尾采样点分别进行二次拟合与线性拟合修正,拟合前与拟合后波头波尾信号的具体数据变化如表 1所示,新序列指修正前波头或波尾的数据与滑动平均滤波后数据相合产生的拟合序列,拟合效果如图 17所示.
波头幅度/mV | 波尾幅度/mV | ||||
修正前 | 新序列 | 修正后 | 修正前 | 新序列 | 修正后 |
1 686 | 1 686 | 1 679 | 1 664 | ||
1 639 | 1 639 | 1 646 | 1 674 | ||
1 620 | 1 620 | 1 621 | 1 681 | ||
1 591 | 1 591 | 1 603 | 1 688 | ||
1 603 | 1 603 | 1 592 | 1 701 | ||
1 594 | 1 688 | 1 688 | 1 705 | ||
1 594 | 1 716 | 1 716 | 1 712 | ||
1 598 | 1 738 | 1 738 | 1 720 | ||
1 618 | 1 728 | 1 728 | 1 728 | ||
1 651 | 1 726 | 1 726 | 1 735 |
图 18展示了经过一级改进型滤波算法处理后的信号,对比图 16中的原始接收信号,滤波处理后的信号平滑度得到了改善,通过峰值点左右两端斜率判断后可被正确识别,不需要进行二级运算,其信号峰值标注于图中.
为了进一步验证改进型滤波算法的可靠性,选择另一受噪声干扰较大的时间段采集到的数据进行处理.接收到的信号如图 19所示.
对该信号进行一级改进型滤波算法后可得到如图 20所示的信号,检测该信号的峰值,并根据峰值左右各5个点计算出的斜率来决定是否需要进行二级运算,峰值点及左右各5个点的信号放大图在该图中进行了标注,这11个点的具体数值由表 2给出.
采样点位置 | 幅度/mV |
26 | 2 925 |
27 | 2 973 |
28 | 3 012 |
29 | 3 063 |
30 | 3 084 |
31(一次处理峰值点) | 3 129 |
32 | 3 114 |
33 | 3 091 |
34 | 3 085 |
35 | 3 096 |
36 | 3 080 |
以峰值点为中心,分别对表 2中左右两端进行“幅度差/时间差”计算,以此判断斜率是否符合“左正右负”的标准,由数值计算得出峰值点左端符合信号幅度变化趋势,而由右端第34个与第35个采样点计算得出的信号幅度变化不符合下降趋势,不能正确判断出该信号是否正确,因此需要进行二次运算. 图 21显示了二次运算后得到的符合要求的信号.
6 结论具有高信噪比的滤波处理算法对于高压脉冲轨道电路中的信号识别非常重要.基于轨道电路模型的仿真和实际验证结果表明,本文提出的改进型滤波算法,相对于最小二乘拟合法、滑动平均滤波法、及端部平滑的滑动平均滤波法来说,其信噪比有了提高,经过滤波后的波形平滑度达到了较好的水平,并且算法复杂度适中,适用于实际工程.
[1] |
何忠贺, 王力, 张玲玉, 等.
城市交通网络正系统模型与稳态信号控制[J]. 信息与控制, 2016, 45(4): 499–505.
He Z H, Wang L, Zhang L Y, et al. Positive system model of urban traffic networks and steady-state signal control[J]. Information and Control, 2016, 45(4): 499–505. |
[2] |
李卫东, 侯丽虹.
基于卫星导航系统的高速列车定位技术研究[J]. 信息与控制, 2016, 45(4): 479–486.
Li W D, Hou L H. Research on high-speed train positioning technology based on satellite navigation system[J]. Information and Control, 2016, 45(4): 479–486. |
[3] | Thurston D F. Broken rail detection:Practical application of new technology or risk mitigation approaches[J]. IEEE Vehicular Technology Magazine, 2014, 9(9): 80–85. |
[4] | Tim D B, Kim V, Robert B. Railway track circuit fault diagnosis using recurrent neural networks[J]. IEEE Transactions on Neural Networks and Learning Systems, 2017, 28(8): 523–533. |
[5] | Allotta B, Dadamio P, Innocenti A, et al. An innovative algorithm for train detection[C]//IEEE International Instrumentation and Measurement Technology Conference (I2MTC). Piscataway, NJ, USA: IEEE, 2015: 1057-1062. |
[6] | Verbert K, Schutter B D, Babuska R. Fault diagnosis using spatial and temporal information with application to railway track circuits[J]. Engineering Applications of Artificial Intelligence, 2016, 56: 200–211. DOI:10.1016/j.engappai.2016.08.016 |
[7] |
赵军.
高压脉冲轨道电路的维护和故障处理[J]. 铁道通信信号, 2012, 48(8): 47–48.
Zhao J. Maintenance and fault treatment of high-voltage impulse track circuit[J]. Railway Signal & Communication, 2012, 48(8): 47–48. DOI:10.3969/j.issn.1000-7458.2012.08.019 |
[8] |
司开波.不对称脉冲轨道电路技术研究[D].西安: 西安建筑科技大学, 2010. Si K B. Asymmetric pulse track circuit technology study[D]. Xi'an: Xi'an University of Architecture and Technology, 2010. http://cdmd.cnki.com.cn/Article/CDMD-10703-2010134985.htm |
[9] |
黄赞武, 魏学业, 刘泽.
冲击干扰下高压不对称脉冲特征频率提取算法[J]. 铁道学报, 2012, 12(34): 71–75.
Huang Z W, Wei X Y, Liu Z. Extracting the characteristic frequency of high asymmetric pulse disturbed by impulse[J]. Journal of the China Railway Society, 2012, 12(34): 71–75. |
[10] | Aslam M S. Maximum likelihood least squares identification method for active noise control systems with autoregressive moving average noise[J]. Automatica, 2016, 69(C): 1–11. |
[11] | Banerjee P, Bansal S. Revisit of moving average technique for smoothing GNSS based timing data[J]. Mapan, 2017, 32(1): 77–85. DOI:10.1007/s12647-016-0200-6 |
[12] |
王建新, 方华威, 段晓毅, 等.
基于滑动平均的能量分析攻击研究与实现[J]. 电子与信息学报, 2017, 39(5): 1256–1260.
Wang J X, Fang H W, Duan X Y, et al. Research and implementation of power analysis based on moving average[J]. Journal of Electronics & Information Technology, 2017, 39(5): 1256–1260. |
[13] |
曾明如.
递推最小二乘法中数据饱和现象的一种消除方法[J]. 江西工业大学学报, 1990, 12(2): 62–70.
Zeng M R. One method of eliminating the "data saturation" phenomenon in the recursive least squares method[J]. Journal of Jiangxi Polytechnic University, 1990, 12(2): 62–70. |
[14] |
廖德春, 廖新浩.
几种常用滤波方法端部效应的比较[J]. 中国科学院上海天文台年刊, 2001(27): 50–56.
Liao D C, Liao X H. Comparison of deformation at ends of the results obtained by some commonly used filters[J]. Annuals of Shanghai Observatory Academia Sinica, 2001(27): 50–56. |
[15] | Victorian Rail Industry Operators Group Standards VRIOGS 012.7.6: High Voltage Impulse Track Circuits[S]. Victorian Rail Industry Operators' Group. 2009. |
[16] | Mariscotti A, Pozzobon P. Determination of the electrical parameters of railway traction lines:Calculation, measurement, and reference data[J]. IEEE Transactions on Power Delivery, 2004, 19(4): 1538–1546. DOI:10.1109/TPWRD.2004.835285 |
[17] | Chen J, Roberts C, Weston P. Fault detection and diagnosis for railway track circuits using neuro-fuzzy systems[J]. Control Engineering Practice, 2008, 16(5): 585–596. DOI:10.1016/j.conengprac.2007.06.007 |
[18] | Mariscotti A, Ruscelli M, Vanti M. Modeling of audiofrequency track circuits for validation, tuning, and conducted interference prediction[J]. IEEE Transactions on Intelligent Transportation Systems, 2010, 11(1): 52–60. DOI:10.1109/TITS.2009.2029393 |
[19] |
王泽文, 邱淑芳, 阮周生.
数值分析与算法[M]. 北京: 科学出版社, 2016: 34-35.
Wang Z W, Qiu S F, Ruan Z S. Numerical analysis and algorithm[M]. Beijing: Science Press, 2016: 34-35. |
[20] |
李庆扬, 王能超, 易大义.
数值分析[M]. 第5版. 北京: 清华大学出版社, 2008: 73-76.
Li Q Y, Wang N C, Yi D Y. Numerical analysis and algorithm[M]. 5th ed. Beijing: Tsinghua University Press, 2008: 73-76. |
[21] |
石文林, 卢先领.
TISO-OEAR模型的分解递推最小二乘辨识方法[J]. 信息与控制, 2016, 45(3): 294–300.
Shi W L, Lu X L. Decomposition-based recursive least squares algorithm for TISO-OEAR model[J]. Information and Control, 2016, 45(3): 294–300. |
[22] | Ashish B, Jyotshana K, Geetam S T, et al. A robust detector using SNR with adaptive threshold scheme in cognitive radio networks[J]. International Journal of Signal Processing, Image Processing and Pattern Recognition, 2016, 9(5): 173–186. DOI:10.14257/ijsip |
[23] | Kweon S, Shin S, Jo S, et al. Reconfigurable high-order moving average filter using inverter-based variable transconductance amplifiers[J]. IEEE Transactions on Circuits & Systems Ⅱ:Express Briefs, 2014, 61(12): 942–946. |
[24] | Jakub N, Mieczyslaw J. A novel method of blind signal detection using the distribution of the bin values of the power spectrum density and the moving average[J]. Digital Signal Processing, 2017, 66: 18–28. DOI:10.1016/j.dsp.2017.04.002 |
[25] |
裴益轩, 郭民.
滑动平均法的基本原理及应用[J]. 火炮发射与控制学报, 2003(1): 21–23.
Pei Y X, Guo M. The fundamental principle and application of sliding average method[J]. Gun Launch & Control Journal, 2003(1): 21–23. |