基于平均集合卡尔曼滤波的叶栅风洞实验工况参数修正方法

1.本发明属于航空航天实验测量领域，涉及一种基于流场数值模拟和平均集合卡尔曼滤波的叶栅风洞实验工况参数修正方法。


背景技术：

2.叶栅风洞实验是测量叶轮机械叶型性能的常用方法，其工况参数主要指来流攻角和来流马赫数：一般通过转动叶栅安装转盘到对应角度以获得来流攻角；来流马赫数则通过调节风洞气源的阀门或者风机转速，使实验件前部的马赫数传感器显示到给定的马赫数。然而在实际的风洞实验过程中，风洞中的气流受到风洞壁面影响，气流在流经马赫数传感器后存在损失和折转，真实吹到叶栅区域的气流马赫数和攻角并非给定的工况参数。在大攻角、高马赫数时，工况偏差一般更会为明显，严重影响了叶栅风洞实验的效果。
3.现有的叶栅风洞实验工况参数修正方法使用流场数值模拟并利用了叶片表面全部的测量数据，但是在叶片吸力面后50％区域一般会发生流动分离，但通常的湍流模型对分离区预测准确性较差，校正结果会被分离区预测偏差污染，影响工况参数校正结果的准确性。
4.所以目前叶栅风洞实验领域缺乏一种准确、可靠的工况参数修正方法。


技术实现要素：

5.本发明所要解决的技术问题是：
6.为了避免现有修正方法的样本数据存在流动分离导致修正错误的问题，本发明提供一种基于流场数值模拟和平均集合卡尔曼滤波的叶栅风洞实验工况参数修正方法。
7.为了解决上述技术问题，本发明采用的技术方案为：
8.一种基于流场数值模拟和平均集合卡尔曼滤波的叶栅风洞实验工况参数修正方法，其特征在于，包括：
9.步骤1：浮动实验工况参数，构建多组来流工况参数样本，对每组工况参数进行流场数值模拟；
10.步骤2：将浮动工况参数及数值流场的叶片弦长前50％区域的测点位置的相应物理量组成状态向量；
11.步骤3：计算样本协方差矩阵，生成子样本噪声向量及噪声协方差矩阵；
12.步骤4：计算子样本的卡尔曼滤波矩阵；
13.步骤5：对每个子样本进行卡尔曼滤波；
14.步骤6：对所有样本的子样本的状态向量进行平均，从平均后的状态向量中提取修正后的工况参数。
15.本发明进一步的技术方案：步骤1具体如下：根据来流参数的调控精度设置浮动范围，采用拉丁超立方方法进行随机抽样，构建来流参数的样本，在湍流模型对每个样本进行叶栅流场数值计算。
16.本发明进一步的技术方案：其特征在于，步骤2具体如下：
17.将浮动的来流参数及其数值流场结果组建样本状态向量x
f,i
：
[0018][0019]
式中，ma
in,i
和αi分别为来流参数浮动后第i个样本的来流马赫数和来流攻角，表示第i个样本的流场中位于实验中前50％叶片弦长区域第j个测点位置的对应被测物理量，j＝1,2,...,m，m表示验中前50％叶片弦长区域测点的数目。
[0020]
本发明进一步的技术方案：所述的被测物理量包括表面等熵马赫数、静压系数、静压。
[0021]
本发明进一步的技术方案：步骤3具体如下：
[0022]
步骤3-1：计算样本状态向量平均值
[0023][0024]
步骤3-2：计算样本协方差矩阵p：
[0025][0026][0027]
步骤3-3：生成噪声向量w
i,l
和噪声矩阵w
l
：其中l＝1,
…
,l，表示第l个子样本，l为子样本的个数，噪声向量的长度于测点数目相同，即为一组零均值的高斯分布的随机数，其分布方差可参照测量误差，将噪声向量组合成噪声矩阵：
[0028]wl
＝[w
1,l
,w
2,l
,
…
,w
n,l
]
[0029]
步骤3-4：计算子样本的噪声协方差矩阵：
[0030][0031]
本发明进一步的技术方案：步骤4所述的卡尔曼滤波矩阵计算式如下：
[0032]kl
＝ph
t
(hph
t
+r
l
)-1
[0033]
式中，h为观测矩阵：
[0034]
h=[0m×2,im×m]
[0035]
用于将状态向量映射为相应的实验测量值，其中0m×2表示m
×
2维的零矩阵，im×m表示m
×
m维的单位矩阵。
[0036]
本发明进一步的技术方案：步骤5具体如下：
[0037]
对集合样本的子样本进行卡尔曼滤波，得到每个子样本校正后结果x
a,i,l
：
[0038]
x
a,i,l
＝x
f,i
+k
l
(y
exp-hx
f,i
+w
i,l
)
[0039]
式中，y
exp
为叶栅风洞实验中前50％叶片弦长区域测得的表面等熵马赫数组成得列向量。
[0040]
本发明进一步的技术方案：步骤6具体如下：
[0041]
对所有滤波后样本进行平均，获得平均后的状态向量
[0042][0043]
向量前两个数据即为修正后来流马赫数和来流攻角。
[0044]
一种计算机系统，其特征在于包括：一个或多个处理器，计算机可读存储介质，用于存储一个或多个程序，其中，当所述一个或多个程序被所述一个或多个处理器执行时，使得所述一个或多个处理器实现上述的方法。
[0045]
一种计算机可读存储介质，其特征在于存储有计算机可执行指令，所述指令在被执行时用于实现上述的方法。
[0046]
本发明的有益效果在于：
[0047]
本发明提供的一种基于流场数值模拟和平均集合卡尔曼滤波的叶栅风洞实验工况参数修正方法，利用流场数值模拟和平均集合卡尔曼滤波方法修正了叶栅风洞的来流马赫数和攻角，修正结果精度高，即便对偏差较大的参数也能够较好的修正，计算量小；本发明提出的平均集合卡尔曼滤波方法相比于已有的集合卡尔曼滤波方法，更简单，无需迭代，鲁棒性更好。具体如下：
[0048]
1、因为在叶片后50％区域一般会发生流动分离，目前的湍流模型预测流动分离都不十分准确，故后50％的数据被湍流模型预测偏差污染，会带来成修正错误。因此，本发明采用前50％弦长区域的数据而没有采用全部数据，而前50％区域是没有流动分离的附着流动，湍流模型一般能够非常准确地预测。
[0049]
2、本发明采用平均集合卡尔曼滤波方法进行数据同化，提出了子样本的概念，即生成合成噪声的时候，对每个样本生成若干个噪声向量，对每个样本构成了若干个子样本，对每个子样本都进行滤波运算，再进行总体平均。相对于现有技术中的集合卡尔曼滤波方法鲁棒性更好。
[0050]
3、采用平均集合卡尔曼滤波方法，计算结果稳定随机性小，鲁棒性好，结果准确，其需要的样本少，计算量少。在集合卡尔曼滤波方法中有两部分随机性，一个是样本的生成，一个是随机噪声，传统方法仅进行一次平均，无法完全消除随机性，为了结果稳定，只能使用较多的样本，平均集合卡尔曼滤波使用了子样本，进行了两个尺度的平均，更好地消除了结果的随机性，所以稳定性更好，较少的样本就能获得的稳定的结果。
附图说明
[0051]
附图仅用于示出具体实施例的目的，而并不认为是对本发明的限制，在整个附图中，相同的参考符号表示相同的部件。
[0052]
图1为本发明的基于流场数值模拟和平均集合卡尔曼滤波的叶栅风洞实验工况参数修正方法流程图；
[0053]
图2为来流参数修正前后流场计算结果与实验测量结果的对比图。
具体实施方式
[0054]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图和实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。此外，下面描述的本发明各个实施方式中所涉及到的技术特征只要彼此
之间未构成冲突就可以相互组合。
[0055]
如图1所示，本发明提供的一种基于流场数值模拟和平均集合卡尔曼滤波的叶栅风洞实验工况参数修正方法，包含以下步骤：
[0056]
步骤1：浮动工况参数，具体浮动幅值依赖于来流参数的调控精度，一种典型的浮动范围可以设置为：浮动来流马赫数浮动
±
0.05，来流攻角浮动
±1°
，可以采用拉丁超立方等方法进行随机抽样，构建来流参数的样本，对每个样本进行叶栅流场数值计算，本发明对流场数值计算过程中湍流模型不敏感，选取适宜的湍流模型即可；
[0057]
步骤2：将浮动的来流参数及其数值流场结果组建样本状态向量x
f,i
：
[0058][0059]
式中，ma
in,i
和αi分别为来流参数浮动后第i个样本的来流马赫数和来流攻角，表示第i个样本的流场中位于实验中前50％叶片弦长区域第j个测点位置的对应被测物理量，j＝1,2,
…
,m，m表示验中前50％叶片弦长区域测点的数目，被测物理量一般可为表面等熵马赫数、静压系数、静压等；这里采用前50％弦长区域的数据而没有采用全部数据是因为在叶片后50％区域一般会发生流动分离，目前的湍流模型预测流动分离都不是十分准确，故后50％的流场预测数据会被湍流模型预测偏差污染，会带来成修正错误，而前50％区域是没有流动分离的附着流动，湍流模型一般能够非常准确地预测。
[0060]
步骤3：计算样本协方差矩阵p，生成子样本噪声向量w
i,l
及噪声协方差矩阵r
l
，具体步骤如下：
[0061]
步骤3-1：计算样本状态向量平均值
[0062][0063]
步骤3-2：计算样本协方差矩阵p：
[0064][0065][0066]
步骤3-3：生成噪声向量w
i,l
和噪声矩阵w
l
：其中l＝1,
…
,l，表示第l个子样本，l为子样本的个数，噪声向量的长度于测点数目相同，即为一组零均值的高斯分布的随机数，其分布方差可参照测量误差，将噪声向量组合成噪声矩阵：
[0067]wl
＝[w
1,l
,w
2,l
,
…
,w
n,l
]
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0068]
步骤3-4：计算子样本的噪声协方差矩阵：
[0069][0070]
步骤4：计算子样本的卡尔曼增益矩阵k
l
：
[0071]kl
＝ph
t
(hph
t
+r
l
)-1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0072]
式中h为观测矩阵，表示将状态向量映射为对应的实验测量值；
[0073]
h＝[0m×2,im×m]
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0074]
其中0m×2表示m
×
2维的零矩阵，im×m表示m
×
m维的单位矩阵；
[0075]
步骤5：对集合样本的子样本进行卡尔曼滤波，得到每个子样本校正后结果x
a,i,l
：
[0076]
x
a,i,l
＝x
f,i
+k
l
(y
exp-hx
f,i
+w
i,l
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)
[0077]yexp
为实验中前50％叶片弦长区域测得的表面等熵马赫数组成得列向量；
[0078]
步骤6：对所有滤波后样本进行平均，获得平均后的状态向量
[0079][0080]
向量前两个数据即为修正后来流马赫数和来流攻角。
[0081]
实施例1：
[0082]
以man ghh 1-s1压气机叶栅为例，其设计马赫数为0.62，设计攻角为0
°
，实验测量了该叶片压力面和吸力面各10个测点的表面等熵马赫数的分布，对其来流马赫数0.62、来流攻角为4
°
的工况参数进行修正，修正参数之前叶片表面等熵马赫数的数值模拟与实验测量结果如附图2所示，本实例的具体实施步骤如下：
[0083]
步骤1：浮动工况参数，具体浮动幅值依赖于来流参数的调控精度，一种典型的浮动范围可以设置为：浮动来流马赫数ma
in,i
＝0.62
±
0.05，来流攻角浮动αi＝4
°±1°
，可以采用拉丁超立方等方法进行随机抽样，i＝1,2,
…
,n，n为样本的数目，本例中n＝16，对每个样本进行叶栅流场数值计算，本发明对流场数值计算过程中湍流模型不敏感，选取适宜的湍流模型即可，本例中选取sst湍流模型；
[0084]
步骤2：将浮动的来流参数及其数值流场结果组建样本状态向量x
f,i
：
[0085]
x
f,i
＝[ma
in,i
,αi,ma
is,i,1
,ma
is,i,2
,
…
,ma
is,i,m
]
t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0086]
式中，ma
in,i
和αi分别为来流参数浮动后第i个样本的来流马赫数和来流攻角，ma
is,i,j
表示第i个样本的流场中位于实验中前50％叶片弦长区域第j个测点位置的叶片表面等熵马赫数，j＝1,2,
…
,m，m表示验中前50％叶片弦长区域测点的数目，本例中m＝12；
[0087]
步骤3：计算样本协方差矩阵p，生成子样本噪声向量w
i,l
及噪声协方差矩阵r
l
，具体步骤如下：
[0088]
步骤3-1：计算样本状态向量平均值
[0089][0090]
步骤3-2：计算样本协方差矩阵p：
[0091][0092][0093]
步骤3-3：生成噪声向量w
i,l
和这噪声矩阵w
l
：其中l＝1,
…
,l，表示第l个子样本，l为子样本的个数，对于本例l＝300，噪声向量的长度于测点数目相同，即为一组零均值的高斯分布的随机数，其分布方差可参照测量误差，对于本例可取0.01，将噪声向量组合成噪声矩阵：
[0094]wl
＝[w
1,l
,w
2,l
,
…
,w
n,l
]
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0095]
步骤3-4：计算子样本的噪声协方差矩阵：
[0096]
[0097]
步骤4：计算子样本的卡尔曼增益矩阵k
l
：
[0098]kl
＝ph
t
(hph
t
+r
l
)-1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0099]
式中h为观测矩阵，表示将状态向量映射为实验测量值，对于本例：
[0100]
h＝[0m×2,im×m]
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0101]
其中0m×2表示m
×
2维的零矩阵，im×m表示m
×
m维的单位矩阵；
[0102]
步骤5：对集合样本的子样本进行卡尔曼滤波，得到每个子样本校正后结果x
a,i,l
：
[0103]
x
a,i,l
＝x
f,i
+k
l
(y
exp-hx
f,i
+w
i,l
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)
[0104]yexp
为实验中前50％叶片弦长区域测得得表面等熵马赫数组成得列向量；
[0105]
步骤6：对所有滤波后样本进行平均，获得平均后的状态向量
[0106][0107]
向量前两个数据即为修正后来流马赫数和来流攻角，本例中，修正后的工况参数为：来流马赫数0.6349，来流攻角为2.3394
°
。
[0108]
修正后预测结果的等熵马赫数曲线和实验测结果的对比图见附图2，可以看出校正来流参数后，除尾缘区域外，流场数值模拟结果与实验结果很好相符，尾缘区域存在偏差是因为湍流模型对分离区预测能力本身有限，但湍流模型对于附着流动预测准确性较高，本发明在进行参数修正时使用了前50％弦长区域的实验测量数据，该区域为附着流动，从而规避了湍流模型预测偏差的影响。可以看出该工况攻角修正幅值较大，采取神经网络之类的代理模型修正方法无法对训练数据范围外的空间进行外推，本发明在校正工况参数大偏离问题上更具有优势。
[0109]
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明公开的技术范围内，可轻易想到各种等效的修改或替换，这些修改或替换都应涵盖在本发明的保护范围之内。

技术特征：
1.一种基于流场数值模拟和平均集合卡尔曼滤波的叶栅风洞实验工况参数修正方法，其特征在于，包括：步骤1：浮动实验工况参数，构建多组来流工况参数样本，对每组工况参数进行流场数值模拟；步骤2：将浮动工况参数及数值流场的叶片弦长前50％区域的测点位置的相应物理量组成状态向量；步骤3：计算样本协方差矩阵，生成子样本噪声向量及噪声协方差矩阵；步骤4：计算子样本的卡尔曼滤波矩阵；步骤5：对每个子样本进行卡尔曼滤波；步骤6：对所有样本的子样本的状态向量进行平均，从平均后的状态向量中提取修正后的工况参数。2.根据权利要求1所述的一种基于流场数值模拟和平均集合卡尔曼滤波的叶栅风洞实验工况参数修正方法，其特征在于，步骤1具体如下：根据来流参数的调控精度设置浮动范围，采用拉丁超立方方法进行随机抽样，构建来流参数的样本，在湍流模型对每个样本进行叶栅流场数值计算。3.根据权利要求2所述的一种基于流场数值模拟和平均集合卡尔曼滤波的叶栅风洞实验工况参数修正方法，其特征在于，步骤2具体如下：将浮动的来流参数及其数值流场结果组建样本状态向量x
f,i
：式中，ma
in,i
和α
i
分别为来流参数浮动后第i个样本的来流马赫数和来流攻角，表示第i个样本的流场中位于实验中前50％叶片弦长区域第j个测点位置的对应被测物理量，j＝1,2,...,m，m表示验中前50％叶片弦长区域测点的数目。4.根据权利要求3所述的一种基于流场数值模拟和平均集合卡尔曼滤波的叶栅风洞实验工况参数修正方法，其特征在于，所述的被测物理量包括表面等熵马赫数、静压系数、静压。5.根据权利要求4所述的一种基于流场数值模拟和平均集合卡尔曼滤波的叶栅风洞实验工况参数修正方法，其特征在于，步骤3具体如下：步骤3-1：计算样本状态向量平均值1：计算样本状态向量平均值步骤3-2：计算样本协方差矩阵p：2：计算样本协方差矩阵p：步骤3-3：生成噪声向量w
i,l
和噪声矩阵w
l
：其中l＝1,
…
,l，表示第l个子样本，l为子样本的个数，噪声向量的长度于测点数目相同，即为一组零均值的高斯分布的随机数，其分布方差可参照测量误差，将噪声向量组合成噪声矩阵：
w
l
＝[w
1,l
,w
2,l
,...,w
n,l
]步骤3-4：计算子样本的噪声协方差矩阵：6.根据权利要求5所述的一种基于流场数值模拟和平均集合卡尔曼滤波的叶栅风洞实验工况参数修正方法，其特征在于，步骤4所述的卡尔曼滤波矩阵计算式如下：k
l
＝ph
t
(hph
t
＝r
l
)-1
式中，h为观测矩阵：h＝[0
m
×2,i
m
×
m
]用于将状态向量映射为相应的实验测量值，其中0
m
×2表示m
×
2维的零矩阵，i
m
×
m
表示m
×
m维的单位矩阵。7.根据权利要求6所述的一种基于流场数值模拟和平均集合卡尔曼滤波的叶栅风洞实验工况参数修正方法，其特征在于，步骤5具体如下：对集合样本的子样本进行卡尔曼滤波，得到每个子样本校正后结果x
a,i,l
：x
a,i,l
＝x
f,i
+k
l
(y
exp-hx
f,i
+w
i,l
)式中，y
exp
为叶栅风洞实验中前50％叶片弦长区域测得的表面等熵马赫数组成得列向量。8.根据权利要求7所述的一种基于流场数值模拟和平均集合卡尔曼滤波的叶栅风洞实验工况参数修正方法，其特征在于，步骤6具体如下：对所有滤波后样本进行平均，获得平均后的状态向量对所有滤波后样本进行平均，获得平均后的状态向量向量前两个数据即为修正后来流马赫数和来流攻角。9.一种计算机系统，其特征在于包括：一个或多个处理器，计算机可读存储介质，用于存储一个或多个程序，其中，当所述一个或多个程序被所述一个或多个处理器执行时，使得所述一个或多个处理器实现权利要求1所述的方法。10.一种计算机可读存储介质，其特征在于存储有计算机可执行指令，所述指令在被执行时用于实现权利要求1所述的方法。

技术总结
本发明涉及一种平均集合卡尔曼滤波的叶栅风洞实验工况参数修正方法，属于航空航天实验测量域。本发明的基于流场数值模拟和平均集合卡尔曼滤波的叶栅风洞实验工况参数修正方法主要步骤包括：对欲修正实验工况的来流马赫数和攻角进行浮动，构建多组来流参数样本并进行叶栅流场数值计算；使用平均集合卡尔曼滤波方法，修正来实验测量时的流马赫数和来流攻角。本发明能够有效修正叶栅风洞实验的来流边界条件，能有效克服湍流模型产生的分离区预测偏差，修正能力强，无需特意选取湍流模型，计算量小，无需迭代，特别适合于叶栅风洞实验参数的快速修正。的快速修正。的快速修正。


技术研发人员：高丽敏 刘锬韬 但玥 李泽锐
受保护的技术使用者：西北工业大学
技术研发日：2023.04.04
技术公布日：2023/7/12