0 引言

现代飞行器的发展对高速、隐身性能的要求越来越高。传统的安装在飞行器前端的空速管等探针式大气数据传感器对飞行器的隐身、机动性能的不利影响较大。嵌入式大气数据传感系统(Flush Air Data Sensing System, FADS System)是依靠嵌入在飞行器头部的若干压力传感器来获得飞行器表面压力分布，并由此压力分布间接得出飞行参数的系统，曾用于X-33、X-34和X-38等飞行器上^[1]。由于压力传感器分布与飞行器表面平齐，相比传统的探针式传感器系统，这种嵌入式传感系统在大马赫数、大迎角下依然能够很好地工作，且有利于提高飞行器隐身性能。

早期对FADS系统的研究主要针对球形头部的飞行器，提出了基于圆球绕流理论解的三点法^[1-3]。后来发展了针对该系统的神经网络算法^[4-9]，通过大量的试验数据训练，适用于非线性系统建模。

模糊逻辑建模方法以往多用于对非定常气动力的建模^[10-12]或者工业控制领域。本文研究了模糊逻辑对FADS系统的建模。以模型表面若干测压点的压力或压力系数作为模糊系统的输入，以迎角、侧滑角、来流速度和海拔高度作为输出，混合使用梯度下降法和最小二乘估计法来识别模糊系统的参数，从而建立针对该系统的模糊逻辑模型。该模型建立完成之后，将它运用于飞行器的传感系统，可以实时地根据若干测压点的压力得到飞行器的飞行状态参数，然后传递给控制系统。

1 风洞试验及CFD仿真 1.1 风洞试验部分

试验在南京航空航天大学NH-2低速风洞进行。风洞试验段截面为3.0 m×2.5 m的矩形，试验模型为某飞翼布局飞行器，表面共布置62个测压点，头部分布46个测压点，内翼左右各分布8个测压点，位置基本对称(图 1、图 2)。测压采用内埋测压管的方式，测压管内径0.3 mm。每个测压管均与压力测量模块相连。

图 3中给出了半侧机身上的测压点分布，除1~8号测压孔外，另半侧机身上的测压孔分布与图中“Z”开头编号的测压孔对称，相对应的编号以“Y”开头，例如“Z18”对应“Y18”。

机头为坐标系原点。x轴位于飞机对称面指向向后，与机身水平线平行；y轴向上为正；z轴指向左机翼为正。模型展长3.5 m，部分关键测压点坐标如表 1所示。

试验采用DTC Initium电子扫描系统采集表面压力。该系统有一台主控制器，可带8个测压模块，每个模块可连接64个测压通道，数据采样频率330 Hz，系统精度0.1%。试验采用1个量程为1 psi的测压模块测量模型表面压力。

试验风速为35 m/s、55 m/s和70 m/s，试验迎角-18°~24°，每隔2°进行试验，侧滑角0°~16°，每隔4°进行试验。最终得到在不同迎角、侧滑角和风速情况下的模型表面62个测压点的压力和压力系数值。

由于地面风洞试验不能模拟不同海拔高度的情况，所以辅助使用CFD仿真来计算在不同海拔高度下的流场。

1.2 CFD仿真部分

流场网格使用四面体/棱柱体混合网格，使用Fluent求解器，湍流模型采用S-A模型。来流速度为35 m/s、70 m/s，迎角为-10°、0°、10°，侧滑角为0°、8°，海拔高度为500 m、1000 m、1500 m、2000 m、3000 m、4000 m、5000 m和6000 m共8个高度。一共96个计算状态。CFD计算得到的压力数据用于对海拔高度的建模。

2 模糊逻辑建模 2.1 模糊结构

本文采用Takagi和Sugeno提出过一种非线性模糊模型^[13]，规则是：If x₁ is A and x₂ is B, then y=f(x₁, x₂)。A和B表示模糊集合。当f(x₁, x₂)是一阶多项式时，称为一阶Sugeno模型。将输入量的个数扩展至k个，if x₁ is Aⁱ₁(x₁), x₂ is Aⁱ₂(x₂)…，x_k is A_kⁱ(x_k)，此时输出表示为：

$ P^{i}=p_{0}^{i}+p_{1}^{i} x_{1}+p_{2}^{i} x_{2}+\cdots+p_{k}^{i} x_{k} $

(1)

A_rⁱ(x_r)表示x_r的隶属函数，Pⁱ叫做内部函数。若对于输入变量x_r划分N_r个隶属函数，则内部函数的个数n=N₁×N₂×…×N_r×…×N_k, i=1~n。

下面以2个输入量为例，说明该模糊结构的建立。

图 4以5层为例展示模糊结构。

第1层表示对于两个输入量x₁、x₂各分配若干个隶属函数，分别是A₁、A₂、B₁、B₂、B₃。

第2层表示对分别属于各输入变量的隶属函数进行交叉乘积计算，得到内部函数的权重值，分别是w₁、w₂、w₃、w₄、w₅、w₆。

第3层表示将各权重值进行归一化处理，得到权重系数w_i。

$ {{\bar w}_i} = \frac{{{w_i}}}{{\sum\limits_{i = 1}^6 {{w_i}} }} $

(2)

第4层表示将权重系数与内部函数相乘。

第5层表示将第4层的输出进行求和运算，得到最终输出y。

$ y = \sum\limits_{i = 1}^6 {{{\bar w}_i}} {P^i} $

(3)

2.2 隶属函数

识别隶属函数参数时，根据不同的识别方法，分别使用了三角形^[10]、梯形和钟形^[14]隶属函数。本文使用到的3种不同形状的隶属函数示意图如图 5~图 7所示。

钟形隶属函数表达式：

$ A_{r}\left(x_{r}\right)=\frac{1}{1+\left[\left(\frac{x_{r}-c_{i}}{a_{i}}\right)^{2}\right]^{b_{i}}} $

(4)

此函数由3个参数a_i、b_i、c_i决定，在进行参数识别的时候利用梯度下降法动态调整这3个参数。

2.3 参数识别 2.3.1 隶属函数参数识别

参考图 4的模糊结构，记第k层具有个#(k)个节点，把第k层的第i个节点的位置记做(k, i)。每个节点的输出O_i^k依赖于之前的输入和参数集合，该部分算法见文献[14]，即：

$ O_{i}^{k}=O_{i}^{k}\left(O_{i}^{k-1}, \cdots, O_{\#(k-1)}^{k-1}, a, b, c, \cdots\right) $

(5)

式中a, b, c, …表示该节点包含的参数。假设共有m组样本数据，对于第j组样本数据，定义对某组样本数据的平方差E_j：

$ {E_j} = \sum\limits_{i = 1}^{\# (L)} {{{\left( {T_{i, j}^L - O_{i, j}^L} \right)}^2}} $

(6)

T_{i, j}^L表示第j组样本数据在第L层结构的第i个分量的目标输出(对于单输出系统，目标输出个数是1)，O_{i, j}^L表示第j组样本数据在第L层结构的第i个分量的实际输出(对于单输出系统，实际输出个数是1)，则所有样本的平方差之和$E = \sum\limits_{j = 1}^m {{E_j}} $。

为求解E对隶属函数参数的偏导数，首先求解∂E_j/∂O。对于第L层结构可以根据E_j的表达式直接推导出来。

$ \frac{\partial E_{j}}{\partial O_{i, j}^{L}}=-2\left(T_{i, j}^{L}-O_{i, j}^{L}\right) $

(7)

对于第L层之前的结构，即对第k层结构(1≤k≤L-1)在(k, i)位置处的节点输出O_{i, j}^k，根据链式求导法则得到：

$ \frac{{\partial {E_j}}}{{\partial O_{i, j}^k}} = \sum\limits_{l = 1}^{\# (k + 1)} {\frac{{\partial {E_j}}}{{\partial O_{l, j}^{k + 1}}}} \frac{{\partial O_{l, j}^{k + 1}}}{{\partial O_{i, j}^k}} $

(8)

根据式(7)和式(8)，可以求出对任意节点处的∂E_j/∂O。

现在来求解E_j对隶属函数参数的偏导数。假设α是其中一个隶属函数的某个参数。

$ \frac{{\partial {E_j}}}{{\partial \alpha }} = \sum\limits_{o* \in S} {\frac{{\partial {E_j}}}{{\partial {O^*}}}} \frac{{\partial {O^*}}}{{\partial \alpha }} $

(9)

S表示节点输出依赖于α的节点集合。因为有m组样本数据，则：

$ \frac{{\partial E}}{{\partial \alpha }} = \sum\limits_{j = 1}^m {\frac{{\partial {E_j}}}{{\partial \alpha }}} $

(10)

相应的，参数更新变化量为$ \Delta \alpha=-\eta \frac{\partial E}{\partial \alpha}$。其中η表示更新速率，可以表示为：

$ \eta=\frac{k}{\sqrt{\sum\limits_{\alpha}\left(\frac{\partial E}{\partial \alpha}\right)^{2}}} $

(11)

k表示在参数空间内梯度变化的步长。在编程计算时可以通过改变k来加快收敛。例如当平方差之和E连续四次迭代降低时，可以将步长变为1.25k，当E出现增加减少来回波动的情况时，将步长变为0.8k。

2.3.2 内部函数参数识别

可以使用梯度下降法或者最小二乘估计法来识别内部函数参数。

模糊输出：

$ \begin{array}{*{20}{l}} {{{\hat y}_j} = }\\ {\frac{{\sum\limits_{i = 1}^n w \left[ {A_1^i\left( {{x_{1, j}}} \right), \cdots , A_r^i\left( {{x_{r, j}}} \right), \cdots , A_k^i\left( {{x_{k, j}}} \right)} \right]{P^i}}}{{\sum\limits_{i = 1}^n w \left[ {A_1^i\left( {{x_{1, j}}} \right), \cdots , A_r^i\left( {{x_{r, j}}} \right), \cdots , A_k^i\left( {{x_{k, j}}} \right)} \right]}}} \end{array} $

(12)

其中内部函数Pⁱ=p₀ⁱ+p₁ⁱx₁+p₂ⁱx₂+…+p_kⁱx_k。w[A₁ⁱ(x_{1, j}), …, A_rⁱ(x_{r, j}), …, A_kⁱ(x_{k, j})]表示各隶属函数的交叉乘积运算，个数n=N₁×N₂×…×N_r×…×N_k。

2.3.2.1 梯度下降法

模糊输出与样本输出的平方差之和：

$ E = \sum\limits_{j = 1}^m {{{\left( {{{\hat y}_j} - {y_j}} \right)}^2}} $

(13)

m是样本个数。定义相关系数：

$ {R^2} = 1 - \frac{{\sum\limits_{j = 1}^m {{{\left( {{{\hat y}_j} - {y_j}} \right)}^2}} }}{{\sum\limits_{j = 1}^m {{{\left( {\bar y - {y_j}} \right)}^2}} }} $

(14)

y表示样本输出的平均数。R²越接近于1，建模精度越高。

待定的内部函数参数的个数是N₁×N₂×…×N_r×…×N_k×(k+1)。

目标模糊逻辑结构要使得E最小，使用梯度下降法更新内部函数参数。

$ p_{r, t+1}^{i}=p_{r, t}^{i}-\alpha_{r} \frac{\partial(E)}{\partial p_{r}^{i}} $

(15)

α_r是收敛因子，在(0, 1)之间。当α_r选取的较小时，建模速度慢，当α_r选取过大时，计算出现溢出。需要不断调整选取合适的值，最大化建模速度。

将E的表达式代入展开，得到：

当r=0时，

$ \begin{array}{*{20}{c}} {p_{0, t + 1}^i = p_{0, t}^i - 2{\alpha _0}\sum\limits_{j = 1}^m {\left[ {\left( {{{\hat y}_j} - {y_j}} \right) \cdot } \right.} }\\ {\left. {\frac{{w\left[ {A_1^i\left( {{x_{1, j}}} \right), \cdots , A_r^i\left( {{x_{r, j}}} \right), \cdots , A_k^i\left( {{x_{k, j}}} \right)} \right]}}{{\sum\limits_{l = 1}^n w \left[ {A_1^l\left( {{x_{1, j}}} \right), \cdots , A_r^i\left( {{x_{r, j}}} \right), \cdots , A_k^l\left( {{x_{k, j}}} \right)} \right]}}} \right]} \end{array} $

(16)

当r=1, 2, 3, …, k时，

$ \begin{array}{*{20}{c}} {p_{r, t + 1}^i = p_{r, t}^i - 2{\alpha _r}\sum\limits_{j = 1}^m {\left[ {\left( {{{\hat y}_j} - {y_j}} \right) \cdot } \right.} }\\ {\left. {\frac{{w\left[ {A_1^i\left( {{x_{1, j}}} \right), \cdots , A_r^i\left( {{x_{r, j}}} \right), \cdots , A_k^i\left( {{x_{k, j}}} \right)} \right]{x_{r, j}}}}{{\sum\limits_{l = 1}^n w \left[ {A_1^l\left( {{x_{1, j}}} \right), \cdots , A_r^l\left( {{x_{r, j}}} \right), \cdots , A_k^l\left( {{x_{k, j}}} \right)} \right]}}} \right]} \end{array} $

(17)

2.3.2.2 最小二乘估计法

模糊结构可以表示为AP=B。

P表示内部函数参数构成的列向量矩阵，元素个数为(k+1)×N₁×…×N_r×…×N_k。A矩阵由样本输入和隶属函数计算得到，行数为样本数m，列数为(k+1)×N₁×…×N_r×…×N_k。B矩阵表示样本输出，元素个数为样本数m。对此一阶线性方程组，通常使用最小二乘法求解。

将P矩阵的最小二乘估计记做P^*，则：

$ P^{*}=\left(A^{\mathrm{T}} A\right)^{-1} A^{\mathrm{T}} B $

(18)

A^TA是非奇异矩阵时，它的逆存在。当矩阵维数较大时，直接计算逆矩阵的效率很低。如果样本数量m < (k+1)×N₁×…×N_r×…×N_k，上述公式无法使用。下面使用卡尔曼滤波算法来计算P^*[13-14]。

$ P_{i+1}=P_{i}+S_{i+1} a_{i+1}\left(b_{i+1}^{\mathrm{T}}-a_{i+1}^{\mathrm{T}} P_{i}\right) $

(19)

$ S_{i+1}=S_{i}-\frac{S_{i} a_{i+1} a_{i+1}^{\mathrm{T}} S_{i}}{1+a_{i+1}^{\mathrm{T}} S_{i} a_{i+1}}, \quad(i=0, 1, \cdots, m-1) $

(20)

a_i^T是矩阵A的第i行向量，b_i^T是矩阵B的第i个元素(若矩阵B的列数不为1，即对应多输出系统，则表示第i行元素)。S_i称作协方差矩阵。最终得到的P_m即为P^*。初始值P₀=0, S₀=γI, γ表示一个极大的正数，I是维数为m×m的单位矩阵。

2.4 模糊结构识别

此模糊结构的输入是62个测压点的压力或压力系数，选取其中k个点的值作为输入，输出为迎角、侧滑角、速度和高度。

将输入记做x₁, …, x_r, …, x_k，将输出记做y_α, y_β, y_v, y_h。以对迎角建模为例，对侧滑角、来流速度和海拔高度的建模相同。

(1) 将所有样本分为两组，一组为训练样本，用于参数识别，一组为测试样本，用于检测建模精度。给定训练样本的R²最小值R_min²。

(2) 初始时，对每个输入选取N_begin个隶属函数。则初始隶属函数个数的结构是(N_begin, …, N_begin, …, N_begin)。然后对其中某一个输入的隶属函数个数加1，其他输入的隶属函数个数保持不变，这样可以建立k个隶属函数个数不同的结构。

(3) 计算所有k个结构的R_train²和R_test²，将其中R_test²的最大值对应的结构作为最优结构，记此时的结构寻找次数N_S=1。

(4) 重复上述寻找过程，当满足R_train²(N_S)>R_min²并且R_test²(N_S)>R_test²(N_S+1)>R_test²(N_S+2)时，结构识别结束。结构识别过程中可能会遇到R_test²比较接近且来回波动的情况，此时可适当放宽要求，当满足R_test²(N_S)>R_test²(N_S+1), R_test²(N_S)>R_test²(N_S+2)时即可。

(5) 得到最优隶属函数个数的结构之后，用所有样本数据重新进行参数识别，得到最终模糊结构参数，建模结束。

3 建模结果分析

试验模型表面共布置62个测压点。选取其中若干个测压点的压力或压力系数作为输入，以迎角、侧滑角、风速和高度作为输出。

对迎角和侧滑角建模时选择压力系数作为输入。试验风速最大为70 m/s。CFD计算时选取的最大高度为6000 m。在此区间内在姿态角不变的情况下，风速和高度的改变只引起压力系数的微小变化，几乎可以忽略。所以对风速和高度的建模选取压力作为输入。

考虑此模型的非线性复杂程度，通常选择不少于3个点作为建模输入。建模点数量增加时，参数识别消耗的时间以指数程度增加，平衡考虑建模的精度和时间效益，建模点数量一般不多于7个。

从62个测压点选择若干个点建模，选取方式很多，可以适当探讨建模点位置选取与建模精度的关系，为选点位置的优化设计提供参考。

3.1 建模结果

以选取4个点(2, 6, Z4, Y4)建模为例，结果如下。2、6编号的点位于模型头部上下区域，压力变化对迎角的改变敏感，Z4、Y4编号的点位于头部左右区域，压力变化对侧滑角的改变敏感。

图 8中横坐标为测压孔2的压力系数C_p。

表 2和图 8~图 12展示了选择4个点(2, 6, Z4, Y4)建模的结果。模糊输出与试验值的对比达到了很高的精度。

在对迎角与侧滑角及风速建模时，此时有300多组样本，在样本较多时，选用三角隶属函数与梯形隶属函数差异不大。对高度进行建模时，样本较少，有96组。在本文使用的样本及样本数量较少的情况下, 梯形隶属函数的精度更高。

对高度建模时，加入了使用变参数钟形隶属函数的建模结果。使用变参数隶属函数建模时，分别使用了梯度下降法和最小二乘估计法识别隶属函数参数与内部函数参数, 从表 2的数据及图 11和图 12的对比可以看出，精度相比三角隶属函数与梯形隶属函数高很多，且建模速度很快。使用卡尔曼滤波的最小二乘估计法更适用于小样本的情况。图 13和图 14展示了模型训练前后针对测压孔2的隶属函数曲线。

3.2 建模点数量与位置的优化选择

以对迎角建模为例，分别选择不同位置的3个点、4个点或5个点用于建模。在建模点数目相同时，使用相同的收敛因子、迭代次数和相同的隶属函数个数和形状(三角隶属函数)，记录精度R²，比较筛选出最优的建模位置。

建模结果如下，表格数据以R²降序排序。

表 3~表 6中选点是编写程序从62个测压孔中随机选取的，精度均很高，相差较小，从R²降序结果观察比对选点的位置区域，考虑到试验时使用的测压仪器在模型不同区域的测量精度不同，通常选取模型头部附近的测压点进行建模。

从表 3~表 5还可以看出，选择5个点的建模精度稍高于4个点的建模精度，4个点的建模精度稍高于3个点的建模精度。不过随着建模点数的增加，建模速度大大降低。表中数据表明选择4个点或5个点的建模精度能够很好地满足工程需求。

计算时发现存在一些精度非常低的选点方式(即处于同一水平面或铅垂面的点)，在实际应用中应当避免。

4 结论

本文研究了对嵌入式大气数据传感系统的模糊逻辑建模方法。以模型表面若干测压点的压力或压力系数作为模糊系统的输入，以迎角、侧滑角、来流速度和海拔高度作为输出，探讨了在使用不同隶属函数和参数识别方法情况下的建模结果，讨论分析了建模点数量与位置的优化选择。结论如下：

1) 针对本文的模糊逻辑系统，使用自适应形状参数的隶属函数的建模精度要高于使用固定形状的三角隶属函数和梯形隶属函数。且在建模时分别使用梯度下降法和最小二乘估计法来识别隶属函数参数和内部函数参数，不论是从精度还是速度上来说都优于单独使用梯度下降法来识别内部函数参数。

2) 在建模点数量与位置的优化选择方面，考虑传感器在模型不同区域的测量精度，通常选择头部附近的测压数据进行建模。对迎角的建模应当避免选择处于同一水平面的测压点，对侧滑角避免选择处于同一铅垂面的测压点。结果表明选择4个点或5个点时的精度足以满足实际工程需求，且建模速度较快。对于复杂外形飞行器，可以适当增加建模点数和隶属函数个数。

受限于风洞试验的高昂成本，本文建立的模糊逻辑模型主要基于低速情况下的流动数据。想要拓宽模型的应用范围，需要更多的不同流动情况下的数据，未来的方向是将模型的研究与验证拓展到更高速的流动场景下。