时域方法

2024-06-11

时域方法（精选10篇）

时域方法篇1

在无线通信系统中, 多径衰弱已经是一个不可避免的问题[1]。在接收端进行信道估计, 提取准确的信道信息, 是克服这个问题的有效手段。但一些常规的方法却不能达到所需要的效果, 而利用训练序列 (PN序列) 良好的循环自相关性[2], 设计有效的方法, 可以达到所需要的目的。

1 导频和数据的设计

训练序列 (PN序列) 在通信系统中获得了广泛的应用。而其中的最大长度线性移位寄存器序列 (m序列) 有非常优越的相关特性, 本文利用的就是这种序列。

长度为N的m序列有如下循环自相关性:

长度为63的m序列的循环自相关结果:

本文中, 导频序列为重复两次的一个m序列, 并在前后各增加1位, 总长度为63*2+2=128个符号。在导频序列之后紧跟长度为128*19位数据, 使导频序列的开销为5%。高速数据通道的数据帧格式如图2所示。

2 基于PN序列的时域信道估计的方法

设信道冲击响应为h (n) , 对于两径信道模型

信号x (n) 经过信道, 在接收端的信号为

所以对于导频段的m序列, 在接收端的信号为

若记m序列的自相关结果为R (n) , 经过相关器后, 可以得到

当n=0时, , 因为, 故

同理

由此, 只要经过后续的相关峰检测和参数提取, 就可以得到相应的信道参数了。

3 主要模块的FPGA实现方法

3.1 相关器

相关器是信道估计中的核心模块, 后续的所有操作都是基于相关结果是否准确的基础上的, 在本文中, 采用的是基于流水线加法器的数字相关器。

数字相关器可以看做是滤波器, 即输入数据与本地码在同步时钟的驱动下, 逐级相乘并累加, 其过程如下图所示。其中, 累加的过程是一个流水线过程。

因为数据与本地码相乘改变的只是符号, 所以用符号判决来代替乘法器, 这样就大大节省了乘法器资源, 并且达到了目的。下图是用Chipscope得到的FPGA中相关器的结果。

从图中可以看出, 相关的结果是四个相关峰, 其中第一个相关峰由于有多径的缘故, 其结果受到部分数据的影响, 准确性不如第三个相关峰, 同样的, 第四个相关峰准确性不如第二个相关峰, 所以, 在峰值搜索的时候, 搜索的是第二和第三个相关峰。

3.2 相关峰搜索和I, Q路信道参数检测

相关峰搜索的目的是得到各个相关峰的位置和时延参数, 在实际的实现中, 其实就是一个状态机。

具体过程如下:

(1) 搜索态 (复位时也在此状态) , 此时, 同时对两路输入数据进行搜索, 其中第二路数据时第一路数据延时一个m序列长度产生的。注意, 在这种情况下, 相关峰是成对出现的, 所以, 当搜索到一对相关峰 (无多径时) 或两对相关峰 (有多径时) , 转入下一个状态, 否则, 一直处于搜索态。

(2) 保护态, 此状态主要是防止偶然的数据错误造成判断错误, 在此状态下, 所进行的工作和同步态一样, 但不输出同步标志和其他相关的参数, 若连续不能搜索到相关峰, 则返回搜索态。

(3) 同步态, 此状态下不再是搜索相关峰, 而是周期性检测预期位置相关峰, 若存在, 则输出同步标志, 并输出相关参数, 否则转入保护态。

当同步后, 就可以得到各个峰值的相对位置, 并且, 在同步的同时, 将各路相关的结果存入一个双口RAM, 我们只要根据各个峰值的相对位置, 就可以得到两路的信道参数了。

4 结语

本文提出了一种利用PN序列实现时域信道估计的方法, 主要介绍了如何在FPGA中实现。由于利用了PN序列的良好的自相关性, 使得时域信号估计具有更高的精度这中方法简单有效准确[3], 另外, PN序列在进行信道估计的同时, 还可以用于频偏估计和定时同步等[4,5]。这种方法简单准确, 而且在实现的过程中占用的资源比较少, 达到了所要求的目的。

参考文献

[1]杨霄鹏, 姚昆, 史浩山等.航空信道仿真研究[J].空军工程大学学报, 2006, 7 (3) :16-17.

[2]Simon Haykin.通信系统 (第4版) [M].北京:电子工业出版社:2003.

[3]姜永权, 马逸新.一种基于PN序列的OFDM时域信道估计方法[J].汕头大学学报, 2006 (11) .

[4]Tufvesson F, Edfors O, Faulkner M.Time and frequencysynchronization for OFDM using PN-sequence preambles[C].Amsterdam:Proceedings of IEEE Vehicular Technol-ogy Conference, 1999:2203-2207.

[5]Schmidl T, Cox D.Robust frequency and timing synchron-ization for OFDM[J].IEEE Trans.on Communications, 1997, 45 (12) :1616-1621

时域方法篇2

学院

专业电子信息工程

班级

姓名

学号

时间

实验一

时域离散信号与系统分析

一、实验目的

1、熟悉连续信号经理想采样后的频谱变化关系，加深对时域采样定理的理解。

2、熟悉时域离散系统的时域特性，利用卷积方法观察分析系统的时域特性。

3、学会离散信号及系统响应的频域分析。

4、学会时域离散信号的MATLAB编程和绘图。

5、学会利用MATLAB进行时域离散系统的频率特性分析。

二、实验内容

1、序列的产生（用Matlab编程实现下列序列（数组），并用stem语句绘出杆图。（要求标注横轴、纵轴和标题）

（1）.单位脉冲序列x(n)=δ（n）

（2）.矩形序列x(n)=RN(n),N=10 201321111053 陈闽焜10.90.80.70.60.50.40.30.20.10-30-20-100n102030RN(n)10.90.80.70.60.50.40.30.20.10-50510***11053 陈闽焜

图1.1 单位脉冲序列图1.2 矩形序列

nδ(n)

(3).x(n)=e（0.8＋3j）n ； n取0－15。

***6420

图1.3 复指数序列的模图1.4 复指数序列的相角

（4）.x(n)=3cos（0.25πn＋0.3π）＋2sin（0.125πn＋0.2π）n取0－15。

201321111053 陈闽焜43210-1-2-3-4-502468n10121416-202468n10121416y(n)

05n1015

图1.4 复合正弦实数序列

（5）.把

2x *** 陈闽焜10-1-2-3-***820

图1.6 y(n)序列杆图

（7）、编一个用户自定义matlab函数，名为stepshf（n0，n1，n2）实现单位阶跃序列u[n－n1]。其中位移点数n1在起点n0和终点n2之间任意可选。自选3个入口参数产生杆图。

201321111053 陈闽焜10.90.80.70.60.50.40.30.20.***1820

M文件子程序如上所列。图1.7 自定义stepshf函数效果举例

2、采样信号及其频谱分析

（1）绘出时间信号x(t)=cos(50πt)sin(πt),时间范围t取0到2秒。

201321111053 陈闽焜10.80.60.40.2x(t)0-0.2-0.4-0.6-0.8-100.20.40.60.811.2t/second1.41.61.82

图2.1 连续信号x(t)的波形图及频谱图

（2）对于连续信号x(t)=500exp（-200nT）sin(50πnT)u（n），n=0，1，2，…，49；

分别求在T=0.5ms和T=1ms以及T=2ms三种情况下的x(t)的序列图和频谱X的幅频响应.观察是否有频谱混叠现象。

150100500-500510******53 陈闽焜25002000***00510***04550

图2.2-a 以T=0.5ms采样的序列及幅频谱图

150100500-500510******53 陈闽焜***00510***04550

图2.2-b 以T=1ms采样的序列及幅频谱图

150100500-500510******53 陈闽焜60040020000510***04550

图2.2-c 以T=2ms采样的序列及幅频谱图

3、系统的单位脉冲响应

求以下差分方程所描述的系统的单位脉冲响应h(n), 长度 0—49共50点

y(n)＋0.2 y(n－1)+0.6y(n－2)=2x(n)－3x(n－1)

8642x 10-3201321111053 陈闽焜 Amplitude0-2-4-6-8051015202530n(samples)354045

图3.1 离散系统单位脉冲响应h(n)

4、计算离散线性卷积

序列x＝[1,-1, 2, 3]与

201321111053 陈闽焜***053 陈闽焜403020010-10-2-10-3-4-2000.511.522.533.500.511.522.533.5

图5.1 系统幅频特性图5.2系统相频特性

三、回答思考题内容

（1）、在分析理想采样序列时，当选择不同采样频率获取数据，其DFT的数字频率是否一样？它们的值所对应的模拟频率是否相同？为什么？

甲醇气体的太赫兹时域波谱研究篇3

关键词：自由感应衰变；太赫兹时域光谱；甲醇

中图分类号： O 433.4 文献标志码： A doi： 10.3969/j.issn.1005-5630.2016.05.014

文章编号： 1005-5630（2016）05-0450-06

引言

甲醇气体的频谱中蕴含着大量的信息，其中涉及复杂的分子物理、量子力学等问题[1-2]。虽然甲醇分子只有6个原子构成，是一种最简单的能够受阻内旋的稍不对称分子之一，但是它的能级结构极为复杂。因此，甲醇在红外范围有丰富的光泵浦转动谱线[3]。转动谱线在物理学中十分重要，太赫兹波与分子发生共振所产生的独特谱线可作为指纹谱并用于物质识别和检测。在太空物质探测，工业生产等领域甲醇气体的检测和标定非常重要，通过太赫兹波段甲醇气体的独特指纹谱可以将它标定为探针，继而研究其他极性气体在该波段的光学特性。电磁波在共振媒体中的传播激发分子从基态跃迁到激发态，这种现象可利用光谱学的方法来检测，如转动光谱法或拉曼光谱法。此外，纯转动能级大多在太赫兹区域较为活跃，近年来最新技术的发展，为超快相干脉冲的产生和在远红外区域检测超短和相干的电磁脉冲的技术开辟了新的途径，使得气体在远红外辐射或太赫兹辐射区域可以被感应和记录[4-6]。许多分子的转动振动模式落在太赫兹波段，该波段独特的光谱特征，使太赫兹时域波谱分析技术成为非常有效的气体检测和识别手段。自由感应衰变（FID）就是出现在时域谱中的一种相干瞬时现象[7-8]，其特点是脉冲激励之后出现快速振荡和相应的呼应，这种相呼应的现象被Harde等称作相称回波[9]。FID是一个反应样品吸收和色散特性的时域表示，可以用于标定分子的指纹[10]。许多气态分子尤其是极性分子在太赫兹波段都有类似的FID特性，如CO，HCl，NH3，CH3OH等。在甲醇气体的太赫兹波谱研究中，Yu等对甲醇气体的自由感应衰变现象进行了讨论[11]。目前，国内外对甲醇气体太赫兹波段时域特征的研究较多，而对频域信息的研究较少，对在太赫兹波段甲醇气体浓度的变化情况没有相关报道。

本文在前人的基础上利用太赫兹快速扫描时域光谱系统对甲醇气体进行了检测，并通过朗博比尔定律计算了其吸收谱线，得出了不同浓度的甲醇气体在0.2～1.6 THz的吸收特性，经过数据拟合分析对自由感应衰变以及色散、离心畸变等现象进行了解释。

1 实验原理

实验装置如图1所示，实验中使用的光源是光纤式钛宝石飞秒激光器，其产生的飞秒激光中心波长为780 nm，重复频率为80 MHz，脉冲宽度为80 fs，输出功率为 160 mW。

激光从激光器发出后，经过分束片分为泵浦光和探测光两路。泵浦光分别经过慢速扫描直线电机、快速扫描直线电机、衰减片、反射镜后到达太赫兹源及光电导天线，并在天线前由透镜将能量聚焦到光电导天线的间隙上，激发太赫兹波，太赫兹波经过四个抛物面镜在样品处聚焦并最终到达太赫兹探测器上；探测光经过衰减片及反射镜到达太赫兹探测器，当探测光与泵浦光到达太赫兹探测器的光程相等时可以探测到太赫兹电场。仪器的频谱分辨率为9 GHz。

在太赫兹光路中，测试样品为一个密封的不锈钢气体样品管，两侧的窗片是2 mm厚的聚苯乙烯材料。通常，在激光传播以及太赫兹波传播过程中会有不可避免的扰动造成仪器误差，如水蒸气这种分子在太赫兹波段有强烈的吸收峰，会对太赫兹波形产生扰动。所以在实验过程中需不断地向实验环境中充入压缩干燥空气，将空气湿度控制在3%以下，尽量消除水蒸气以及其他环境因素造成的测量误差。

实验过程中，空载的气体样品管经真空分子泵抽至10 Pa左右的中真空状态，测试此时的空样品管并作为参考信号。150 μL液体甲醇样品从气管上橡胶活塞进样口注入，在相对压强的作用下甲醇分子气化充满整个气室。

2 实验结果与讨论

图2为空载以及加载甲醇气体样品时测得的太赫兹时域波形和频谱，在图2（a）中，总扫描时间是110 ps，主脉冲出现在34 ps处，聚苯乙烯窗片造成的反射峰出现在55.52 ps处，系统的动态范围大约为60 dB，1.6 THz以前的数据高于噪声线，属于可靠数据。图2（c）为太赫兹脉冲经过甲醇气体后的时域波形，可以发现除主脉冲以及窗片反射峰以外，有明显的3个相称回波出现，分别在52.10、72.8、93.6 ps 且幅度逐渐衰减。甲醇气体的太赫兹时域谱通过快速傅里叶变换计算得到功率谱，并将纵坐标尺度变换为lg坐标，清晰

可见甲醇气体在0.2～1.6 THz区域有大量的吸收峰，在后面会进一步分析。

甲醇分子属于极性分子，分子构型属于非对称分子结构，当一个亚皮秒的太赫兹脉冲在甲醇气体分子间传播时会激发分子的转动跃迁，使得分子以特定的方式转动或振动并重新辐射出自由感应衰变信号。从图2（c）中可以看到在主脉冲后面每隔一段时间会出现一个类似于主脉冲形状的时域波形，这是由于甲醇分子的转动跃迁谱线间隔大致相等，导致在时域谱上会出现时域波形的解构与重构。

图3中的4幅图分别是甲醇太赫兹时域信号的主脉冲、第1、第2以及第3的相称回波。从主脉冲波形可以看到，甲醇气体的浓度变化会使得主脉冲形成相移，从图3（b）、（c）、（d）中看出，在相称回波的波包区域甲醇气体浓度变化对幅值变化有较大影响，而对其他的时间窗口影响较小，可以确定甲醇气体在波包处相互作用形成共振模式的几率较大，而其他部分则处于杂乱的运动过程，未形成明显的振动模式。每一个重构的波形即是一个波包，代表局限在空间的某有限区域内的波动，而在其他区域的部分非常微小，波包可以分解为一组不同频率、波数、相位、波幅的正弦波；在任意时刻，这些正弦波在空间的某有限区域形成相长干涉，在其他区域会形成相消干涉。可以看出太赫兹在甲醇气体中传播时重构的时域波形会迅速衰减，这是由于电磁波在甲醇气体的传播过程中色散对信号振幅有影响。对于每一个波包，其长度都在增加，这是因为分子转动过程中引入了离心畸变而导致的。随着脉冲的激励，经过一段时间后自由感应衰变的现象逐渐消失，代表了能量的逐渐耗散，这符合能量守恒定律。

nlc202309082338

图4是对图3中的时域重构波形振幅的统计，图中的误差代表了波包截取的操作误差和多次试验的系统误差。3条实线是利用函数y=Aex/t+B对3个甲醇气体的太赫兹相称回波进行了拟合，可以看出甲醇气体在不同浓度下与太赫兹相互作用时其相称回波符合指数形式的能量衰减。

由朗博比尔定律可知，一束单色光在通过一定吸收介质后，由于介质吸收了一部分光能，透射光的强度减弱。吸收介质的浓度愈大，介质的厚度愈大，则光强度的减弱愈显著，其关系为

A=lg（Er/Es）（1）

式中：A为吸收系数；Er、Es分别为参考信号功率谱与样品信号功率谱。由此公式计算出的吸收谱线如图5 （a）所示。为了直观选取了4个频率的吸收峰进行表征，图5 （b）中实线是用y=AeR0x+B进行拟合，随着气体浓度的升高吸收峰的强度逐渐增加，其他吸收峰均符合这一特征。

经过统计，甲醇气体吸收谱线的平均间隔为48 GHz，经数据分析发现这一频域间隔与时域间隔符合关系Δt-1 ∝Δv，从而进一步解释了太赫兹脉冲激发甲醇分子跃迁，转动谱线直接反映了跃迁能级的变化，相对应的这种跃迁能级变化造成了太赫兹时域脉冲在甲醇系列气体中传播时发生自由感应衰变的现象。在频域上所有的间隔不全相等，相较于对称结构的分子如CO，HCl等具有等间隔转动谱线，非对称结构的极性分子由于刚性转轴周围原子空间质量分布的不均匀，会在转动振动时出现重心偏移，同时受到离心畸变以及多普勒展宽的影响也导致吸收谱线成非线性分布。

3 结论

利用太赫兹时域光谱技术测试了气态下甲醇在太赫兹波段的时域波形，观测到了一系列的相称回波，并利用吸收谱、分子结构以及分子动力学理论进行了解释。由于甲醇分子的非对称结构导致转动谱线大致相隔一个固定的值，但每条谱线的间隔略有差异。由于离心畸变，多普勒展宽以及色散的影响导致在时域脉冲激励后呈现出相称回波之间幅度逐渐递减，波包长度逐渐增加，能量逐渐耗散的现象。本文对利用太赫兹时域光谱技术进行气体检测、识别具有参考价值。

参考文献：

[1] BORDEN A，BARKER E F.The infra-red absorption spectrum of methyl alcohol[J].The Journal of Chemical Physics，1938，6（9）：553-563.

[2] KOEHLER J S，DENNISON D M.Hindered rotation in methyl alcohol[J].Physical Review，1940，57（11）：1006-1021.

[3] MORUZZI G，WINNEWISSER B P，WINNEWISSER M B，et al.Microwave，infrared，and laser transitions of methanol atlas of assigned lines from 0 to 1258 cm-1[M].CRC Press，1995.

[4] 陈麟，高春梅，徐嘉明，等.铝质亚波长圆孔阵列与牛眼结构的太赫兹波传输特性[J].光学仪器，2013，35（6）：1-6.

[5] 彭滟，周云燕，陈麟.太赫兹波的古依相移研究[J].光学仪器，2014，36（6）：489-491.

[6] 苏云云.液态一元醇的太赫兹时域光谱研究[J].光学仪器，2014，36（6）：499-503.

[7] ABRAGAM A.The principles of nuclear magnetism[M].Oxford：Oxford University Press，1963.

[8] GARSIDE B.Optical resonance and two-level atoms[J].Optica Acta：International Journal of optics，1975，22（12）：1041-1042.

[9] HARDE H，KEIDING S，GRISCHKOWSKY D.THz commensurate echoes：periodic rephasing of molecular transitions in free-induction decay[J].Physical Review Letters，1991，66（14）：1834-1837.

[10] MITTLEMAN D M，JACOBSEN R H，NEELAMANI R，et al.Gas sensing using terahertz time-domain spectroscopy[J].Applied Physics B：Lasers and Optics，1998，67（3）：379-390.

[11] YU B L，YANG Y，ZENG F，et al.Direct observation of coherent rotational excitation，dephasing and depopulation of methanol and its isotopes using THz pulse radiation[J].Applied Physics Letters，2005，86（10）：101108.

时域方法篇4

复杂背景下低信噪比的点目标检测一直是红外图像检测与处理研究的热点和难点。由于目标仅占很少几个像素,而且随着现代隐身技术的发展,使得军事目标辐射强度越来越小,背景杂波强度往往比传感器噪声和被检测目标信号还大,图像信噪比极低。因此,提高目标的信噪比成为小目标检测的关键。目前,提高目标信噪比的方法主要是通过抑制杂波背景来实现的,如自适应背景预测[1]、中值滤波器[2]、形态学方法[3]、小波变换[4]等。这些方法在背景具有缓变结构时效果很好,但背景起伏较大时具有较高的虚警率。此外,此类方法都是单帧处理,并不能结合弱目标的运动信息抑制强杂波背景。

鉴于此,以提高弱目标的信噪比为主要目的,利用多帧条件下目标的运动特性,本文提出了一种新的空时域滤波方法。

1 空时域滤波框架

本文给出的空时域滤波结构如图1所示。首先,以第一帧为参考帧图像,对各帧图像进行全局运动估计,补偿由于摄像头位置变化引起的背景偏移;然后,对运动补偿后的各帧图像进行方差加权信息熵滤波,对滤波后图像采取双向隔帧差分的时域滤波;最后通过检测差分图像中的“凸包”检测小目标,抑制背景和噪声。图1中,Fk表示第k帧图像,F′(k)表示经过方差加权信息熵滤波后的第k帧图像,B1为后向差分图像,B2为前向差分图像。

2 基于方差加权信息熵的空域滤波

从数值意义上来说,红外图像的方差可以在统计意义上表述此图像中整体像素与其均值的偏离程度,但它无法提供任何关于图像灰度分布的空间信息。为了更好地描述图像空域特性,本文将方差和信息熵结合起来描述红外图像的空域特性,实验结果验证了此方法的有效性。

2.1 一种有效的红外图像特性描述指标:加权信息熵[5]

信息熵是描述图像中所含信息量的有效途径[6]。对于一幅具有n个灰度级的红外图像,设ps是图像中灰度值s出现的概率,S是由此红外图像中所含像素值所组成的集合,则S的概率空间可表述为

其中:0≤ps≤(1s=1,0,...,n-)1;∑s=0n-1ps=1。

对于我们所讨论的具有256个灰度级的红外小目标图像,其图像信息熵可定义为

规定当ps=0时,pslogps=0。图像的信息熵能有效地表达出图像中所包含的平均信息量,但它也忽视了图像中灰度信息本身的重要性,无法真实地反映主观上对背景复杂程度所做的判断。由此,为了体现高灰度部分对于小目标检测的贡献,将式(2)中概率ps对应的灰度值s作为权重,即可将式(1)修改为

上式可称为红外图像的加权信息熵(Weighted Information Entropy,WIE),由于它兼顾了图像的灰度信息和人眼的主观判断,所以也就提供了一种定量描述不同红外背景特性的有效途径。

2.2 融合加权信息熵和图像方差信息的改进指标:方差加权信息熵

基于加权信息熵的红外背景特性描述指标在复杂背景下的小目标检测中表现出了良好性能。但是,当待处理的红外图像的像素数量增大时,基于加权信息熵的红外背景特性描述指标的鲁棒性将受到影响。红外图像的加权信息熵能将图像的灰度信息和人们关于背景中灰度成分分布的主观判断结合起来,但红外图像的方差仍然是一种反映该图像中灰度值整体变化程度的重要指标。结合信息熵和方差,将式(4)修改为

上式中,s为红外图像的整体灰度均值。将式(5)命名为图像的方差加权信息熵(Variance WIE)。根据对方差加权信息熵的定义,设计方差加权信息熵滤波器。通过计算红外小目标图像中各像素点处的局部方差加权信息熵值,可以得到一个局部方差加权信息熵图像,称此方法为方差加权信息熵滤波。

3 基于凸包检测的前后向隔帧差分时域滤波

对单帧图像进行方差加权信息熵滤波能够大幅度提高目标的信噪比,这只是将目标检测出来的一个先决条件。结合目标的多帧运动特性,考虑算法的实时性,本文提出一种基于凸包检测的前后向隔帧差分时域滤波方法。因此对图像序列作时域差分一方面可以大大减弱起伏背景[7],另一方面也可以突出目标的运动信息,此外还可以去除一些固定的噪声点[8]。时域差分另一个突出的优势就是运算简单,算法效率较高。

3.1 前后向隔帧差分

常用时域差分方法存在以下弊端:一是当目标运动缓慢时,相邻两帧目标存在重叠,差分图像提取目标运动信息不充分;二是当目标运动较快时,差分图像可以提取到运动目标,但是检测出的目标比真实的小目标大很多,对目标的定位也不准确。针对以上缺点,本文采取双向隔帧差分相乘的时域差分方法。

参照图1,后向差分图像为B1=f(x,y,k)-f(x,y,k-n);前向差分图像为B2=f(x,y,k)-f(x,y,k+n)。n与目标速度估计有关,当目标速度较大时可适当增大,一般情况下取1~3。当前帧的运动目标增强图像为

通过前后向隔帧差分相乘后,当前帧运动小目标位置得到有效增强。

3.2 凸包检测抑制背景和噪声

小目标在红外图像中呈现亮斑特性,因此在通过上述隔帧前向、后向差分以后,差分图像B1和B2(未取绝对值)中沿目标运动的方向必然会形成一种紧邻的亮暗差分对。这种紧邻的明暗差分对反映了小目标短时运动连续性质,而红外背景中的随机噪声其位置特性都是不固定的,产生这种差分对的可能性很小,即使在强噪声水平下产生了这种明暗差分对,也可利用前后向差分图像中差分对的位置关于较亮的目标位置对称的性质,剔除强噪声。称这种前后向差分图像中,关于亮点位置对称的紧邻的明暗差分对为前后向差分图像的“凸包”特性,采取如下的基于差分图像“凸包”检测的噪声抑制方法:

1)在差分图象Dtime(x,y,k)每个像素的局部邻域内搜索其极大值,邻域大小(2nvmax+1)×(2nvmax+1),vmax为其可能的两帧间最大速度,n为隔帧数。

2)以搜索到的局部极大值位置(x0,y0)为中心,在图像B1和B2中对应的八邻域方向,nvmax范围内搜索局部极小值,在B1和B2中搜索得到的位置分别为(x1,y1),(x2,y2),如果(y1-y0)/(x1-x0)=(y2-y0)/(x2-x0),且B1和B2中极小值数值之差不超过其较小值的30%,可认为局部极大值点(x0,y0)满足凸包特性,属于小目标点。保留Dtime(x,y,k)中以(x0,y0)为中心的局部区域内的像素值,认为该区域为小目标区域,在此区域不再进行凸包性质检测。

3)若某点像素不满足2)中所述的凸包特性,则认为该点为噪声点,予以剔除。

4 实验结果与分析

初步估计目标的运动速度在1.5 pixel/f,选择红外视频序列中的第15帧,第17帧和第19帧图像进行实验,即隔帧差分的隔帧数为1帧。图像分辨力为240×320,Variance WIE滤波采用3×3的窗口,实验结果如图2所示。图2(a1)∼(a4)分别为第15帧原始图像,中值滤波后图像,标准差滤波后图像,本文Variance WIE滤波后图像,图2(b1)为第15帧和第17帧未进行运动补偿的差分图像,图2(b2)为对第15帧和第17帧进行运动补偿、Variance WIE滤波后的后向差分图像,图2(b3)为对第17帧和第19帧进行运动补偿、Variance WIE滤波后的前向差分图像,图2(b4)为图2(b2)和图2(b3)的差分增强图像,图2(b5)为对图2(b4)的凸包检测图像。

实验运行环境:CPU双核2.4 GHz;内存2 G;仿真软件Matlab7.01。产生(b5)图像耗时0.476 s。比较图2(a1)∼(a4)的具体滤波效果可知,方差加权信息熵滤波对于红外图像中小目标对象信号的强度的提升比标准差滤波和中值滤波要明显的多,如表1所示。从图2(b2)∼(b5)可以看出,经过本文滤波后的差分图像,目标信号强度很高,再经过对差分增强图像的凸包检测,小目标被凸现出来。

5 结论

针对红外弱小目标图像的特点,本文提出了一种新的空时域滤波方法检测小目标。空域滤波采用方差加权信息熵滤波方法,时域滤波采用实时性较好的前后向隔帧差分方法,最后对前后向隔帧差分的增强图像进行凸包检测,提取出小目标。实验结果表明该滤波方法能极大地提高目标的信噪比,实时性好。但是,文中主观估计目标运动速度,这会影响隔帧数的选择,从而影响检测概率。下一步将重点研究如何估计目标的运动速度,从而提高检测概率。

参考文献

[1]CHENG De-jie,LI Xiao-feng,LI Zai-ming.Motion Analysis and Morphological Approach for Optical Aerial Conflict Alarm in Trafffic Alert and Collision Avoidance System[C]//The 6th International Conference on ITS Telecommunications,Chengdu,China,June 21,2006:314-317.

[2]Barnett J.Statistical analysis of median subtraction filtering with application to point target detection in infrared background[J].SPIE(S0277-786X),1989,1050:10-18.

[3]YE Bin,PENG Jia-xiong.Application of order morphology filtering on detection of small target and point target[J].SPIE(S0277-786X),2001,4554:94-99.

[4]Davidson G,Griffiths H D.Wavelet detection scheme for small target in sea clutter[J].IEEE Electronics Letters(S0013-5194),2002,38(19):1128-1130.

[5]YANG Lei,ZHOU Yong,YANG Jie,et al.WIE Based Infrared Images Processing[J].IEEE Electronics Letters(S0013-5194),2006,42(15):857-859.

[6]Brink A D,Pendock N E.Minimum cross-entropy threshold selection[J].Pattern Recognition(S0031-3203),1996,29(1):179-188.

[7]LIA Cheng-fang,ZHANG Fei,SHI Li-na.Effective Target Tracking Algorithm from IR Image Sequences[J].Optical Engineering(S0091-3286),2001,45(3):1021-1029.

液态一元醇的太赫兹时域光谱研究篇5

关键词：太赫兹光谱；一元醇；微孔金属片

中图分类号： TN 29文献标志码： Adoi： 10.3969/j.issn.10055630.2014.06.007

引言太赫兹波（Terahertz wave）是介于微波和红外辐射之间，振荡频率为0.1～10 THz（1 THz=1012 Hz）的一种电磁辐射。近几十年超快飞秒激光技术的迅速发展，为THz脉冲的产生提供了稳定、可靠的激光光源，使得太赫兹技术作为一种新的光谱手段在实际应用方面取得了很大的进展。许多分子之间弱的相互作用如氢键、范德华力、生物大分子的偶极旋转和振动跃迁都处于THz频带范围，这使得太赫兹技术在生物、医药检测、化学材料分析等方面有着广泛的应用[15]。醇类物质是常见的一类重要有机物质，在食品、化工、生物制药等领域有广泛的应用，同时由于醇分子的自身特点也常被用来作为研究氢键网络的理想模型。1996年Kindt等[6]研究了极性液体包括水和一些醇以及氨水等在远红外波段（2～50 cm-1）的介电常数等介电性质，并通过德拜模型进一步提取其内部的微观弛豫过程，提示太赫兹时域光谱技术可用于研究液态物质微观结构和动力学。2005年Woods等[7]利用远红外光谱检测甲醇液体，反映了其太赫兹频段的低频振动与液体内部氢键结构以及分子间相互作用有关。2007年，Jepsen等[8]采用反射式太赫兹时域光谱技术检测了一些含醇液体样品在0.1～1 THz的介电信息，提示该技术可在一系列商业饮品中鉴别出其所含醇、糖的浓度。2010年Yomogida等[9]报道了一系列不同结构一元醇在太赫兹波段的介电响应及温度效应，发现其介电弛豫包括三个部分：介电弛豫过程的高频阶段、1.2～1.5 THz处的振动模式、低频分子间的振动模式。其中宽的振动模式峰不随温度变化，其介电性质与一元醇的分子结构有着密切的关系。随着太赫兹技术应用研究的不断发展，结合传感器件对液体样品检测的需求十分突出。最近，基于表面等离子体效应的超材料结构已作为传感器应用于液体样品检测方面。2009年，Tian等[10]采用微孔金属片器件通过表面共振增强效应成功区分出甲醇中氢的同位素。2012年，Hasebe等[11]研究表明特定结构的带孔金属网栅可有效区分微量的单、双链DNA分子。2009年，Mendis等[12]在可激发出TE（transverse electric wave）模式的平行板波导中镶嵌矩形槽形成高灵敏度的光学共振器，并通过监测折射率变化有效区分出微量的辛烷等一系列烷类物质。本文以甲醇、乙醇、正丙醇一元醇体系为研究模型，利用太赫兹时域光谱技术并结合特定结构的微孔金属器件，对样品进行检测，提取一系列如吸收系数、折射率、介电常数等光学参数，并对经微孔阵列金属器件产生的共振吸收进行分析和讨论。

1实验实验所用的甲醇（纯度>99.5%）、乙醇（纯度>99.7%）、正丙醇（纯度>99.5%）均购于国药科技有限公司。实验采用透射式时域太赫兹光谱装置（Terahertz time domain spectroscopy，THzTDS）[1315]，该系统采用光电导天线发射和接收THz波，Toptic飞秒激光器用于产生和检测太赫兹波。激光器的中心波长为780 nm，输出功率为150 mW，脉冲宽度为90 fs，重复频率为80 MHz。系统的发射器、接收器均为低温生长的砷化镓光电导天线，系统在0.4 THz处的信噪比超过105∶1，太赫兹的光斑半径约为5 mm。实验中，通入干燥空气，使湿度降至3%以下，以减小环境中水蒸气对太赫兹信号的吸收，同时环境温度维持在23.0～25.0 ℃。实验中利用液体样品厚度差的方式处理数据，通过选用不同厚度（400 μm、600 μm）的石英样品池可以减小样品容器对实验造成的误差[16]。2实验结果与分析利用太赫兹时域光谱装置，可以获得样品的时域光谱信息，再经傅里叶变换后进行样品光学参数的提取，其参数表达式为α（ω）=2[lnI（d1，ω）-lnI（d2，ω）]d2-鉴于表面等离子体效应的超材料结构在液体检测方面的敏感性，本文利用特定结构的微孔金属器件，对不同结构的一元醇液体样品进行检测，所采用的微孔金属片结构如图2所示。该微孔金属片为铝制材料，长宽为100 mm×100 mm、厚度250 μm，分布有超过5 000个亚波长小孔，其中小孔直径为700 μm，六角晶格周期是1 130 μm[1921]。其工作原理为：当入射波照射到金属孔表面时，在金属介质界面可以激发出表面等离子体激元（SPPs）效应，透射谱中将会出现共振峰。该现象可以采用Drude模型对SPPs效应做进一步分析，SPPs的波矢和频率之间存在以下关系[22]kSPPs=ωcεdεmεd+εm=ωcεm-ε2mεd+εm（5）式中，εm、εd分别为金属与样品的介电常数，ω为电磁波的角频率，c为自由空间的光速。当周围样品的介电常数发生改变，则对应的金属样品界面的SPPs波数发生改变，共振峰发生漂移。本实验中不同的醇产生不同的共振峰，并且随着介电常数的增加，共振频点发生蓝移。实验中，通过装有液体的石英样品池后的太赫兹信号作为参考信号Ein（ω），通过装有液体及金属片的石英比色皿后的信号作为样品信号Eout（ω），透过率t（ω）=Eout（ω）/Ein（ω）。图3为微孔金属片存在情况下所获得的甲醇、乙醇及正丙醇的透过谱，其中产生的共振峰频点分别为0.133 THz、0.138 THz、0.143 THz。由此可见，采用微孔金属片区分样品优点主要是可以产生共振峰，通过样品介电常数的不同可以看出其对应共振频点发生蓝移。同时，和常规太赫兹时域透射光谱所探测到的无特征吸收光谱相比，共振峰的体现使检测更显直观，所使用的样品量也大大减少，此外，在数据处理方面也更简单方便。因此该类器件作为传感器在生物化学领域有着潜在的应用价值。

3结论本文利用太赫兹时域光谱技术结合特定结构的微孔金属器件，对一元醇液体样品进行了检测，从中提取出了吸收系数、折射率、介电常数等一系列重要的光学参数。研究表明，太赫兹时域光谱对液体样品的极性和结构特点非常敏感，结合微结构器件在液态化学和生物样品检测方面具有潜在的应用价值。

参考文献：

[1]FERGUSON B，ZHANG X C.Materials for terahertz science and technology[J].Nature Materials，2002，1（1）：2633.

[2]WOODWARD R M，COLE B E，WALLACE V P，et al.Terahertz pulse imaging in reflection geometry of human skin cancer and skin tissue[J].Physics in Medicine and Biology，2002，47（21）：3853.

[3]DU S Q，LI H，XIE L，et al.Vibrational frequencies of antidiabetic drug studied by terahertz timedomain spectroscopy[J].Applied Physics Letters，2012，100（14）：143702.

[4]YAMAMOTO K，KABIR M H，TOMINAGA K.Terahertz timedomain spectroscopy of sulfurcontaining biomolecules[J].JOSA B，2005，22（11）：24172426.

[5]朱亦鸣，高春梅，陈麟.基于平行板波导的双槽谐振腔的特性研究[J].光学仪器，2014，36（4）：323327.

[6]KINDT J T，SCHMUTTENMAER C A.Farinfrared dielectric properties of polar liquids probed by femtosecond terahertz pulse spectroscopy[J].The Journal of Physical Chemistry，1996，100（24）：1037310379.

[7]WOODS K N，WIEDEMANN H.The influence of chain dynamics on the farinfrared spectrum of liquid methanol[J].The Journal of Chemical Physics，2005，123（13）：134506.

[8]JEPSEN P U，MOLLER U，MERBOLD H.Investigation of aqueous alcohol and sugar solutions with reflection terahertz timedomain spectroscopy[J].Optics Express，2007，15（22）：1471714737.

[9]YOMOGIDA Y，SATO Y，NOZAKI R，et al.Comparative dielectric study of monohydric alcohols with terahertz timedomain spectroscopy[J].Journal of Molecular Structure，2010，981（1）：173178.

[10]TIAN Z，HAN J G，LU X C，et al.Surface plasmon enhanced terahertz spectroscopic distinguishing between isotopes[J].Chemical Physics Letters，2009，475（1/3）：132134.

[11]HASEBE T，KAWABE S，MATSUI H，et al.Metallic meshbased terahertz biosensing of singleand doublestranded DNA[J].Journal of Applied Physics，2012，112（9）：09470210947027.

[12]MENDIS R，ASTLEY V，LIU J B，et al.Terahertz microfluidic sensor based on a parallelplate waveguide resonant cavity[J].Applied Physics Letters，2009，95（17）：171113.

[13]ULLAH R，LI H，ZHU Y M.Terahertz and FTIR spectroscopy of ‘Bisphenol A’[J].Journal of Molecular Structure，2014，1059：255259.

[14]CHEN L，GAO C M，XU J M，et al.Observation of electromagnetically induced transparencylike transmission in terahertz asymmetric waveguidecavities systems[J].Optics Letters，2013，38（9）：13791381.

[15]CHEN L，XU J M，GAO C M，et al.Manipulating terahertz electromagnetic induced transparency through parallel plate waveguide cavities[J].Applied Physics Letters，2013，103（25）：25110512511054.

[16]刘丹，吴胜伟，马明旺，等.水和二氧六环水溶液的THz波段介电性质[J].核技术，2010，33（9）：688692.

[17]INDIRA T，PARTHIPAN G，ASWATHAMAN H，et al.Dipole moment studies of complexes of alcohols with ethyl bromide[J].Journal of Molecular Liquids，2009，150（1/3）：2224.

[18]BERTERO N M，TRASARTI A F，APESTEGUA C R，et al.Solvent effect in the liquidphase hydrogenation of acetophenone over Ni/SiO2：A comprehensive study of the phenomenon[J].Applied Catalysis A：General，2011，394（1/2）：228238.

[19]XU J M，CHEN L，XIE L，et al.Effect of boundary condition and periodical extension on transmission characteristics of terahertz filters with periodical hole array structure fabricated on aluminum slab[J].Plasmonics，2013，8（3）：12931297.

[20]CHEN L，ZHU Y M，ZANG X F，et al.Mode splitting transmission effect of surface wave excitation through a metal hole array[J].Light：Science & Applications，2013，2（3）：e60.

时域方法篇6

关键词：智能车,路径跟踪,控制策略,滚动时域估计

0 引言

随着智能控制技术的发展, 其研究成果广泛应用于诸如家庭、医疗、救灾等许多领域。在这些应用中, 有些需要人工操作机器进行控制, 还有些由机器自动完成作业任务。智能技术的不断成熟及应用减轻了工人的劳动强度, 在危险或不便由人员参与的恶劣环境中更可以替代工人工作。以智能车为依托的智能控制技术已经成为该领域内的一个研究热点[1,2,3]。其中, 智能车路径跟踪技术又是智能控制研究的一个重要分支。许多文献设计实现了小车对预设道路的准确跟踪[4,5], 而对循迹跟踪精度及其快速性的研究更是得到了广泛的关注[6]。

模糊策略等许多优化控制方法被用来优化提高小车循迹的准确性。为了在保障准确性的基础上优化小车的循迹速度, 实现在线处理并优化控制量, 本研究基于模型预测的滚动优化原理, 以智能车的路径跟踪问题作为研究对象, 阐述目标道路信息的获取、处理、识别过程。

具有机器视觉的智能车路径跟踪控制系统能够更全面地获取道路信息, 利用各种有效图像处理算法对路径进行识别, 提高路径识别的准确性。由于该系统拓宽了智能车的视野, 更能对更远的路径提前做好判断, 提高了路径跟踪的平均速度。目标道路信息由安装在小车顶部的CCD摄像头获取[7,8], 研究者基于视觉图像来进行分析识别, 并在此基础上运用有效的控制算法对智能小车进行控制, 对小车的转向和速度进行调整, 使得智能小车能准确快速地对道路进行跟踪。

基于滚动优化原理的滚动时域估计方法随着模型预测控制研究的不断深入得到广泛关注[9]。它的基本思想是将估计问题转化为固定时长的优化问题, 简化了计算, 使得在线处理系统约束问题成为可能。

本研究首先建立智能小车的非线性约束动态数学模型, 然后基于滚动优化时域估计方法提出智能小车循迹控制策略。

1 智能小车的数学模型

本研究首先建立智能小车的动态数学模型, 如图1所示。

X, Y—小车的位置;θ—小车的方向角;v (t) , ω (t) —小车的即时前进速度和转动角速度

根据以上参数, 该小车的数学模型为:

为了方便对小车进行数字化智能控制, 笔者利用前向差分近似计算方法, 并取采样时间为T, 可得到小车系统式 (1) 的离散数学模型为:

接下来分析系统的状态空间表达式。

由于:

应用三角函数加法定理, 并引入参数:η1=sinθ, η2=cosθ, α1 (ω) =cos (Tω) , α2 (w) =sin (Tω) , 由离散系统式 (2) , 可得:

进一步取系统变量:

可得到智能小车离散模型的状态空间表达式为:

其中:

在系统的状态方程式 (7) 中, 输入控制量为小车前进速度v和转动角速度ω。显然, 行列式A与速度v和角速度ω有关。若假设小车前进速度v恒定, 则只有转动角速度为控制变量, 此时行列式A仅与角速度ω有关, 状态空间模型可用下式表示:

2 滚动时域估计控制策略

卡尔曼滤波常被用来设计估计器, 但是, 在存在约束条件的目标道路跟踪问题中, 卡尔曼滤波方法往往难以达到期望精度, 会导致估计结果与实际情况不符[10]。

本研究采用滚动时域优化估计方法来解决小车对目标道路的跟踪优化问题。这种方法具有处理约束和滚动优化的特点, 能够将估计问题转化为带约束的优化问题, 从而使问题得以解决[11]。另外, 滚动时域估计算法还避免了计算量随时间不断增大的缺点, 使得在线处理优化问题成为可能。

根据滚动时域估计方法的基本思想, 系统的预定跟踪目标设定为r, 预测区间长度为N, 系统输出的预测输出值为y (k+1) , y (k+2) , ⋯, y (k+N) 。本研究将当前时刻开始的N个控制输入量记为u (k) , u (k+1) , ⋯, u (k+N-1) 。这些值同时计算得到, 但是只有当前控制输入u (k) 被用于系统的实际控制。每一步的控制量的计算方法都一样。计算输入控制量的目标是为了使预测区间内的预测输出值接近目标值。这种优化方法被称为滚动时域估计方法。

坐标系的选择和滤波算法是目标跟踪研究领域的两个主要问题。本研究仅考虑智能小车实现平面跟踪的情形。在时刻k, 系统 (1) 的输出, 即小车位置为:

跟踪的目标路线位置为:

式中:Xwl, Ywl—广域坐标系中目标道路位置的坐标。

在离散模型式 (9) 中, 系统在时刻k, k+1, ⋯, k+N-1的对应输入量为wi, i=1, 2, ⋯, N, 相应的预测输出为:

因此, 系统状态方程为:

其中, 小车转动角速度ω需要满足约束:

其中, 约束条件的参数D, δ由小车转动电机的性能决定。

为了使得输出预测值与目标值之间的误差达到最小, 实现智能小车对目标道路的精确跟踪, 本研究引入目标函数:

本研究假设系统的估计初值x (0) 是零均值的正态分布变量, 求解具有约束条件 (14) 的最小优化问题:

即可得到系统优化输入列:

进而实现小车与目标道路之间误差的最小化。

3 目标道路的获取方法

智能小车追踪的目标是目标道路的位置, 下面对目标道路的获取方法进行说明。

该试验中, 目标道路的位置不是预先给定的。研究者利用固定在小车顶部的摄像机采集道路信息, 并以摄像机采集到的图像的其左上角为原点建立屏幕坐标系Os-xsys, 此时, 目标道路的位置坐标为 (xsl , ysl) 。再换个角度, 从三维空间中俯瞰的角度来分析要跟踪的目标道路。

行走路线的摄像头坐标如图2所示。

本研究以摄像头所在位置为坐标原点, 建立三维坐标系Oc-XcYcZc, 在这个坐标系中, 原来的路线坐标 (xsl , ysl) 被重新定义为 (Xcl , Ycl , Zcl) , 二者之间的关系为:

其中:

式中:Xoff, Yoff, Zoff, WD, WS—摄像机硬件参数;d—视觉误差值, 这个值越大表示被拍摄物体离镜头越近, 相反地, 它越小就表示被拍摄物体离镜头越远;Ycl—摄像机与路线地面的高度差。

该试验中, 假定行车道路是平坦的, 即Ycl=c为给定常数, 那么, 道路位置为:

在式 (20) 中, 目标道路的坐标表示其在三维坐标系Oc-XcYcZc中位置。因为该坐标系是以小车为原点建立的, 目标道路的位置并非其绝对坐标, 而是以小车为基准的相对位置坐标。接下来, 笔者设法将坐标 (Xcl , Ycl ) 转换成广域坐标系中的道路位置 (Xwl , Ywl ) 。广域坐标系是指三维空间中以任意点为坐标原点的坐标系。该试验中, 笔者取小车的出发点为原点建立广域坐标系Ow-XwYw, 如图2所示。那么, 二者之间的变换公式为:

式中:Xwl , Ywl—广域坐标系中目标道路的位置;X , Y—由式 (9) 求得的小车位置坐标;θv—两个坐标系坐标轴之间的夹角。从摄像头坐标到广域坐标的变换图如图3所示。

式 (11) 中的目标道路位置坐标值通过式 (21) 在小车行进过程中依次计算获取。

4试验结果

试验中, 用方向盘控制小车的转角w。方向盘转角u和小车转动角度的关系可由试验数据得到, 即:ω=0.275 4 u。

假设小车匀速前进, 设定速度v=150 mm/s。受约束条件式 (14) 的限制, 方向盘的最大转动角度选取u1=-30°和u2=30°两个值。此时, 系统行列式为:

目标道路坐标值:

试验小车跟踪其前方约25 cm处开始的目标道路, 其结果如图4所示。

输出预测值, 即小车预定位置坐标由式 (22, 23) , 取N=3, 按照式 (12) 计算得到。

由最小化目标函数式 (17) 计算小车方向角的控制输入量, 结果如图5所示。

从图4中可以看出, 小车的行进轨迹与识别路径相一致, 实现了对预定路线的跟踪控制。需要注意的是图4中小车轨迹的起始位置和停车位置和识别路径有误差, 这是由于摄像机安装在车顶, 其采集得到的目标道路必然超前于小车本身的位置。

5 结束语

针对搭载摄像头的智能小车, 本研究建立了小车的俯视平面动态模型, 并基于滚动时域控制方法提出了对目标道路跟踪的优化控制策略。通过摄像头采集道路信息具有更大的信息获取量, 拓宽了智能小车的视野, 提高了其路径识别的准确性。而基于滚动时域估计的预测控制方法, 将估计问题转化为固定时长的优化问题, 也提高了其对目标道路预测的快速性。

试验结果表明, 该方案对小车实现智能目标道路跟踪是合理可行的, 并且具有较高的跟踪精度。但是试验中的目标道路仅仅是一段平滑圆弧, 并未对急转弯等极端情形下小车的跟踪性能进行考察, 这将是接下来的研究内容;另外, 小车行走过程中的平滑度问题也需在今后的研究中进一步深入讨论。

参考文献

[1]卢韶芳, 刘大维.自主式移动机器人导航研究现状及其相关技术[J].农业机械学报, 2002, 1 (2) :112-116.

[2]尉成果, 赖欢, 王太江.轮式移动机器人轨迹跟踪控制算法的研究[J].机电工程, 2012, 29 (6) :730-732.

[3]LAMIRAUX F, LAUMON J P.Smooth motion planning forcar-like vehicles[J].IEEE Transactions on Robotics And Automation, 2001, 17 (4) :498-502.

[4]董雷刚, 崔晓微, 张丹.智能循迹电动小车的设计[J].大庆师范学院学报, 2010, 30 (6) :8-10.

[5]龙鸿儒, 何汉武, 胡兆勇, 等.虚拟环境下基于非水平路面的智能自主汽车运动研究[J].机电工程技术, 2011, 40 (8) :13-15.

[6]许晓飞, 夏培容.智能车速度的模糊控制最优循迹策略研究[J].北京信息科技大学学报, 2009, 24 (3) :84-86.

[7]刘崧, 戚小平, 钟双英.CCD摄像机原理及应用[J].中国有线电视, 2005 (14) :1417-1419.

[8]张云洲, 师恩义, 吴成东.基于CCD的智能车导航系统研究[J].东北大学学报:自然科学版, 2009, 30 (2) :162-165.

[9]王曌, 刘志远, 裴润.一种滚动时域估计方法在目标跟踪中的应用[J].机器人, 2001, 23 (7) :594-597.

[10]WANG Xiao-fei, ZOU Zao-jian, WANG Yang, et al.Path following control of underactuated ships based on unscent ed Kalman filter[J].Journal of Shanghai Jiaotong Uni versity (Science) , 2010, 15 (1) :108-113.

时域方法篇7

基于OFDM系统的信道估计方法众多,G3-PLC利用前导序列估计信道信息,通常可采用传统的LS(Least Squares)、LMMSE(Linear Minimum Mean Square Error)算法,或基于DFT(Discrete Fourier Transform)变换的估计算法。LS算法简单,但在低信噪比时性能不理想。LMMSE算法具有良好的性能,但统计信道相关性和矩阵求逆的复杂度较高,通常采用奇异值分解(Singular Value Decomposition,SVD)降低求逆复杂度[3]。基于变换域的DFT算法将LS频域估计值变换到时域,通过时域置零实现降噪,该算法复杂度低于LMMSE算法,但其性能对阈值的选取较为敏感[4]。

本文针对G3-PLC系统,通过预置电力线信道最大时延估计其时域特性,具有降噪效果,但相应的时频变换矩阵不满秩[5],无法求逆。在预置最大信道时延后,Simeone[5]提出采用子空间跟踪算法迭代获得时延空间正交基,将传统LS频域估计值在子空间上投影,达到信噪分离的效果,该方法需要多次矩阵相乘,实现复杂度较高。本文通过修正变换矩阵,解决了采用LS准则存在的矩阵不可逆问题,并取得了良好的性能。同时进一步提出一种无需信道统计特性的时域LMMSE算法,大幅降低了复杂度。仿真结果表明,两种改进算法复杂度均较低,性能良好,对系统实现具有参考意义。

1 系统模型

1.1 电力线信道

目前的文献中,窄带电力线通信的信道模型尚未形成统一标准。根据Zimmermann[6]提出的电力线信道多径模型,其传递函数可表示为

其中,N为路径总数;gi、di表示路径i的加权系数和长度;α0、α1表示衰减参数,由传输线特性决定;β是指数衰减因子,通常介于0.5~1之间;εr是电力线介电常数;c0表示光速。

本文在进行系统算法仿真中,即采用了上述信道模型,并叠加以高斯白噪声。

1.2 G3-PLC系统

G3-PLC[7]采用OFDM调制技术,图1为系统物理层收发框图。系统采样率为400 k Hz,FFT点数为256,支持35.9~90.6 k Hz频段通信,对应256个载波中36个子载波(编号23~58)。G3-PLC的前导序列由8个SYNCP符号和1.5个SYNCM符号组成,系统利用以上9.5个OFDM符号进行信道估计。

在G3-PLC中,已知36个幅值相同、相位各异的调制信息P(k),23≤k≤58,经系统OFDM调制得训练序列p(n),0≤n≤255。p(n)和-p(n)分别对应前导的SYNP、SYNM符号。训练序列p(n)经过信道,接收端FFT解调,得到OFDM符号Y(k)

式中,W(k)是功率为σw2的噪声;H(k)即电力线信道频域特性。

2 信道估计算法

2.1 算法分析

设H(k)对应有限长时域冲激响应h(n),两者关系为

假设信道最大时延为L,则上式用向量表示为

其中,HK对应K个可用子载波;Fsub表示变换矩阵

其中,WNkn=e-j2πkn/N。

根据式容易得到信道频域LS估计值,相应的估计误差为

其中,Δ可等效为功率σw2的噪声。将式代入得

上式即预置时延后的时频变换关系。相对于系统的时间分辨率,这种时延间隔过于紧密,导致Fsub不满秩[5],是一个病态系数阵,无法直接应用LS准则估计时域h。

2.2 修正的时域LS算法

为有效应用LS准则,本文采取了修正变换矩阵Fsub的方法。设Fsub的奇异值分解形式为

其中,U1和V1是K×K和(L+1)×(L+1)的酉矩阵,不妨设K≥L+1,则

其中,σ1≥σ2≥…≥σr≥σL+1>0,因Fsub呈病态,奇异阵S包含若干接近于零的奇异值σi。当σi较小时,直接采用LS准则的均方误差将较大,虽不是无偏的,但并非为最优解[8]。本文采用截断奇异值法[8],设置截断参数T,将<T的奇异值直接用0代替,得修正的奇异阵Sr

其中,σr满足≥T。相关文献给出了L曲线法、极小化均方误差法等计算T,这里Fsub确定,可通过仿真分析选取最优截断参数T。

综上,经过修正,相应的时域估计值为

进一步将时域估计值变换到频域,即

将式(8)和式(11)代入上式可化解得

其中,Ir×r表示r阶单位阵;U1∑U1H具有共轭对称性。

这种修正的时域LS算法中,Fsub由可用子载波k的分布和信道最大时延L确定,所以实际可预先计算出矩阵系数U1∑U1H存储在本地,所需存储空间为(K2+K)/2。接收端先计算出频域LS估计值,再和本地矩阵相乘,即得算法改进值,其中包含K2次复乘,K2-K次复加运算,其复杂度较低。

2.3 简化的时域LMMSE算法

为进一步提高估计性能,对式采用LMMSE准则估计h

其中

整理得

其中,Rhh表示信道时域自相关阵。实际中,统计Rhh和矩阵求逆的复杂度均较高。

假设各径相互独立,则Rhh为对角阵。而本文将Rhh设为单位阵,则上式可简化为

当调制子载波等间距时,FsubFHsub为共轭对称的Toeplitz矩阵,可采用Levinson-Durbin递归算法降低矩阵求逆的复杂度。这里采用SVD[3]进一步化解,则有

其中,U2是酉矩阵;Λ为对角阵,对角元素满足λ1≥λ2≥…≥λK。

综上,简化的信道时域估计值可化解为

进一步将时域估计值变换到频域,并化解

图2是该改进算法的实现流程图。这种简化的时域LMMSE算法,将时域自相关阵简化为单位阵,即在信道统计特性未知的情况下,假定各径统计均匀分布。根据图2流程,简化算法包含2K2+K次复乘,2K2-2K次复加运算。传统的LMMSE算法,不考虑计算信道二次相关性的复杂度,仅矩阵求逆复杂度为o(K3)。因此相比之下,简化算法复杂度大幅降低。

3 仿真结果及分析

针对G3-PLC系统,仿真中采用BPSK调制,假设信道最大时延L为22个采样点。为了最优化修正时域LS算法,在不同截断参数T下,对系统均方误差(MSE,Mean Square Error)性能进行了仿真对比,如图3所示。比较各曲线,当T=1时,修正算法性能总优于T<1;当T>1时,低信噪比时性能较好,但在高信噪比时性能变差。综合整体性能与稳定性分析,选取T=1作为修正算法的最优截断参数。

基于以上条件,将改进算法和传统的频域LS算法、LMMSE算法、DFT算法进行了性能对比,包括MSE和BER(Bit Error Rate)性能,仿真结果如图4和图5所示。其中,LMMSE算法的协方差矩阵由真实的H(f)计算,仿真得到理想LMMSE值。而基于变换域的DFT算法,仿真中截取时域两端窗长为最大时延的3/4和1/4,中间补零[9]。

图4表示各算法的MSE性能,传统的频域LS算法性能最差,DFT算法性能介于LS和LMMSE算法之间。由图4可看出,通过预置信道时延,修正的时域LS算法相比频域LS算法性能提高了7 d B。对于G3-PLC系统,可用子载波数K为36,所以修正算法只要额外1 296次复乘,1 260次复加,即可求得改进值,从而在较低的复杂度下大幅提高了估计性能。对于简化的LMMSE算法,算法复杂度为2 628次复乘,2 520次复加,而传统LMMSE算法复杂度为46 656次复乘。根据图4所示,相比理想LMMSE估计值,由于信道二次相关性被简化近似,简化算法仍存在一定的性能损失,但复杂度大幅降低。进一步,对比两种改进算法,由于简化算法中考虑了噪声的影响,所以其性能优于修正的时域LS算法,特别在低信噪比时,MSE性能提高1~2 d B。

图5给出了G3-PLC系统下,不同信道估计算法的BER性能。相比传统的频域LS算法,修正的时域LS算法性能提升约3 d B,而简化算法与理想LMMSE算法之间的差距<1 d B。仿真表明,两种改进算法性能良好。

4 结束语

时域方法篇8

时域有限差分方法[finite-difference time-domain (FDTD) method]自1966年由Yee[1]提出来后, 一直是计算电磁学常用的方法之一, 并广泛应用于电磁散射、天线的分析与设计、雷达截口的计算等电子工业与国防工业。但是, 由于FDTD方法是条件稳定的, 即在二维情形下, 时间和空间步长分别为Δt, Δx, Δy;必须满足Courant-Friedrichs-Lewy (CFL) 稳定条件: cΔt≤[1/ (Δx) 2+1/ (Δy) 2]-1/2。其中c是在介质中的光速, ε, μ是介电常数和磁导率。为了解决上述问题, 由文献[2,3]提出了交替方向隐式时域有限差分方法 (ADI-FDTD) 方法, 并用Fourier方法证明了这种格式是无条件稳定的。近年来[4], 将分裂算子方法与FDTD方法结合提出新颖而且简单的分裂算子的时域有限差分方法 (S-FDTD) , 与二维ADI-FDTD相比都是无条件稳定的二阶格式, 但是S-FDTD格式采用算子分裂技术, 所以格式简单, 计算时间短, 并且在模拟一种散射问题时, S-FDTD比ADI-FDTD更精确。麦克斯韦方程的高阶时域有限差分方法已有较多的研究工作[4,5,6,7,8,9], 但是这些方法都没有使用分裂算子技术。

本文将高阶差分与S-FDTD方法相结合提出一种新的方法, 即在分裂后的方程基础上对空间采取四阶中心差分, 从而得到分裂 (splitting) 的高阶 (high-order) 的时域有限差分 (FDTD) 格式SHO-FDTDⅠ;并在此格式的基础上通过增加扰动项减小分裂误差给出它的修正格式SHO-FDTDⅡ, 给出了具体的计算步骤, 然后用Fourier方法证明了这两种格式是无条件稳定的, 并给出这两种格式的数值弥散关系式和数值弥散误差的计算。与条件稳定的高阶FDTD相比, 新方法是无条件稳定的;与高阶的ADI-FDTD相比, 新格式更简单, 易于编程和实现。

1 二维麦克斯韦方程的SHO-FDTD格式

考虑以下二维横向电波:

$\frac{\partial E_{x}}{\partial t} = \frac{1}{ε} \frac{\partial Η_{z}}{\partial y}$ (1.1)

$\frac{\partial E_{y}}{\partial t} = - \frac{1}{ε} \frac{\partial Η_{z}}{\partial x}$ (1.2)

$\frac{\partial Η_{z}}{\partial z} = \frac{1}{μ} (\frac{\partial E_{x}}{\partial y} - \frac{\partial E_{y}}{\partial x})$ (1.3)

其中 $\vec{E} = (E_{x} (x, y, t), E_{y} (x, y, t))$ 表示电场, $Η_{z} = Η_{z} (x, y, t)$ 表示磁场, ε和μ分别是介电常数和磁导率, Ω=[0, a]×[0, b], t ∈ (0, T]。其边界满足理想导体边界条件: $E_{x} (x, 0, t) = E_{x} (x, b, t) = E_{y} (0, y, t) = E_{y} (a, y, t) = 0$ 。初始条件: Ex0=Ex0 (x, y) , Ey0=Ey0 (x, y) , Hz0=Hz0 (x, y) 。

为了书写的方便, 我们仅讨论常系数的情形, 这里所用的方法容易推广到变系数的情形。对空间区域Ψ以及时间区域[0, T]采取如下剖分:

其中Δx和Δy分别是沿x轴方向和沿y轴方向的空间离散步长, Δt是时间步长, I, J, N为整数。

对于任意一个网格函数F (t, x, y) , 引入以下记号:

F $_{α, β}^{m}$ =F (mΔt, αΔx, βΔy) , $δ_{t} F_{α, β}^{m} = \frac{F_{α, β}^{m + \frac{1}{2}} - F_{α, β}^{m - \frac{1}{2}}}{Δ t},$

$δ_{x}^{2} F_{α, β}^{m} = \frac{1}{Δ x} (\frac{1}{24} F_{α - \frac{3}{2}, β}^{m} - \frac{9}{8} F_{α - \frac{1}{2}, β}^{m} + \frac{9}{8} F_{α + \frac{1}{2}, β}^{m} - \frac{1}{24} F_{α + \frac{3}{2}, β}^{m})$ ,

$δ_{y}^{2} F_{α, β}^{m} = \frac{1}{Δ y} (\frac{1}{24} F_{α, β - \frac{3}{2}}^{m} - \frac{9}{8} F_{α, β - \frac{1}{2}}^{m} + \frac{9}{8} F_{α, β + \frac{1}{2}}^{m} - \frac{1}{24} F_{α, β + \frac{3}{2}}^{m})$ 。

用E $_{v_{α, β}}^{m}$ (v=x, y) 表示电场Ev (tm, xα, yβ) 的近似值。H $_{z_{α, β}}^{m}$ 表示磁场Hz (tm, xα, yβ) 的近似值。将分裂算子的时域有限差分方法[9]与高阶差分方法[4]相结合, 提出如下格式:

第一步:

$\frac{E_{y_{i, j + \frac{1}{2}}}^{n + 1} - E_{y_{i, j + \frac{1}{2}}}^{n}}{Δ t} = \frac{1}{2 ε} δ_{x}^{2} = (Η_{z_{i, j + \frac{1}{2}}}^{n + \frac{1}{2}} + Η_{z_{i, j + \frac{1}{2}}}^{n})$ (1.4a)

$\frac{Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}} - Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n}}{Δ t} = - \frac{1}{2 μ} δ_{x}^{2} (E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} + E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n})$ (1.4b)

第二步:

$\frac{E_{x_{i + \frac{1}{2}, j}}^{n + 1} - E_{x_{i + \frac{1}{2}, j}}^{n}}{Δ t} = \frac{1}{2 ε} δ_{x}^{2} = (Η_{z_{i + \frac{1}{2}, j}}^{n + \frac{1}{2}} + Η_{z_{i + \frac{1}{2}, j}}^{n})$ (1.5a)

$\frac{Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} - Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}}}{Δ t} = \frac{1}{2 μ} δ_{y}^{2} (E_{x_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} + E_{x_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n})$ (1.5b)

此格式称为分裂算子 (splitting) 的高阶 (high-order) 时域有限差分方法, 记为SHO-FDTDⅠ。其中, 边界条件为: $E_{x_{i + \frac{1}{2}, 0}}^{n} = E_{x_{i + \frac{1}{2}, j}}^{n} = E_{y_{0, j + \frac{1}{2}}}^{n} = E_{y_{i, j + \frac{1}{2}}}^{n} = 0$ 。初始值:

求解步骤为:先解第一步, 由式 (1.4b) 得:

$Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}} = Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n} - \frac{Δ t}{2 μ} δ_{x}^{2} (E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} + E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n})$ (1.6)

将式 (1.6) 代入式 (1.4a) 整理得:

显然这是一个七对角的方程组, 系数矩阵元素都是常数而且在每一个时间层上都不变。因此可由一些求解线性方程组的方法如:超松弛迭代法, 共轭梯度法等, 求出 ${E_{y_{i, j + \frac{1}{2}}}^{n + 1}}$ 。然后, 将 $E_{y_{i, j + \frac{1}{2}}}^{n + 1}$ 代入式 (1.6) 显式求出 $Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}}$ 。式 (1.5) 的计算与第一步相似, 但要用到第一步算出的 $Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}}$ 。方程 (1.7) 中的第一个方程和最后一个方程的系数需要调整, 可根据文献[5]中的方法进行, 这里略去。

检查格式SHO-FDTDⅠ的截断误差发现关于时间它的精度不高。为了提高精度, 引入一个扰动项, 得到修正格式SHO-FDTDⅡ如下:

第一步:

第二步与SHO-FDTDⅠ中的第二步相同.

这种格式的初边值条件与式 (1.4) 、式 (1.5) 的初边值条件完全相同, 求解方法也与前一格式相同。为了了解格式SHO-FDTDⅡ的逼近精度, 由式 (1.5a) 、式 (1.5b) 、式 (1.8a) 、式 (1.8b) 消去中间项 $Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}}$ , 可得到SHO-FDTDⅡ的等价格式:

$δ_{t} E_{x_{i + \frac{1}{2}, j}}^{n + \frac{1}{2}} = \frac{1}{2 ε} δ_{y}^{2} (Η_{z_{i + \frac{1}{2}, j}}^{n + 1} + Η_{z_{i + \frac{1}{2}, j}}^{n}) (1.9 a)$

$\begin{array}{l} δ_{t} E_{y_{i, j + \frac{1}{2}}}^{n + \frac{1}{2}} = - \frac{1}{2 ε} δ_{x}^{2} (Η_{z_{i, j + \frac{1}{2}}}^{n + \frac{1}{2}} + Η_{z_{i, j + \frac{1}{2}}}^{n}) + \\ \frac{Δ t}{4 μ ε} δ_{x}^{2} δ_{y}^{2} (E_{x_{i, j + \frac{1}{2}}}^{n + 1} - E_{x_{i, j + \frac{1}{2}}}^{n}) (1.9 b) \end{array}$

$\begin{array}{l} δ_{t} Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}} = \frac{1}{2 μ} [δ_{y}^{2} (E_{x_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} + E_{x_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n}) - \\ δ_{x}^{2} (E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} + E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n})] (1.10) \end{array}$

格式SHO-FDTDⅠ的等价格式与式 (1.9a) 、式 (1.9b) 、式 (1.10) 相似, 只要将式 (1.9) 最右端后一项中的“-”改为“+”。由此可以看出, SHO-FDTDⅡ的摄动误差 (二阶) 比SHO-FDTDⅠ的摄动误差 (二阶) 要高一阶, SHO-FDTDⅡ是关于时间2阶, 关于空间是4阶的式 (2.4) 格式, SHO-FDTDⅠ则是式 (1.4) 格式。

2 稳定性分析与数值弥散分析

用Fourier方法分析这两种格式的稳定性。假定格式的差分解具有下列形式。

$\vec{E}_{α, β}^{n} = \vec{E}_{0} ξ^{n} e^{- i (k_{x} α Δ x + k_{y} β Δ y)}$ ;

H $_{z_{α, β}}^{n}$ =Hz0ξne-i (kxαΔx+kyβΔy) 。

其中 $i = \sqrt{- 1}$ , $\vec{E}_{0} = (E_{x_{0}}, E_{y_{0}})^{Τ}$ 和Hz0是振幅, $\vec{k} = (k_{x}, k_{y})$ 是波矢量, $| k | = \sqrt{k_{x}^{2} + k_{y}^{2}}$ , ξ是增长因子。下面将E $_{α, β}^{n}$ , H $_{z_{α, β}}^{n}$ 代入式 (1.5a) 、式 (1.9) 、式 (1.10) , 化简得:

$(ξ - 1) E_{x_{0}} - i \frac{Δ t}{ε} (ξ + 1) b_{y} Η_{z_{0}} = 0$ (2.1)

$\begin{array}{l} (ξ - 1) E_{y_{0}} + \frac{(Δ t)^{2}}{ε μ} (ξ - 1) a_{x} b_{y} E_{x_{0}} + \\ i \frac{Δ t}{ε} (ξ + 1) a_{x} Η_{z_{0}} = 0 (2.2) \end{array}$

$\begin{array}{l} (ξ - 1) Η_{z_{0}} - i \frac{Δ t}{μ} (ξ + 1) b_{y} E_{x_{0}} + \\ i \frac{Δ t}{μ} (ξ + 1) a_{x} E_{y_{0}} = 0 (2.3) \end{array}$

其中 $a_{x} = \frac{(\frac{1}{12} \sin \frac{3}{2} k_{x} Δ x - \frac{9}{4} \sin \frac{1}{2} k_{x} Δ x)}{2 Δ x}, b_{y} = \frac{(\frac{1}{12} \sin \frac{3}{2} k_{y} Δ y - \frac{9}{4} \sin \frac{1}{2} k_{y} Δ y)}{2 Δ y}$ 。由于向量 (Ex0, Ey0, Hz0) 是非零向量, 所以关于Ex0, Ey0, Hz0的方程组的系数矩阵的行列式的值为零。整理得:

(ξ-1) (d0ξ2+2d1ξ+d0) 2=0 (2.4)

式 (2.4) 中

$d_{0} = 1 + \frac{(Δ t)^{2}}{ε μ} (a_{x}^{2} + b_{y}^{2}) + \frac{(Δ t)^{4}}{ε^{2} μ^{2}} (a_{x} b_{y})^{2};$

$d_{1} = - 1 + \frac{(Δ t)^{2}}{ε μ} (a_{x}^{2} + b_{y}^{2}) + \frac{(Δ t)^{4}}{ε^{2} μ^{2}} (a_{x} b_{y})^{2}$ 。

解方程式 (2.4) 得: $ξ_{1} = 1, ξ_{2} = d_{0}^{- 1} (- d_{1} + i \sqrt{d_{0}^{2} - d_{1}^{2}}) 。 ξ_{3} = d_{0}^{- 1} (- d_{1} - i \sqrt{d_{0}^{2} - d_{1}^{2}})$ 。显然方程三个根的模都是1, 所以格式SHO-FDTDⅡ是无条件稳定的。对于格式SHO-FDTDⅠ, 同理可得:

c3ξ3+c2ξ2+c1ξ+c0=0 (2.5)

式 (2.5) 中

$\begin{array}{l} c_{3} = 1 + \frac{(Δ t)^{2}}{ε μ} (a_{x}^{2} + b_{y}^{2}) + \frac{(Δ t)^{4}}{ε^{2} μ^{2}} (a_{x} b_{y})^{2}, \\ c_{2} = - 3 + \frac{(Δ t)^{2}}{ε μ} (a_{x}^{2} + b_{y}^{2}) + 3 \frac{(Δ t)^{4}}{ε^{2} μ^{2}} (a_{x} b_{y})^{2}; \end{array}$

$\begin{array}{l} c_{1} = 3 - \frac{(Δ t)^{2}}{ε μ} (a_{x}^{2} + b_{y}^{2}) + 3 \frac{(Δ t)^{4}}{ε^{2} μ^{2}} (a_{x} b_{y})^{2}, \\ c_{0} = - 1 - \frac{(Δ t)^{2}}{ε μ} (a_{x}^{2} + b_{y}^{2}) + \frac{(Δ t)^{4}}{ε^{2} μ^{2}} (a_{x} b_{y})^{2} 。 \end{array}$

方程 (2.5) 的根的表达式非常复杂 (为了简单, 这里省略) , 通过对根的表达式的研究发现 $| ξ | = 1 + Ο (Δ t)$ , 因此格式SHO-FDTDⅠ是耗散的无条件稳定的。为了直观的感受 $| ξ |$ , 下面给出它在不同情况下的变化曲线图。设kx=kcos φ。ky=ksin φ。其中k, φ。是向量k的圆柱坐标, 再设 $Ν_{λ} = \frac{λ}{h}, ω = c k, ω$ 是频率, λ是波长, Δx=Δy=h, 则Nλ是一个波长内的节点个数。

令 $S = \frac{c Δ t}{h}$ , 由CFL的定义可知

$c Δ t \sqrt{\frac{1}{Δ x^{2}} + \frac{1}{Δ y^{2}}} = \frac{c Δ t}{h} \sqrt{2} = \sqrt{2} S$ 。因此, 不妨称S为CFL数。由方程 (2.4) 或式 (2.5) 及其系数的表达式可知, 增长因子ξ是关于S, φ, Nλ的函数。例如, 可推出 $\sin \frac{3}{2} k_{x} Δ x = \sin (\frac{3 π}{Ν_{λ}} \cos φ)$ 。方程 (2.5) 的根的表达式比较复杂。下面我们用Matlab求出它的近似解, 并画出不同情况下它的根的变化曲线。

图1是在S=0.35, Nλ=10情况下给出了根 $| ξ |$ 随角度φ的变化曲线。上边的曲线是方程复数根的模, 下边的曲线是实数根的绝对值。很容易看出复数根的模比1大, 但是 $| ξ | = 1 + Ο (Δ t)$ , 其中, $Δ t = S \times \frac{1}{Ν_{λ}} = 0.035$ 。从图1中显然可以看出曲线位于: y=1±0.035两条直线之间, 这就表明SHO-FDTDⅠ格式是无条件稳定的。

图2和图3分别是在S=1.5, φ=35°和Nλ=10, φ=65°情况下, 方程 (2.5) 的根的模分别随着Nλ和S的变化情况, 从这两个图中同样可以看出SHO-FDTD格式是无条件稳定的。

2.2 数值弥散关系

假设麦克斯韦方程的差分解为: $\vec{E}_{α, β}^{n} =$ $\vec{E}_{0} e^{i (k_{x} α Δ x + k_{y} β Δ y - ω n Δ t)}$ , H $_{z_{α, β}}^{n}$ =Hz0ei (kxαΔx+kyβΔy-ωnΔt) , 其中, ω, kx, ky, $\vec{E}_{0}$ 和Hz0与2.1节的符号相同。将E $_{α, β}^{n}$ , H $_{z_{α, β}}^{n}$ 代入SH-FDTDⅡ的等价格式 (1.5a) 、式 (1.9) 、式 (1.10) , 类似于2.1节对增长因子ξ的推导, 可得:

$\begin{array}{l} \sin^{2} (\frac{1}{2} ω Δ t) = (c Δ t)^{2} [a_{x}^{2} + b_{y}^{2} + (c Δ t)^{2} a_{x}^{2} b_{y}^{2}] \\ \cos^{2} (\frac{1}{2} ω Δ t) (2.6) \end{array}$

式 (2.6) 中 $c^{2} = \frac{1}{ε μ}$ , 式 (2.6) 就是SHO-FDTDⅡ的数值弥散关系式。

同理可得SHO-FDTDⅠ格式的数值弥散关系式:

$\begin{array}{l} \sin^{2} (\frac{1}{2} ω Δ t) = (c Δ t)^{2} [a_{x}^{2} + b_{y}^{2} + (c Δ t)^{2} a_{x}^{2} b_{y}^{2} \times \\ \cos (\frac{1}{2} ω Δ t) \sin^{- 1} (\frac{1}{2} ω Δ t)] \times \\ \cos^{2} (\frac{1}{2} ω Δ t) (2.7) \end{array}$

由于 $\lim_{x \to 0} a_{x} = \lim_{x \to 0} a_{x}^{2} = \frac{k_{x}^{2}}{4}, \lim_{x \to 0} b_{y} = \lim_{x \to 0} b_{y}^{2} = \frac{k_{y}^{2}}{4}$ ;所以当Δx, Δy, Δt趋于0时, 数值弥散关系式 (2.6) 、式 (2.7) 趋向于理想的数值弥散关系: ω2=c2[k $_{x}^{2}$ +k $_{y}^{2}$ ]。此外, 比较式 (2.6) 、式 (2.7) 容易看出SHO-FDTDⅠ和SHO-FDTDⅡ的格式在时间上的精度分别是一阶的和二阶的, 并且SHO-FDTDⅡ比SHO-FDTDⅠ的数值弥散误差要小得多。

2.3 数值弥散误差

下面通过试验求出格式SHO-FDTDⅠ和格式SHO-FDTDⅡ的数值弥散误差, 并与理论结果进行对比分析。令ξ=eiωΔt是方程式 (2.1) 、式 (2.4) 的根。ω=ωR+iωI, 则ξ=e-ωIΔt[cos (ωRΔt) +isin (ωRΔt) ]从而, $\tan (ω_{R} Δ t) = \frac{Ι_{m} (ξ)}{R_{e} (ξ)}$ ; 其中Im (ξ) , Re (ξ) 分别表示ξ的虚部与实部。波的数值相速vp与光速c的比的数值弥散误差: $\frac{v_{p}}{c} = \frac{ω_{R}}{k c} = \frac{1}{c k Δ t} \arctan \frac{Ι_{m} (ξ)}{R_{e} (ξ)} = \frac{Ν_{λ}}{2 π S} \arctan \frac{Ι_{m} (ξ)}{R_{e} (ξ)}$ , 其中S是CFL数, Nλ是一个波长内的节点数, k是波长值。

图4—图6给出了Vp/c在不同情况下的变化曲线。图4给出了SHO-FDTDⅠ与SHO-FDTDⅡ的Vp/c在S=0.35, Nλ=10的条件下随着角度φ的变化图, 从图4中可以看出SHO-FDTDⅡ格式的Vp/c更接近于1。图5给出了两种格式在S=2.4, φ=120°的条件下Vp/c随着Nλ的变化图, 从图5中可以看出SHO-FDTDⅡ格式的数值弥散误差小于SHO-FDTDⅠ的误差。图6是两种格式在φ=65°, Nλ=10的条件下Vp/c随着S的变化图, 从图6中可以看出随着S的增大, 数值弥散误差变得越来越大, 但是在所有的情况下, SHO-FDTDⅡ格式的数值弥散误差比SHO-FDTDⅠ格式的误差要小得多。

参考文献

[1] Yee K S.Numerical solution of initial boundary value problems invol-ving Maxwell’s equations in isotropic media.IEEE Trans Antennasand Propagation, 1966;14:302—307

[2] Zheng F, Chen Z, Zhang J.Toward the development of a three-dimen-sional unconditionally stable finite—difference time—domain method.IEEE Trans Microwave Theory Tech, 2000;48:1550—1558

[3] Namiki T.AnewFDTD algorithmbased on alternatingdirection implicitmethod.IEEE Trans Microwave Theory Tech, 1999;47:2003—2007

[4] Turkel E, Yefet A.Fourth order compact method for the Maxwell e-quations with discontinuous coefficients.Applied Numerical Mathe-matics, 2000;33:125—134

[5] Shang J S.Higher-order compact difference schemes for time depend-ent Maxwell’s equations.Journal of Computational Physics, 2003;153:312—333

[6] Xu L J, Yuan N C.高阶ADI-FDTD算法的数值色散分析.电子信息学报, 2005;27 (3) :1662—1665

[7] Xie Z Q, Zhang B, Chan C H.An explicit fourth-order staggered finitedifference time-domain method for Maxwell’s equations.Journal ofComputational and Applied Mathematics, 2002;147:75—98

[8] Manry C W, Broschat S L, Scneider J B.Higher-order FDTD meth-ods for large problems.J Applied Computational Electromagnetics So-ciety, 1995;10:17—29

时域方法篇9

目前,对于汽油机爆震燃烧产生的机理、爆震的检测、控制、预测等,已有较深入的认识[1];通过测量燃烧室噪声、测量机体震动,或者通过光学观测等手段进行汽油机的爆震检测,已有相关研究探索[2,3,4]。相对而言,基于缸内燃烧压力信号进行爆震的识别,因具有直观、实用、成熟、准确和可靠等特点,在工程中仍然作为最常用方法应用于汽油机爆震识别及爆震强度的评价[5,6]。

通常爆震识别和爆震强度的定量评价主要是基于角基的缸压信号,它要求同时采集缸压信号和曲轴转角信号。曲轴转角信号的采集要求直接安装到曲轴自由端的高精度光电解码器,经要求电源供给的脉冲倍增器输出转角信号和用于TDC修正的触发信号,在工程应用中存在传感器安装困难、测试准备时间长等问题。特别是对于结构紧凑的小型高速汽油机而言,由于没有直接的曲轴自由端,光电解码器很难找到合适和安全的安装位置。基于角基缸压信号的爆震识别及其强度评价也不便于车载状态下的试验检测。

针对上述问题,本文对基于时域缸压信号的爆震检测进行了研究,形成了时域条件下爆震识别和爆震强度评价的动态方法;提出了动态窗口域方法用于准确识别爆震发生和合理评价不同工况条件下的爆震强度;根据时域测量信号的连续性特征,建立了时域缸压信号爆震检测分析系统;对时基和角基爆震识别和强度评价的一致性进行了试验验证。

1 时域缸压信号爆震识别的原理

汽油机爆震直接表现为缸内压力曲线的高频振荡,这种高频压力振荡随时间增长而衰减。对缸压信号分别进行高通滤波和低通滤波,低通滤波信号可以认为是除去高频爆震信号后的压力信号,高通滤波信号可以认为是滤掉低频压力信号得到的高频爆震信号(包含高频噪声信号),爆震压力信号频率范围通常在10～15 kHz[1]。爆震识别的基本原则是对爆震能量的评估,而爆震能量通过爆震振动的累积效应进行评估。这种原则在时域条件下仍然存立。

时域缸压信号的爆震识别将西门子VDO角域爆震识别的基本思想进行扩展,把低通滤波压力信号最大压力峰值对应的时刻t0设置为爆震始点,将爆震始点t0前Δt时间范围设置为参考窗口,爆震始点t0后Δt时间范围设置为爆震窗口,如图1a所示。在参考窗口和爆震窗口内分别对整流后的高频爆震压力信号进行时基积分,将爆震窗口内的压力积分值除以参考窗口内的压力积分值,即得到爆震因子:

undefined

式中,KF为时基爆震因子;pk、pf分别为爆震窗口域和参考窗口域的高通压力。KF可从整流后的高通滤波压力的积分曲线(图1b)上得到:

undefined

式中,E1为至参考窗口始点时爆震能量的累积积分;E0为至爆震始点时爆震能量的累积积分;E2为至爆震窗口终点时爆震能量的累积积分。若KF超过参考值(经验值为2),则认为是一个爆震事件。由于爆震的出现具有随机性,可用一段时间内出现爆震的循环次数占该段时间内总循环次数的百分比(KH)来衡量爆震强度:KH<1 %为无爆震;1 %～5 %为爆震边缘(轻微爆震);5 %～10 %为轻爆震;>10 %为中度爆震以上。

上述算法有两个优点:爆震预测的起始点是变化的,由低通滤波的最高压力决定;发动机不同工况下背景噪声的不同,引入了背景噪声参考窗口,可消除背景噪声对爆震识别的影响。

2 时域爆震窗口域的动态识别

基于时域缸压信号的爆震识别,需要确定参考窗口和爆震窗口的时间宽度Δt。相对不同的发动机转速工况,则Δt相差很大。不过,缸内燃烧过程与曲轴转角之间却具有直接的相关性,故可将时间坐标近似转化到曲轴转角坐标。

通过程序识别采集得到了一组连续的时域缸压低通滤波信号的压力波峰数目,压力波峰数目准确表征发动机的循环数,四冲程发动机可以通过式(3)近似计算发动机转速:

undefined

式中,n为发动机转速,r/min;f为采样频率,Hz;z为采样点数;m为波峰数目。当m足够大时,通过式(3)可以获得精确的发动机转速。在此基础上,可通过式(4)将时间坐标近似与曲轴转角坐标互换:

基于角域的爆震识别方法采用固定窗口宽度,通常取值Δθ=30 °CA。实际上,不同的窗口宽度对爆震识别及其强度评价却具有重要影响,如图2所示。由全负荷、转速分别为4 000、5 000、6 000、7 000 r/min四种工况下,连续采样时域缸压信号进行爆震识别和爆震强度KH的计算可见:采用不同的窗口宽度得到的爆震强度KH存在很大的差异;当窗口宽度<30 °CA时,4 000 r/min转速工况下的爆震强度要大于7 000 r/min转速工况;而窗口宽度>30 °CA时,得到了相反结论。

为了合理评价不同工况条件下的爆震强度,引入爆震能量相对标识参数ΔEI。ΔEI表征高通压力信号在爆震窗口和参考窗口内的累积积分能量之差,其定义为:

理论上可知,当爆震窗口的关闭时刻与爆震过程的实际结束时刻相同时,爆震事件的ΔEI达到最大,如果再延长爆震窗口的宽度,爆震事件的ΔEI理论上应当基本维持不变。

图3a为随机选取的6 000 r/min转速工况中A、B、C、D爆震事件的窗口能量差ΔEI和爆震因子KF随着窗口宽度变化的曲线。其中,1、2、3、4四条曲线分别为A、B、C、D爆震事件的ΔEI变化;而1′、2′、3′、4′四条曲线分别为A、B、C、D爆震事件的KF变化。由图3a可见:随着窗口宽度的增大,则爆震事件的ΔEI逐渐增大;但当窗口宽度增大到一定程度时,ΔEI达到最大并基本保持恒定。对于不同的爆震事件,由于循环的变动,ΔEI达到最大时所对应的窗口宽度是不相同的。图3b为随机选取的6 000 r/min转速工况中,A、B、C、D非爆震事件的ΔEI和KF随窗口宽度的变化。由图3b可见:非爆震事件的爆震窗口内无爆震能量,因此非爆震事件的ΔEI的变化范围很小;C、D非爆震事件的ΔEI出现负值,表明爆震窗口的累积能量小于参考窗口的累积能量。随窗口宽度的改变非爆震事件的ΔEI的变化体现噪声信号的随机特点。

对于明显的非爆震事件,可以采用固定的窗口宽度来计算KF值;而对于可能的爆震事件,爆震窗口的宽度应当根据ΔEI的变化进行动态识别。具体方法为:首先采用不用窗口宽度计算爆震事件的ΔEI,然后取其最大值的96 %～98 %对应的窗口宽度作为爆震窗口域并计算KF。时基爆震窗口域的动态识别能够更准确地进行爆震识别和爆震强度的评估。

3 时域缸压信号的爆震检测分析系统

基于时域缸压信号的爆震检测分析系统主要由原始信号处理、窗口域确定、爆震参数计算和信号处理显示等部分组成,如图4所示。系统连续时域信号的分析处理,如图5所示。

1.测量信号输入 2.低通滤波 3.高通滤波 4.去除噪声 5.爆震始点识别 6.采样点数计算 7.爆震窗口信号提取 8.参考窗口信号提取 9.△EI计算 10.窗口域识别 11.KF计算 12.爆震循环数计算 13.KH计算

原始信号处理即通过带通滤波器从原始缸压信号中将爆震产生的高频振荡压力信号分离出来,随之进行整流和去噪处理,同时将动态数据转化为数组。原始缸压信号的采样频率应当足够大以便能准确地采集爆震信号,通常设置为100 kHz以上。由于原始缸压信号高通滤波后的爆震信号中仍然混有高频干扰信号,考虑到高频干扰信号的振幅相对较小,可通过滤去幅值较小的信号进行去噪处理。带通滤波器频带范围的设置应视爆震产生的高频振荡压力信号的频率范围设定,不同情况下爆震信号频带会有所不同,常用的频带范围取值10～15 kHz,高通滤波和去噪处理由程序中的子模块3和子模块4实现。为了确定爆震始点,首先对原始信号低通滤波并把动态数据转化为数组,然后通过波峰检测器确定波峰的数目m(循环次数)和爆震始点t0。其中t0是通过最高压力数据点在数组中的位置表示的;t0的确定在程序中是通过子模块2和子模块5实现的。

爆震因子计算程序用到了4个for循环、1个while循环和1个条件结构。子模块7和子模块8的for循环为最里层的for循环,分别用于获取爆震窗口和参考窗口的高通压力信号数据,即首先通过公式(3)、(4)得到Δt,并求得Δt时间内包含的数据点数N(N=Δt×f),然后通过for循环内数组索引模块从爆震始点t0向前和向后各取z个数据,得到两个窗口对应的高频振荡压力信号数据。第二层for循环计算不同窗口宽度的爆震能量相对标识参数ΔEI,即通过子模块9对最里层两个for循环确定的高通压力信号进行数值积分并求两者之差,其中不同的窗口宽度是通过改变最里层两个for循环的循环次数实现的。while循环和条件结构确定合适的窗口宽度,即先在转换曲轴转角20～60 °CA区间内采用不同窗口宽度计算求出ΔEI的最大值,然后参照其最大值的96 %～98 %所对应的窗口宽度确定爆震窗口域,并计算该发动机循环的爆震因子KF;对于明显的非爆震事件以转换曲轴转角30 °CA窗口宽度对应的KF作为该发动机循环的爆震因子。最外层的for循环用于发动机连续循环爆震因子的计算,对其中的每个发动机循环首先求取不同窗口宽度对应的KF,然后通过索引得到已确定的合适窗口域对应的KF,即为该发动机循环的爆震因子。最里层的两个for循环的循环次数为Δt0时间段内的数据点数N,其他for循环的循环次数通过索引自动确定。

爆震因子KF的显示是通过for循环索引功能把每循环的KF存放在数组里,然后完成发动机循环数-KF曲线,如图5b所示。KF的大小表征该循环的爆震程度,而爆震强度KH则通过发生爆震的循环数占测试总循环的百分比来表征,KH计算通过子模块12和子模块13实现。

4 时基动态窗口域方法的试验验证

采用角基VDO算法和时基动态窗口域方法进行了爆震识别及其强度评价的对比测试。测试过程中采用Kistler 6052B1微型高温压电晶体压力传感器测量缸内瞬态压力,其自振频率为130 kHz,灵敏度为-200 pc/MPa,热冲击误差小于0.05 MPa,工作温度至400 ℃;采用Kistler 4045A5压阻式压力传感器监测进气瞬时压力,采用Kistler 4075A5压阻式压力传感器监测排气瞬时压力,并设置水冷系统对压阻式传感器进行冷却;上止点位置则采用Kistler 2629电容式上止点位置传感器进行标定。缸内瞬时压力信号和进排气瞬时压力信号分别经过Kistler 5011B和Kistler 4618A0电荷放大器输入DEWETRON-2010数据采集系统。

角基爆震测试采用Kistler 2613B曲轴转角信号发生器,曲轴转角信号与缸压信号同步采集,采样频率取值0.2 °CA,每个工况采集120个循环数据。时基爆震测试则直接对缸内压力连续采样,采样频率为200 kHz,随机提取120个循环进行数据分析。

图6a、图6b、图6c分别示例了表征各种爆震倾向的154FMI、156FMI和171FMI小型高速汽油机的爆震强度曲线。由图6可见:时基和角基的爆震识别和强度评价具有相当的一致性。这表明了时基动态窗口域方法应用于无曲轴自由端的小型高速汽油机的爆震检测,已取得了良好的工程应用效果。

5 结论

(1) 爆震识别的基本原则是对爆震能量的评估,这种原则在时域条件下仍然成立,并可将角域爆震检测的动态方法扩展应用于时域条件下爆震检测,但采用不同的爆震窗口域对爆震强度的评价会产生明显不同的结论。

(2) 爆震事件的爆震能量相对标识参数ΔEI随窗口宽度的增加而逐渐增大,但当窗口宽度增加到一定程度时ΔEI会达到最大并基本保持恒定;据此现象进行爆震窗口域的动态判定可准确识别爆震发生和合理评价不同工况下的爆震强度。

(3) 采用西门子VDO角基算法和时基动态窗口域方法对表征各种爆震倾向的小型高速汽油机进行了爆震检测。测试结果表明两种方法具有良好的一致性。

摘要：针对时基缸压信号的爆震检测,引入爆震能量相对标识参数ΔEI表征高通爆震信号在爆震窗口和参考窗口内的累积积分能量之差。随机选取若干爆震事件和非爆震事件,计算各事件的ΔEI和爆震因子KF随窗口宽度的变化。计算结果表明:爆震事件的ΔEI随窗口宽度的增加而逐渐增大,但当窗口宽度增加到一定程度时,ΔEI会达到最大并基本保持恒定;据此现象进行爆震窗口域的动态判定,可准确识别爆震发生和合理评价不同工况下的爆震强度。分别采用西门子VDO角基算法和时基动态窗口域方法对表征各种爆震倾向的154FMI、156FMI和171FMI小型高速汽油机进行了爆震检测。测试结果表明:基于两种方法的爆震识别和强度评价具有良好的一致性。

关键词：内燃机,爆震识别,爆震强度,信号处理,燃烧测试

参考文献

[1] Heywood J B. Internal combustion engine fundamentals [M]. New York: McGraw-Hill Inc., 1988.

[2] Ollivier E, Bellettre J, Tazerout M, et al. Detection of knock occurrence in a gas SI engine from a heat transfer analysis[J]. Energy Conversion and Management,2006,47(7-8):879-893.

[3] Antoni J, Daniere J, Guillet F. Effective vibration analysis of IC engines using cyclostationarity. part I. a methodology for condition monitoring [J]. Journal of Sound and Vibration,2002,257(5):815-837.

[4]Bianca M,Vaglieco B M.Knock investigation by flame andradical species detection in spark ignition engine for different fu-els[J].Energy Conversion and Management,2007,48(9):2897-2910.

[5]Rogers D.Introducing engine combustion pressure sensor tech-nology[J].Engineering Technology,2006,9(2):45-49.

时域方法篇10

对分布式发电微网系统暂态问题的研究需要根据不同的研究目的对系统进行建模,在此基础上,获取元件的模型参数,采用准确、快速的数值方法仿真计算得到所需结果。在此过程中,元件模型参数的获取通常较为困难,需要和现场物理试验的结果相互校正;对于元件建模,由于暂态过程的时间尺度特性,元件在不同的频率范围内需采用不同的模型加以描述[1,2],如在低频范围内(0.1~3000 Hz)线路模型可用集中参数表示,而在高频范围内则需要采用分布参数线路模型,有时还需要计及元件的频率相关特性;对于仿真计算,则更多地体现在处理各种复杂问题的能力和数值求解的精度与数值稳定性方面,强调程序的建模与仿真能力。

分布式发电微网系统暂态仿真侧重于分布式发电单元、混合发电单元以及微网中各种快速变化的暂态过程的详细仿真,同时兼具含分布式电源及(或)微网的大电网的仿真能力,特别强调仿真结果的准确性和完整性,因此在系统层面采用详细的元件模型对包括电网、电力电子装置、分布式电源及各种控制器进行建模,结合电力系统电磁暂态仿真与电路仿真的基本理论与方法,可以捕捉频率范围从几百k Hz到工频之间系统中的电气量和非电气量的动态过程。

本文基于一般的暂态分析建模需求,从模型实现角度介绍了分布式发电单元中各部分元件的暂态仿真建模方法,根据模型的建模需求描述了暂态仿真程序的框架设计思路,并详细介绍了暂态仿真程序的对象与接口设计。

1 分布式发电系统暂态仿真建模

从结构上看,一个完整的分布式发电单元主要包括分布式电源、电力电子变流器、滤波器、控制器以及电网和负荷,如图1所示。从程序实现上看,各类元件需要在电气及控制系统中分别建模求解。

1.1 分布式电源建模

分布式电源是分布式发电单元的重要组成部分,对不同类型的分布式电源,其动态模型差异较大:一些分布式电源模型以静态特性为主,如光伏电池[3,4,5]、燃料电池[6];另一些则以动态特性为主,如风机[7]、微型燃气轮机[8,9]。其共性在于具有较强的非线性特性,其中大部分属于连续的非线性特性。针对于各种分布式电源模型的实现,可以通过内置于电气系统的完整模型提供给用户,也可以通过控制系统的基本元件由用户自定义完成。对于完整的内置模型建模,一方面由于目前的各种分布式电源模型尚未成熟,很多种类的分布式电源模型还在不断发展之中,内置模型不利于模型的更新与扩展,降低了程序的灵活性;另一方面,在电气系统中各种非线性元件的求解主要基于预测校正法和补偿法,前者的数值稳定性较差,而后者则必须用分布参数线路分隔多个非线性元件[10]。通过上述分析,本文采用控制系统的基本元件对各种分布式电源进行自定义建模。很多分布式电源的数学模型是以传递函数形式表示的,因此自定义的建模方式更易于模型的实现,如图2所示的单轴形式的微型燃气轮机模型。此外,由于控制系统采用基于牛顿迭代的联立求解方法,比预测校正法的数值稳定性更好。

利用控制系统元件对分布式电源建模时,电气系统中测量元件(电压表、电流表等)得到的电气量作为控制系统的输入,而将控制系统的输出以受控源(受控电压源、受控电流源等)的形式模拟分布式电源在电气系统中的外特性。图3给出了包含永磁同步电机模型的微型燃气轮机发电单元原动部分的建模示意图,其中PLL为锁相环[11],PMSM为基于传递函数的永磁同步电机模型,MT为图2给出的微型燃气轮机模型,[T(θ)]和[T(θ)]-1分别为派克(Park)变换和反变换。其他类型的分布式电源(光伏电池、燃料电池、蓄电池等)的建模方法与此类似,这里不再赘述。

在分布式发电单元中,相对于机械传动过程、热力学动态以及化学和电化学反应而言,单元中电场与磁场的相互作用通常具有相对较小的时间常数,根据经验可认为单元中机械变量、热力学变量等反映慢动态过程的参数来不及发生变化,因此可针对不同的研究目的对分布式电源模型进行适当化简,提高仿真程序的计算速度。如在燃料电池发电单元中,由于内部气体分压力的时间常数较大[12],在暂态仿真建模时可近似认为它们为常值。

1.2 电力电子装置建模

除少数的分布式电源可通过同步或异步电机直接并网外,大多数类型的分布式电源需要通过电力电子装置与工频交流电网(或负荷)接口[13],用以解决不同电压等级、不同频率以及交直流系统间的能量传递问题。当前,各种形式的整流器、逆变器及斩波电路等在分布式发电单元中均获得了广泛的应用。

电力电子装置通常由电力电子开关组成。暂态仿真建模时,一般需要计及电力电子装置详细的动态过程以便进行谐波分析等系统级分析。详细模型是利用基本的元件通过拓扑连接实现对电力电子电路的建模,部分商业软件则通过内部封装向用户提供典型的电力电子装置模型。在面向系统级仿真时,包括二极管、晶闸管、IGBT等在内的各种电力电子元件可采用开关模型或双电阻模型表示。在低频场景下的应用中,如果不需要考虑电力电子装置输出的谐波成分,可基于状态空间平均模型对电力电子装置进行化简[14]。

1.3 控制器建模

分布式发电单元的控制从功能上可分为3个层面,它包括分布式电源的控制、电力电子装置自身的控制以及网络层面的电压与频率调节。从功能实现角度而言,这些控制功能多是通过对电力电子装置的控制实现的。光伏发电单元中的最大功率点跟踪MPPT(Maximum Power Point Tracking)就是一种典型的分布式电源控制;电力电子装置的控制主要是根据不同的控制目的采用不同的控制方法,如比例积分控制、滞环控制等,通过单环或双环的控制结构以达到所需的暂态和稳态控制效果;对于网络层面的控制,则根据电网运行的需要对分布式电源运行方式进行调节,可在并网逆变器处采用恒功率控制、U/f控制或下垂控制等控制策略[15]实现分布式电源和电网的交互。

对分布式发电单元中的各种控制器建模宜采用基本控制元件的组合来实现并在控制系统中计算求解。分布式发电单元的控制器模型中常含有较多不连续的非线性环节,也称硬非线性(hard nonlinearity)环节,如比较环节、选择环节等。此时,描述元件输入/输出关系的特性方程会随着输入(输出)量的变化而不同,计算求解时需要检测上述环节的控制逻辑并在元件特性变化时重新形成雅可比矩阵,这同电气系统中的开关模型在本质上是一致的,在应用伪牛顿法求解控制系统时需要特别注意。图4为基于伪牛顿法在一个步长内控制系统求解的基本流程。

1.4 网络及负荷建模

在电网侧,分布式发电单元需要对各种常见的电气元件进行建模,包括线路、变压器、滤波器及负荷等。分布式电源在运行时以接入中低压配网为主,需要考虑配电系统中参数不对称、负载不平衡甚至是接地方式的影响。本文的暂态建模仿真方法能够完全满足配电系统三相详细建模的需求,同时也具有特殊情况下单相建模的能力,能够充分反映不同网络结构、不同接地方式与不同负载水平下对分布式发电单元动态过程的影响。在一般情况下,由于配电线路的供电范围较短,采用集中参数表示的PI型等效电路或串联阻抗模型是合适的[10,16],其电气参数可由线路的几何参数计算得到[17]。在分布式发电单元中,除了作为馈线处的配电变压器实现电压调整功能外,一些情况下分布式电源的并网变压器也用以抑制注入系统的谐波,因此需要着重考虑变压器接线方式的影响。对系统中负荷的建模与传统电磁暂态仿真的负荷模型是一致的,可以采用串联阻抗和(或)异步电机模拟。对于线性阻抗支路、变压器及电机等常规电气元件的建模方法在文献[10,18]中已有介绍,此处不再赘述。

2 总体架构设计

基于前述暂态仿真算法与分布式发电单元暂态仿真的建模需求,本文采用面向对象思想,以C++语言实现了分布式发电微网系统暂态仿真程序TSDG(Transient Simulator for Distributed Generation and microgrid)。从技术层面看,分布式发电微网系统暂态仿真程序可分为3个部分,如图5所示,包括核心计算资源、基础仿真计算及高级功能需求。

2.1 计算资源层

位于最底层的是程序的核心计算资源,它具有线性稀疏方程组求解、非线性方程组求解、线性插值、特征值计算及各种数值积分等暂态仿真程序所需的多种基本数学问题的计算求解能力。部分简单的功能可根据计算需求编程实现,复杂的功能可借助现有成熟可靠的商业软件包来实现,如本文采用开源免费的sparse库[19]实现对线性稀疏方程组的求解。计算资源层与它上层的仿真计算层进行接口,用以处理仿真计算层的各种计算要求,并将计算结果的数据返回给仿真计算层。

2.2 仿真计算层

仿真计算层主要完成各种基本的仿真计算功能,包括稳态计算、暂态仿真、频率扫描等,其核心是暂态仿真功能,一个完整的暂态时域仿真流程如图6所示,它由4个部分组成。

a.读取数据文件。它包括仿真参数和元件参数两部分,前者应提供仿真步长、仿真时间、系统缺省频率等参数,后者应包括元件的基本信息(元件类型与元件名称)、拓扑连接关系(始末节点号)、元件基本参数、数据输出控制参数4个部分。

b.系统划分与初始化。由于控制系统与电气系统的元件特征不同,应在初始化前加以识别;此外,在采用控制系统元件对分布式电源及其控制器建模时,不同的分布式电源及其控制器之间是自然解耦的,也应通过拓扑连接关系加以识别。而对于初始化过程,应分为2步进行:首先,根据读取的数据文件完成模型参数的初始化、元件等效电导的计算及相关准备工作;其次,由稳态平衡点信息对模型中的电气量进行初始化,如历史量的赋值等。

c.暂态时域仿真。它涉及元件模型与系统解算2个层面,主要是对前文介绍的暂态仿真算法的实现。从系统层面看,它包括电气系统与控制系统的交替求解。电气系统解算侧重于电力电子开关模型的动作逻辑,控制系统解算则侧重于各种非线性特性的精确求解。对于元件的功能设计将在下面的详细设计中加以介绍。

d.结果输出与显示。它可以根据不同的数据格式提供灵活和可扩展的数据输出并显示。

仿真计算层通过对计算资源的调用完成仿真计算功能,并将仿真结果提供给在它之上的高级分析功能层。

2.3 高级分析功能层

高级分析功能层位于整个程序架构的最顶端,它主要根据仿真计算的结果进行各种高级应用分析,如谐波分析、参数优化等。高级分析功能层可以通过接口调用仿真计算层的仿真功能,也可以直接调用计算资源层的各种计算资源,完成不同需求的计算功能,如快速傅里叶变换等。这部分的功能实现不是本文的重点。

总体而言,相对较少的计算资源可以完成丰富的仿真计算功能,从而实现不同目的的各种类型的应用。从图5中可以看到,越往下其需求就越是相同或相似。协作和耦合是从较高层到较低层进行的,要避免从底层到高层的耦合。

3 详细设计

根据前面介绍的暂态仿真功能需求,可以从元件(Element)、算例(Case)及解算器(Solver)3个层面对模型和算法进行抽象,如图7所示。

3.1 对象设计

对于元件建模,文献[20]提出了基于支路的建模设计思想,它适用于电气系统的元件建模,但对控制元件则不适用。本文通过对元件模型的抽象,设计了元件基类,考虑到电气系统与控制系统元件的不同特征,分别派生出电气元件基类与控制元件基类,在此基础上对各种电气元件与控制元件进行建模。图8给出了TSDG中对象之间的相互关系。

从图8中可以看出,元件通过继承可以具有电气元件或控制元件的基本属性,在此基础上实现不同元件的行为特征。它们通过聚合关系构成了电气系统类与控制系统类,对不同属性的元件分别进行管理;此外,电气系统类(CElectrical System)、控制系统类(CControl System)及仿真信息类(CSim Config)通过组合实现算例类(CCase)。模型与仿真算法类(CSolver)通过算例类进行接口。

3.2 接口设计

对象通过接口与外部进行交互,接口提供了类的外部视图,隐藏了它的结构和内部行为。暂态仿真程序设计时,电气元件和控制元件在系统层面上各自具有很多相同的行为特征,可以通过接口设计实现对象的多态性,从而简化程序设计。以电气系统为例,通过在电气系统类中对所有电气元件的遍历可以实现系统层面的各种功能调用,图9给出了电气元件基类、电气系统类与电气算法类的一部分接口设计。在电气元件基类的声明中可将部分成员函数置为纯虚函数,要求派生的电气元件必须给出这些成员函数的实现,主要是考虑到各元件的上述实现必不相同。非纯虚的成员函数则没有这个限制,如有的元件模型(电阻、开关等)不参与注入电流源列向量的形成,因此,可不必给出其成员函数formlaug()的实现。

此外,在一些情况下,由第三方提供的程序(库)也需要通过接口设计封装不同的软件实现,为程序提供清晰、统一的接口,如TSDG中采用的sparse库原本是基于C语言实现的,使用时需要重新进行接口的设计与封装。

4 结语

【时域方法】推荐阅读：