CN109308341A - 基于加权交叠的可快速计算湍流功率谱密度的方法 - Google Patents

Abstract

本发明涉及一种基于加权交叠的可快速计算湍流功率谱密度的方法，依次包括如下步骤：1)设定热线风速仪采样频率，采集测点的脉动风速，测得有限长度的湍流功率谱密度信号；2)定义测量信号中相邻两点采样的时间间隔T，所需分析信号频带宽度表示为ω至Nmaxω，3)对测量信号进行加矩形窗分段处理，处理后的数据段记为u(t)，每个数据段中采样个数为N，N=1/(nωT)且取整，4)求出最新时刻湍流功率谱密度信号中频率为nω的余弦和正弦的幅值，再求出此时频率为nω的幅值、功率和功率谱密度，5)将各数据段平均后得到频率为nω的功率谱密度；6)依次计算出各频率下的整个数据段的功率谱密度并作图。该方法可减少计算步骤，同时兼顾降低噪声，改善方差特性。

Description

基于加权交叠的可快速计算湍流功率谱密度的方法

技术领域

本发明涉及一种计算湍流功率谱密度的方法，特别涉及一种基于加权交叠的可快速计算湍流功率谱密度的方法，属于实验流体力学技术领域。

背景技术

湍流运动是由众多尺度不同的涡团旋转运动叠加而成，空间上某点的湍流脉动风速是由许多不同尺度的涡在该点处形成的各种频率的脉动叠加而成，因此湍流湍动能可以看作是不同频率涡的贡献。空间上某点的湍流脉动风速由许多不同尺度的涡在该点形成的各种频率的脉动风速叠加而成，因此湍流湍动能可以看作是不同频率涡的贡献。湍流功率谱密度是指湍流脉动动能在频率上的分布密度，用来描述湍流中不同尺度涡的动能对湍流脉动动能的贡献。

目前，湍流功率谱的算法很多，有很多经验谱用来描述大气运动脉动风特性。例如：Davenport谱、Harris谱、Kaimal谱、Simiu谱、Karman谱等等。其中Karman谱是VonKarman在1948年建立的湍流功率谱具有广泛的应用，且Karman谱可较为精确的描述出风洞中均匀湍流特性。其在某处顺风向风谱表达式为：

水平向、垂直向Karman谱表达式为：

式中nω为脉动风频率，S(nω)为湍流功率谱密度，为脉动风速的方差，表示湍流场中特征涡旋的平均尺度，下标h取值为u、v、w分别表示顺风向、水平向与垂直向；那么湍流功率谱密度可以通过脉动风速的时间相关函数经傅里叶变换得到，功率谱估计利用离散傅里叶算法(周期图法)，这样的方法计算量巨大，实时性很难保证，且在某一时刻采样信号的计算与下一时刻采样信号的计算是相互独立，计算得出的功率谱噪声很大。

目前常用的是加窗平均周期法(Welch法)，即对脉动风速数据进行分段，然后对每一段数据加一个窗函数进行预处理(常用的窗函数有矩形窗、汉宁窗等)，然后将每一段数据的自相关函数进行傅里叶变化(FFT)得到谱计算，将每一段数据估计结果的平均值作为整段数据功率谱计算的结果。计算步骤复杂，噪声仍然较高，方差特性差。

发明内容

本发明的目的在于，克服现有技术中存在的问题，提供一种基于加权交叠的可快速计算湍流功率谱密度的方法，可以减少计算步骤，同时兼顾降低噪声，改善方差特性。

为解决以上技术问题，本发明的一种基于加权交叠的可快速计算湍流功率谱密度的方法，依次包括如下步骤：1)设定热线风速仪采样频率，利用热线风速仪采集测点的脉动风速，测得有限长度的湍流功率谱密度信号；2)有限长度的湍流功率谱密度测量信号中，定义相邻两点采样的时间间隔为T，采样频率为1/T，功率谱密度所需分析信号频带宽度表示为ω至N_maxω，且ω为所需分析信号的最低频率，N_max为最高频率与最低频率之间比值；3)对测量信号进行加矩形窗分段处理，处理后的数据段记为u(t)，t为时间变量，其中每个数据段中采样个数为N，每段采样个数N与所需分析的信号频率nω之间存在如下关系：

式中n为1至N_max之间的正整数，int为将一个数值向下取整的函数；4)在第一个数据段中，利用式(2)进行计算，计算完成N步后，则u(t)向后滑移N 个数，就可求得首段湍流功率谱密度所需分析信号中频率为nω的正弦和余弦部分的幅值，利用式(3)将数值转化为信号频率为nω的功率谱密度；

式中，P_n为最新时刻湍流功率谱密度信号中频率为nω的余弦部分的幅值，Q_n为最新时刻湍流功率谱密度信号中频率为nω的正弦部分的幅值，且 nω∈[ω～Nmaxω]；P_n′为之前时刻湍流功率谱密度信号中频率为nω的余弦部分的幅值，Q_n'为之前时刻湍流功率谱密度信号中频率为nω的正弦部分的幅值，第一次迭代计算时，P_n′和Q_n'取值为零；Nc为经处理后的某段功率谱密度信号的最新采样点，u(N_CT)表示为数据段中最新采样点的数值；u((N_c-N)T)表示数据段中采样点为N_C之前第N个点的数值，当N_C处于整个有限长度的测量信号中前N个位置时，u((N_c-N)T)取值为零；F_n为最新时刻湍流功率谱密度信号中频率为nω的幅值，pow_n为最新时刻湍流功率谱密度信号中频率为nω的功率，S_n为最新时刻湍流功率谱密度信号中频率为nω的功率谱密度；5)利用式(2)、式(3)依次计算出下一个数据段中频率为nω的功率谱密度，并将各数据段中频率为nω的功率谱密度进行平均，得到整个数据段中频率为nω的功率谱密度；6)返回步骤3)，依次计算出从ω至N_maxω中各频率下的整个数据段的功率谱密度；7)将得到的整段信号中频率从ω至 N_maxω的功率谱密度数值在对数坐标图上标出，得到测量信号的湍流功率谱密度图。

相对于现有技术，本发明取得了以下有益效果：采样的两个相邻时刻的数据段是有重合的，后面的一个时刻的样本是把前一个时刻的样本的头部舍弃，并将最新的样本添加到最后一个数值得来。本发明在计算时对数据进行加窗处理，分段求功率谱之后再进行平均，无需同时计算出全部N个频谱值，而且对于输入序列的长度可调，具有较强的灵活性，分段加窗的处理及分段重叠的思路，体现了加权交叠的设计思想。本发明对原加窗平均周期法进行了改进，可快速便捷地实现流场中湍流信号的计算，同时降低估算方差的噪声，改善分辨率，分段数据的交叠又减小了每一段的不相关性，使得方差的减小不会达到理论计算的程度；起到了很好改善方差特性的效果。

作为本发明的优选方案，步骤3)中相邻的矩形窗可以部分重叠，利用式 (4)进行计算：

式中，相邻两段采样重复点数为m，u((N_c+i-N)T)表示数据段中采样点为 N_C之前第N-i个点的数值，且i取值为0至m-1之间的整数。窗函数部分重叠能够降低估算方差的噪声，改善分辨率，分段数据的交叠又减小了每一段的不相关性，使得方差减小，更好地改善了方差特性，能够更直观的对流场中湍流功率谱密度进行分析。

作为本发明的进一步优选方案，步骤3)中对测量信号加矩形窗分段处理时，所使用的矩形窗函数的表达式如式(5)：

式中，ω(t)为宽度为N的矩形窗函数，利用ω(t)截断测量信号得到u(t)。数据处理时存在信号相当长的情况，频谱分析时需把信号分成有限长度的数据段，相当于对信号进行了加矩形窗的处理，为了减少频谱能量泄漏，采用选用矩形窗函数对信号进行截断，窗口内的系数乘以1，窗口外的系数乘以0，矩形窗可以使得能量相对集中在主瓣当中，就可以保证信号和真实的频谱较为接近。

作为本发明的进一步优选方案，步骤4)的式(2)和式(4)中功率谱密度计算时所用sin、cos函数应分别满足式(6)、式(7)，

其中j表示为N_C在整个有限长度的湍流功率谱密度测量信号中所代表的位置，取值为从0开始的整数。sin、cos函数可按照上述公式预先处理好，数据处理时可直接代入进行计算，无需实时计算，大幅提高了湍流功率谱密度的生成速度，具有重要的工程意义。

附图说明

下面结合附图和具体实施方式对本发明作进一步详细的说明，附图仅提供参考与说明用，非用以限制本发明。

图1为本发明的算法流程框图。

图2为本发明具体实施方式的一种基于加权交叠的可快速计算湍流功率谱密度的方法与传统周期图法的功率谱密度分布对比图。

图3为本发明实施例中在不同格栅距离下的功率谱密度分布图。

图4为本发明实施例中在不同湍流强度下功率谱密度分布图。

图5为本发明实施例中在不同风速工况下功率谱密度分布图。

具体实施方式

如图1所示，本发明基于加权交叠的可快速计算湍流功率谱密度的方法，依次包括如下步骤：

1)设定热线风速仪采样频率，利用热线风速仪采集测点的脉动风速，测得有限长度的湍流功率谱密度信号；

2)有限长度的湍流功率谱密度测量信号中，定义相邻两点采样的时间间隔为T，采样频率为1/T，功率谱密度所需分析信号频带宽度表示为ω至N_maxω，且ω为所需分析信号的最低频率，N_max为最高频率与最低频率之间比值；

3)对测量信号进行加矩形窗分段处理，处理后的数据段记为u(t)，t为时间变量，其中每个数据段中采样个数为N，每段采样个数N与所需分析的信号频率nω之间存在如下关系：

式中n为1至N_max之间的正整数，int为将一个数值向下取整的函数；

4)在第一个数据段中，利用式(2)进行计算，计算完成N步后，则u(t)向后滑移N个数，就可求得首段湍流功率谱密度所需分析信号中频率为nω的正弦和余弦部分的幅值，利用式(3)将数值转化为信号频率为nω的功率谱密度；

式中，P_n为最新时刻湍流功率谱密度信号中频率为nω的余弦部分的幅值，Q_n为最新时刻湍流功率谱密度信号中频率为nω的正弦部分的幅值，且 nω∈[ω～N_maxω]；P_n′为之前时刻湍流功率谱密度信号中频率为nω的余弦部分的幅值，Q_n'为之前时刻湍流功率谱密度信号中频率为nω的正弦部分的幅值，第一次迭代计算时，P_n′和Q_n'取值为零；Nc为经处理后的某段功率谱密度信号的最新采样点，u(N_CT)表示为数据段中最新采样点的数值； u((N_c-N)T)表示数据段中采样点为N_C之前第N个点的数值，当N_C处于整个有限长度的测量信号中前N个位置时，u((N_c-N)T)取值为零；F_n为最新时刻湍流功率谱密度信号中频率为nω的幅值，pow_n为最新时刻湍流功率谱密度信号中频率为nω的功率，S_n为最新时刻湍流功率谱密度信号中频率为nω的功率谱密度；

5)利用式(2)、式(3)依次计算出下一个数据段中频率为nω的功率谱密度，并将各数据段中频率为nω的功率谱密度进行平均，得到整个数据段中频率为nω的功率谱密度；

6)返回步骤3)，依次计算出从ω至N_maxω中各频率下的整个数据段的功率谱密度；

7)将得到的整段信号中频率从ω至N_maxω的功率谱密度数值在对数坐标图上标出，得到测量信号的湍流功率谱密度图。

为兼顾方差与分辨率，步骤3)中相邻的矩形窗可以部分重叠，利用式(4) 进行计算：

式(4)中，相邻两段采样重复点数为m，u((N_c+i-N)T)表示数据段中采样点为N_C之前第N-i个点的数值，且i取值为0至m-1之间的整数。

步骤3)中所使用的矩形窗函数的表达式如式(5)：

式(5)中，ω(t)为宽度为N的矩形窗函数，利用ω(t)截断测量信号得到 u(t)。

步骤4)的式(2)和式(4)中功率谱密度计算时所用sin、cos函数应分别满足式(6)、式(7)，

其中j表示为N_C在整个有限长度的湍流功率谱密度测量信号中所代表的位置，取值为从0开始的整数。

在测量信号中任何时刻获取的湍流功率谱密度测量信号可表示为各频率分量累加：

u(st)表示为任意一段湍流功率谱密度测量信号，s取值为0至N-1整数，那么P_n和Q_n可如式(9)表示：

式(9)中，u(kT)表示为数据段中位置为k的采样点的数值，k的取值为0 至N-1的整数，那么只要根据采样值计算出P_n和Q_n，湍流功率谱信号中频率为nω的幅值F_n、功率pow_n、功率谱密度S_n即可以表示成式(3)。

实施例一

以湍流功率谱密度所需分析信号的最低频率ω的计算为例，对式(9)进行改进可得：

式(10)中，N_c为经处理后的某段功率谱密度信号的最新采样点，显然可以看出最新时刻的P₁和Q₁和前一时刻的P₁’、Q₁’之间存在递归关系，如式(11) 所示：

P₁为最新时刻湍流功率谱密度信号中频率为ω的余弦部分的幅值，Q₁表示为最新时刻湍流功率谱密度信号中频率为ω的正弦部分的幅值，P₁'表示为之前时刻湍流功率谱密度信号中频率为ω的余弦部分的幅值，Q₁'为之前时刻湍流功率谱密度信号中频率为ω的正弦部分的幅值，u(N_CT)表示为数据段中最新的采样点的数值，u((N_c-N)T)表示数据段中采样点为N_c之前第N个点的数值，当Nc处于整个有限长度的测量信号中前N个位置时，u((N_c-N)T)取值为零；式(11)中功率谱密度计算时所用sin、cos函数应分别满足式(6)。

为兼顾方差与分辨率，计算时各段数据会有部分重叠，重叠个数用m表示，即新的一组采样数据至N_c+N-m结束，m的选择可依据所需的分辨率要求与数据长度综合考虑，可将式(11)改写成如式(12)所示：

式(12)中，相邻两段采样重复点数为m，u((N_c+i-N)T)表示数据段中最新采样点之前第N-i个点的数值，且i取值为0至m-1整数，式(12)中功率谱密度计算时所用sin、cos函数应分别满足式(7)。

因此，在已知之前时刻的正弦部分的幅值和余弦部分幅值的基础上，只需简单的递归计算就可以得到当前时刻的新值。计算时可先对采样个数为N 的首段利用式(11)完成迭代计算作为初始化，当u(t)向后滑移N个数后，即计算完成N步后，可利用式(3)求得首段湍流功率谱密度所需分析信号频率为ω的总幅值，随后可利用式(12)、式(3)对下一段湍流功率谱密度迭代计算，再将各段计算所得到的功率谱密度进行总平均，可得到整段信号中频率为ω的功率谱密度，依照上述方法将所需分析信号频带宽度即频率从ω至N_maxω依次求出，在对数坐标图上标出即可得到整段信号中频率从ω至N_maxω的湍流功率谱密度图。

实施例二

定义相邻两点采样的时间间隔为T＝0.0002s，采样频率1/T即为5kHz，所需分析信号的最低频率ω取值为2Hz，Nmax取值为750，nω∈[2～1500]， n为1至750之间的正整数，对式(10)进行改进可得：

N_c为经处理后的某段功率谱密度信号的最新采样点，显然可以看出前一时刻的P₁’、Q₁’和前最新时刻的P₁和Q₁之间存在递归关系，如式(7)所示：

表示为最新时刻的湍流功率谱密度信号中频率为2Hz余弦部分的幅值，表示为最新时刻的湍流功率谱密度信号中频率为2Hz正弦部分的幅值，表示为之前时刻的湍流功率谱密度信号中频率为2Hz的余弦部分的幅值，为之前时刻湍流功率谱密度信号中频率为2Hz的正弦部分的幅值，第一次迭代计算时，取值为零。u(N_CT)表示为数据段中最新的采样点的数值， u((N_c-2500)T)表示数据段中采样点为最新采样点之前第2500个点的数值，当Nc处于整个有限长度的测量信号中前N个位置时，u((N_c-N)T)取值为零， j表示为Nc在整个有限长度的湍流功率谱密度测量信号中所代表的位置，取值为从0开始的整数。

初始化各变量后，利用式(14)进行迭代计算，当u(t)向后滑移2500个数后，即计算完成2500步后，利用式(3)求得此段湍流功率谱密度所需分析信号频率为2Hz的正弦部分和余弦部分的幅值。实际应用中，为兼顾方差与分辨率，计算时各段数据需部分重叠，假定重复率50％，重叠个数m为1250，即新的一组采样数据至N_c+2500-1250结束，可将式(14)改写成如式(15)所示：

式(15)中，u((N_c+i-2500)T)表示数据段中最新采样点之前第2500-i个点的数值，且i取值为0至m-1整数。按式(15)进行离散傅里叶的计算之前，可先利用式(14)对首段(采样个数为2500)完成求和运算作为初始化，再利用式(3)求得首段湍流功率谱密度计算，随后对下一段湍流功率谱密度计算时，随后可利用式(15)、式(3)对下一段湍流功率谱密度迭代计算，再将各段计算所得到的功率谱密度进行总平均，可得到整段信号中频率为ω的功率谱密度，依照上述方法将频率从2至1500Hz的功率谱密度依次求出，即得到整个信号的湍流功率谱密度。

本发明的计算方法与传统周期图法的功率谱密度分布对比如图2所示，其中图2a、图2b分别为风洞中加入被动格栅制造湍流场且风速为20m/s时，在距离格栅截面中心测点的2m位置处，本发明的算法与传统的利用离散傅里叶算法测得的顺流向湍流功率谱密度。可以看出采用本发明的方法后，可以明显改善功率谱密度的分辨率，顺风向在大于30Hz后功率谱密度急速下降，能量涡团产生的频率都在100Hz以下，且实际测量的功率谱密度与Karman 经验谱吻合很好，验证了本发明方法的正确性。

本发明的计算方法在不同工况下功率谱密度对比如图3、图4、图5所示，其中图3为本发明的计算算法在风洞中风速为20m/s时格栅后1m、2m处测点顺流向湍流功率谱密度分布，图4为本发明的计算算法在风速为20m/s时格栅后2m位置处在10.5％、19％湍流强度工况下顺流向湍流功率谱密度分布，图5为本发明的计算算法在15m/s、25m/s风速下湍流强度为10.5％时顺流向湍流功率谱密度分布。从图3、图4、图5皆可看出采用本发明的算法后，实测谱与Karman谱吻合较好，且相同条件下距离格栅距离越远功率谱密度函数越小；流场湍流强度越大，功率谱密度越大；来流风速越大，功率谱密度越大；符合湍流实际变化规律，说明本发明的计算方法适用于不同工况下风洞均匀湍流场功率谱密度的测量。

下表为本发明与传统快速傅里叶算法(FFT)的计算次数对照表，

	快速傅里叶(FFT)	本发明
			加法次数	0.5*N*log<sub>2</sub>N	N
乘法次数	N*log<sub>2</sub>N	2N

从上表的计算次数可以看出，本发明的计算方法远比传统快速傅里叶算法简便。

以上所述仅为本发明之较佳可行实施例而已，非因此局限本发明的专利保护范围。除上述实施例外，本发明还可以有其他实施方式。凡采用等同替换或等效变换形成的技术方案，均落在本发明要求的保护范围内。本发明未经描述的技术特征可以通过或采用现有技术实现，在此不再赘述。

Claims (4)

1.一种基于加权交叠的可快速计算湍流功率谱密度的方法，其特征在于，依次包括如下步骤：

1)设定热线风速仪采样频率，利用热线风速仪采集测点的脉动风速，测得有限长度的湍流功率谱密度信号；

2)有限长度的湍流功率谱密度测量信号中，定义相邻两点采样的时间间隔为T，采样频率为1/T，功率谱密度所需分析信号频带宽度表示为ω至N_maxω，且ω为所需分析信号的最低频率，N_max为最高频率与最低频率之间比值；

3)对测量信号进行加矩形窗分段处理，处理后的数据段记为u(t)，t为时间变量，其中每个数据段中采样个数为N，每段采样个数N与所需分析的信号频率nω之间存在如下关系：

式中n为1至N_max之间的正整数，int为将一个数值向下取整的函数；

4)在第一个数据段中，利用式(2)进行计算，计算完成N步后，则u(t)向后滑移N个数，就可求得首段湍流功率谱密度所需分析信号中频率为nω的正弦和余弦部分的幅值，利用式(3)将数值转化为信号频率为nω的功率谱密度；

式中，P_n为最新时刻湍流功率谱密度信号中频率为nω的余弦部分的幅值，Q_n为最新时刻湍流功率谱密度信号中频率为nω的正弦部分的幅值，且nω∈[ω～Nmaxω]；P_n'为之前时刻湍流功率谱密度信号中频率为nω的余弦部分的幅值，Q_n'为之前时刻湍流功率谱密度信号中频率为nω的正弦部分的幅值，第一次迭代计算时，P_n'和Q_n'取值为零；Nc为经处理后的某段功率谱密度信号的最新采样点，u(N_CT)表示为数据段中最新采样点的数值；u((N_c-N)T)表示数据段中采样点为N_C之前第N个点的数值，当N_C处于整个有限长度的测量信号中前N个位置时，u((N_c-N)T)取值为零；F_n为最新时刻湍流功率谱密度信号中频率为nω的幅值，pow_n为最新时刻湍流功率谱密度信号中频率为nω的功率，S_n为最新时刻湍流功率谱密度信号中频率为nω的功率谱密度；

5)利用式(2)、式(3)依次计算出下一个数据段中频率为nω的功率谱密度，并将各数据段中频率为nω的功率谱密度进行平均，得到整个数据段中频率为nω的功率谱密度；

6)返回步骤3)，依次计算出从ω至Nmaxω中各频率下的整个数据段的功率谱密度；

7)将得到的整段信号中频率从ω至Nmaxω的功率谱密度数值在对数坐标图上标出，得到测量信号的湍流功率谱密度图。

2.根据权利要求1所述的基于加权交叠的可快速计算湍流功率谱密度的方法，其特征在于，步骤3)中相邻的矩形窗可以部分重叠，利用式(4)进行计算：

式中，相邻两段采样重复点数为m，u((N_c+i-N)T)表示数据段中采样点为N_C之前第N-i个点的数值，且i取值为0至m-1之间的整数。

3.根据权利要求1所述的基于加权交叠的可快速计算湍流功率谱密度的方法，其特征在于，步骤3)中对测量信号加矩形窗分段处理时，所使用的矩形窗函数的表达式如式(5)：

式中，ω(t)为宽度为N的矩形窗函数，利用ω(t)截断测量信号得到u(t)。

4.根据权利要求1所述的基于加权交叠的可快速计算湍流功率谱密度的方法，其特征在于：步骤4)的式(2)和式(4)中功率谱密度计算时所用sin、cos函数应分别满足式(6)、式(7)，

其中j表示为N_C在整个有限长度的湍流功率谱密度测量信号中所代表的位置，取值为从0开始的整数。