目前建模与仿真技术已成为人们认识和改造现实世界的重要手段.由于仿真是一种基于模型的活动, 仿真模型是否可信成为用户十分关注的问题.验证是仿真模型可信度评估的重要步骤^[1], 包含概念模型验证和结果验证.仿真结果验证最直接而有效的方法是, 在相同输入条件下度量仿真输出与参考输出数据之间的一致性程度.然而, 针对复杂系统建立的仿真模型往往具有不确定性、多元异类(动、静态)输出, 且各输出变量间可能存在相关性, 此时, 若仍采用传统仿真结果验证方法将导致验证结果不准确.因此, 考虑相关性及不确定性的多元输出仿真结果验证是需要重点研究的问题.

由于仿真模型和参考系统的输入及模型参数通常含有不确定性, 加之仿真模型运行和实际实验过程中引入的不确定性因素和误差, 导致仿真模型和参考系统的输出为随机变量或不确定的时间序列^[2].考虑不确定性的影响, 静态输出结果验证方法的研究多数集中在概率框架, 形成了以参数估计^[3-5]、假设检验^[6-7]、贝叶斯因子^[8-10]、证据距离^[11]、概率分布差异法^[12-13]为代表的5种解决方案.其中, Oberkampf等针对参考数据稀疏的情况, 采用插值和回归分析的方法估计参考输出的均值和标准差, 并与仿真输出的相应统计量进行对比, 得到置信区间形式的验证结果^[5]; 同时, 假设检验和贝叶斯因子方法在静态仿真输出结果验证中的应用日趋完善, Jiang等将贝叶斯区间假设检验应用于模型的分等级评估中^[14]; 考虑到固有和认知不确定性的影响, 文献[11]采用证据理论对动静态输出进行描述, 并引入证据距离度量仿真和参考输出的一致性; Ferson等提出了概率分布差异与u-pooling相结合的方法, 用于处理稀疏参考数据情况下的单输出仿真结果验证问题^[12], 该方法以其原理简单、可操作性强等优点得到了广泛应用.

考虑不确定性的同时, 复杂仿真模型可能存在多个输出变量的情况, 且各变量间可能存在函数或相关关系.在单变量静态输出结果验证方法的基础上, 针对多元输出仿真结果验证方法的研究取得了一定的进展.例如, Rebba等最先提出了假设检验、贝叶斯因子与协方差相结合的多元输出结果验证方法, 并引入了非正态验证数据转化为正态数据的方法以满足假设检验的条件^[15]; Jiang等将区间贝叶斯因子方法进行推广, 将其应用至多元静态输出结果验证问题中^[16]; Zhan等提出了基于概率主成分分析(Probabilistic principal component analysis, PPCA)与贝叶斯因子相结合的方法, 用于解决带有不确定性和相关性的多元动态输出结果验证问题^[17].同时, Li等将概率积分变换(Probability integral transformation, PIT)与概率分布差异法相结合, 将多变量累积概率分布转化为单变量概率分布的形式, 采用概率分布差异法计算仿真和参考输出累积概率分布的差异^[18]; Zhao等分别计算仿真和参考输出与相应总体分布的马氏距离, 进而得到仿真和参考输出马氏距离的累积概率分布, 并应用概率分布差异法计算两者的差异^[19].

从单变量验证到多变量验证, 各变量间的关系是研究的重点, 现有验证方法存在验证结果不够准确和全面的问题.利用传统单变量验证(或结合多种预处理方法)对各变量进行分别验证再综合, 一是对带有相关性的多个验证结果进行加权综合导致最终验证结果不够准确; 二是未考虑多变量间相关性的验证将导致验证结果不全面.此外, 对复杂仿真模型进行结果验证, 其输出变量间的关系不够明确.此时需首先明确输出变量间的关系(独立/函数/相关关系)再进行验证, 现有的多变量验证方法仅适用于变量间关系已知的情况.同时, 多变量验证方法均利用协方差矩阵度量变量间的相关性, 这对非线性等其他相关关系将不再适用, 导致多变量间相关性度量不够准确.

为解决上述问题, 提出基于变量选择和概率分布差异相结合的多元输出仿真结果验证方法, 对具有不确定性的多元异类输出进行联合验证.第1节对多元输出结果验证问题进行描述与分析, 指出现有方法存在的问题; 第2节给出多元静、动态输出的相关性分析及变量选择方法; 第3节提出基于数据特征提取和联合概率分布差异的多元输出仿真结果验证方法; 第4节给出应用实例与对比实验结果; 第5节给出本文结论.

1.   问题描述与分析

用$ {S} $表示系统, $ {S}_{{s}} $和$ {S}_{{r}} $分别表示仿真模型和参考系统, 用$ {\pmb{A}_{{s}}} = \{\pmb{{a}}_{{s1}} $, $ {\pmb{a}}_{{s2}}, \cdots, {\pmb{a}}_{{s}{p}}\} $和$ {\pmb{A}_{{r}}} = \{{\pmb{a}}_{{r1}} $, $ {\pmb{a}}_{{r2}}, \cdots, {\pmb{a}}_{{r}{p}}\} $分别表示仿真模型和参考系统的输入变量集, $ {p} $为输入变量个数, $ {\pmb{Y}\!_{{s}}} = \{{\pmb{y}}_{{s1}} $, $ {\pmb{y}}_{{s2}}, \cdots, {\pmb{y}}_{{s}{m}} $}和$ {\pmb{Y}\!_{{r}}} = \{{\pmb{y}}_{{r1}} $, $ {\pmb{y}}_{{r2}}, \cdots, {\pmb{y}}_{{r}{m}} $}分别表示仿真模型和参考系统的输出变量集, m为输出变量个数, 多元异类输出集$ {\pmb{Y}}\!_{{s}} $、$ {\pmb{Y}}\!_{{r}} $中的静态输出表示为随机变量, 动态输出表示为多个时间序列集合的形式.假设$ {\pmb{y}}_{{s}{i}} $、$ {\pmb{y}}_{{s}{j}} $分别为仿真模型的某一动态和静态输出, 则有

式中, $ i, j\in [1, m] $, 且$ i\ne j $; $ \it N $为时间序列的长度; $ {{t}_{1}}, {{t}_{2}}, \cdots , {{t}_{N}} $表示时间序列的时刻点; 考虑不确定性的影响, 需要进行多次仿真和实际实验, $ \it n $为重复实验次数.用$ C\left( {{{\pmb{Y}}}\!_{{s}}}, {{{\pmb{Y}}}\!_{{r}}} \right) $表示在$ {{{\pmb{A}}}_{{s}}} = {{{\pmb{A}}}_{{r}}} $时, $ {{{\pmb{Y}}}\!_{{s}}} $相对于$ {{{\pmb{Y}}}\!_{{r}}} $的一致性程度, 且$ C\left( {{{\pmb{Y}}}\!_{{s}}}, {{{\pmb{Y}}}\!_{{r}}} \right)\in \left( 0, 1 \right] $.当$ {{{\pmb{Y}}}\!_{{s}}} $与$ {{{\pmb{Y}}}\!_{{r}}} $完全一致, 则有$ C\left( {{{\pmb{Y}}}\!_{{s}}}, {{{\pmb{Y}}}\!_{{r}}} \right) = 1 $; 当$ {{{\pmb{Y}}}\!_{{s}}} $相对于$ {{{\pmb{Y}}}\!_{{r}}} $的一致性程度越差, 表示仿真模型越不可信, 则有$ C\left( {{{\pmb{Y}}}\!_{{s}}}, {{{\pmb{Y}}}\!_{{r}}} \right)\to 0 $^[1].

假设$ {\pmb{Y}}\!_{{sJ}}^{{}}\in {{\bf {R}}^{{{n}_{{s}}}\times m}} $与$ {\pmb{Y}}\!_{{rJ}}^{{}}\in {{\bf {R}}^{{{n}_{{r}}}\times m}} $为多元仿真模型和参考系统的静态输出变量, $ {{n}_{{s}}} $、$ {{n}_{{r}}} $表示仿真和实际实验的重复运行次数.针对带有相关性的多元静态输出结果验证方法主要有:

1) 基于假设检验和马氏距离相结合的方法.文献[16]给出基于似然比检验和马氏距离相结合的验证方法, 得到最终一致性结果.

2) 基于主成分分析的方法^[17].对$ {\pmb{Y}}\!_{{sJ}}^{{}} $和$ {\pmb{Y}}\!_{{rJ}}^{{}} $进行降维, 去除变量间相关性, 得到新的输出变量$ {\pmb{Y}}\!_{{sJ}}^{{new}} = \left[ y_{{sJ}1}^{{new}}, y_{{sJ2}}^{{new}}, \cdots, y_{{sJ}\eta }^{{new}} \right] $和$ {\pmb{Y}}\!_{{rJ}}^{{new}} = \left[ y_{{rJ}1}^{{new}}, y_{{rJ2}}^{{new}}, \cdots, y_{{rJ}\eta }^{{new}} \right] $, $ \eta \le m $为主成分的个数, 而后采用现有静态输出验证方法对若干主成分进行逐一验证并综合得到最终验证结果.

3) 基于概率分布差异的方法^[18].分别计算m维$ {\pmb{Y}}\!_{{sJ}}^{{}} $和$ {\pmb{Y}}\!_{{rJ}}^{{}} $的联合累积概率分布(Cumulative distribution function, CDF)函数并作差, 获得仿真和参考输出数据的差异, 得到$ \left[ 0, +\infty \right) $的误差度量结果.

针对带有相关性的多元动态输出结果验证问题, 常用方法为基于数据特征和主成分分析相结合的方法^[17].首先分别提取动态输出数据的特征矩阵, 而后采用基于主成分分析的多元静态输出验证方法获得最终验证结果.针对上述多元输出仿真结果验证方法进行分析, 存在以下问题需要进一步研究:

1) 复杂仿真模型常存在多元输出变量间的相关或独立关系未知的情况, 目前方法均是在变量关系已知的前提下进行研究, 存在一定局限性;

2) 利用主成分分析获取的多元输出变量的主成分是线性变换后的结果, 被提取主成分所代表的变量含义不够明确, 同时对多元输出变量进行降维将导致验证信息丢失, 使验证结果不够准确和全面;

3) 采用协方差矩阵度量变量相关性, 需假设变量样本服从正态分布, 且仅能描述多元输出变量间的线性关系, 无法度量变量间非线性等其他相关关系, 进而导致变量间相关性度量不准确;

4) 基于联合概率分布差异法可直接度量多元静态输出变量间的差异, 需要已知变量间的独立或相关关系, 同时, 处理多元动态输出存在局限, 得到的差异度量结果无法刻画仿真模型的可信度.

为解决上述问题, 可采用基于变量选择和概率分布差异相结合的多元输出仿真结果验证方法, 考虑不确定性的影响, 对选取到具有相关性的多变量进行联合验证.首先, 引入变量选择方法分别对$ {{{\pmb{Y}}}\!_{{s}}} $、$ {{{\pmb{Y}}}\!_{{r}}} $进行相关性分析, 提取相关变量子集(又称相关变量组, 子集中各变量是相关的, 各子集中变量数的和为输出变量总数), 进而得到多个独立的变量子集; 同时, 提取相同变量子集中多变量的数据特征, 对于静态输出选取数据本身作为变量特征, 对于动态输出选取距离、形状以及频谱特征; 而后计算变量子集中多个变量关于某特征的联合CDF差异, 并将其转化为可信度; 最后将多个变量子集关于若干数据特征的一致性与多个动态输出均值曲线的一致性进行综合得到仿真模型可信度.

2.   多元输出相关变量选择方法

为明确复杂仿真模型中多元输出变量间的独立或相关关系, 引入数据挖掘领域的相应方法对多变量进行相关性分析, 进而提取相关变量子集.本文仅考虑同种类型(静态或动态)输出变量间的相关性, 利用分形维数和互信息方法分别对静、动态输出变量进行相关性分析.

2.1   基于分形维数的静态相关变量选择方法

对随机变量的相关性分析集中于Pearson相关系数, 它仅能度量变量的线性关系, 并对变量间强相关性较敏感, 其结果受奇异值的影响较大, 无法适应具有非线性、不确定性以及非正态分布的数据集.其他一些相关系数, 如Kendall系数、Spearman系数等虽可以描述非线性相关关系, 但却不能完整地刻画变量间的相关性结构.此外, 数据挖掘领域常用的变量选择方法, 如奇异值分解法(Singular value decomposition, SVD)、主成分分析法(Principal component analysis, PCA)、基于神经网络的方法(Neural networks, NN)、基于k-邻近方法(K-nearest neighbor, KNN)、基于决策树的方法(Decision tree, DT)、基于贝叶斯网络的方法(Bayesian network, BN)以及基于分形维数的方法(Fractal dimension, FD)等, 具有不同的特点, 对其进行对比分析如表 1所示.

表 1  常用变量选择方法对比

Table 1  Comparison of general variable selection methods

变量选择方法	是否为原变量集的子集	是否支持非线性相关关系	个体决策所占比例	是否需要训练样本集	运行速度与变量个数的关系
SVD	否	否	是	否	线性增长
PCA	否	否	是	否	线性增长
KNN	是	是	是	是	指数增长
DT	是	是	是	是	指数增长
BN	是	是	是	是	指数增长
FD	是	是	是	否	线性增长

下载: 导出CSV 

| 显示表格

由表 1可知, SVD和PCA方法得到的变量子集失去了其原有的含义, 且只能对具有线性相关性的变量集进行分析; 而基于机器学习的方法需要训练样本集作为支撑, 其运行速度受到变量个数的影响较大, 导致变量个数较多时运行速度较慢; 而基于分形维数的方法不仅能够度量线性相关性, 还能度量非线性等其他相关关系, 具有不需要训练样本集和运行速度快等优点.因此, 本文引入基于分形维数^[20]的方法对多元输出变量进行分析, 提取相关变量子集.假设$ {\pmb{Y}}\!_{{sJ}}^{{}} $、$ {\pmb{Y}}\!_{{rJ}}^{{}} $为仿真和参考多元静态输出变量集, 以$ y_{{sJ}}^{i} $, $ i\in \left[ 1, m \right] $为例给出$ {\pmb{Y}}\!_{{sJ}}^{{}} $的相关变量子集提取步骤如下.

步骤1. 根据自相似性原理计算$ y_{{sJ}}^{i} $的局部固有维度$ pD\left( \cdot \right) $:

式中, r表示将$ y_{{sJ}}^{i} $划分成$ {{2}^{\upsilon }} $个相等大小区间的长度, $ \upsilon $为划分深度, $ {{C}_{a, i}} $表示$ y_{{sJ}}^{i} $中落入第i个区间的样本个数;

步骤2. 设$ c = 1 $, 移除$ {\pmb{Y}}\!_{{sJ}} $中$ pD(y_{{sJ}}^{i})<\xi $的变量$ y_{{sJ}}^{i} $, $ \xi $为预定义的固有维度阈值, 排除的变量为独立变量, 并按照$ pD\left(\cdot \right) $大小将$ y_{{sJ}}^{i} $进行降序排列, 形成新变量集$ {\pmb{Y}}\!_{{sJ}}^{\prime} $, 其变量个数为$ {m}' $;

步骤3. 计算$ pD\left( \left\{ y_{{sJ}}^{1} \right\} \right), pD\left( \left\{ y_{{sJ}}^{1}, y_{{sJ}}^{2} \right\} \right), \cdots $, 直到$ \left| pD\left( \left\{ y_{{sJ}}^{1}\cdot y_{{sJ}}^{k} \right\} \right)-pD\left( \left\{ y_{{sJ}}^{1}\cdot y_{{sJ}}^{k-1} \right\} \right) \right|<\xi \cdot pD\left( y_{{sJ}}^{k} \right) $, $ k = 1, 2, \cdots, {m}' $;

步骤4. 若$ k = {m}' $且$ \left| pD\left( \left\{ y_{{sJ}}^{1}\cdot y_{{sJ}}^{k} \right\} \right)- \right. $ $ \left.pD\left\{ \left( y_{{sJ}}^{1} \cdot y_{{sJ}}^{k-1} \right\} \right) \right|\ge \xi \cdot pD\left( y_{{sJ}}^{k} \right) $, 则算法结束;

步骤5. 设相关性变量超集$ \xi S{{G}_{c}} = \left\{ y_{{sJ}}^{1}, \cdots , \right. $ $ \left. y_{{sJ}}^{k} \right\} $, 并提取$ \xi S{{G}_{c}} $中的相关变量子集$ \xi {{G}_{c}} $和相关变量基$ \xi {{B}_{c}} $, 具体算法见文献[20], 并设循环变量$ j = k+1 $;

步骤6. 若$ \left| pD\left( \xi {{B}_{c}}\bigcup \left\{ y_{{sJ}}^{j} \right\} \right)- \right.\left. pD\left( \xi {{B}_{c}} \right) \right|<\xi \cdot pD\left( y_{{sJ}}^{j} \right) $, 则执行下一步, 否则转至步骤8;

步骤7.        对于$ \xi {{B}_{c}} $中的每个变量$ y_{{sJ}}^{b} $, 若$ \Big| pD\left( \xi {{B}_{c}} \bigcup \left\{ y_{{sJ}}^{j} \right\} \right)-pD\Big( \left( \xi {{B}_{c}}-\left\{ y_{{sJ}}^{b} \right\} \right)\bigcup $ $ \left\{ y_{{sJ}}^{j} \right\} \Big) \Big|<\xi \cdot pD\left( y_{{sJ}}^{b} \right) $和$ \Big| pD\left( \xi {{B}_{c}}\bigcup \left\{ y_{{sJ}}^{b} \right\} \right)- $ $ pD\left( \left( \xi {{B}_{c}}-\left\{ y_{{sJ}}^{b} \right\} \right) \bigcup \left\{ y_{{sJ}}^{j} \right\} \right) \Big|\ge \xi \cdot pD\left( y_{{sJ}}^{j} \right) $同时成立, 则将$ y_{{sJ}}^{j} $加入$ \xi {{G}_{c}} $;

步骤8. 执行$ j\leftarrow j+1 $, 若$ j>{m}' $, 则转至下一步, 否则转至步骤6;

步骤9. 移除$ {\pmb{Y}}_{{sJ}}^{\prime} $中$ \xi {{G}_{c}}-\xi {{B}_{c}} $的变量, 并输出相关变量子集$ \xi {{G}_{c}} $和相关变量基$ \xi {{B}_{c}} $;

步骤10. 执行$ c\leftarrow c+1 $, 并转至步骤3.

通过上述步骤提取$ {\pmb{Y}}\!_{{sJ}} $和$ {\pmb{Y}}\!_{{rJ}} $的相关变量子集$ {\pmb{G}}_{{sJ}}^{i} $、$ {\pmb{G}}_{{rJ}}^{j} $如下.

式中, $ {{\beta }_{{s}}} $、$ {{\beta }_{{r}}} $分别为提取$ {\pmb{Y}}\!_{{sJ}}^{{}} $和$ {\pmb{Y}}\!_{{rJ}}^{{}} $相关变量子集的个数, $ {{m}_{{s}i}} $、$ {{m}_{{r}j}} $分别为$ {\pmb{G}}_{{sJ}}^{i} $、$ {\pmb{G}}_{{rJ}}^{j} $中变量的个数, 且有$ {{m}_{{s}1}}+{{m}_{{s}2}}+\cdots +{{m}_{{s}i}} = {{m}_{{r}1}}+{{m}_{{r}2}}+\cdots $ +$ {{m}_{{r}j}} = m $.

2.2   基于互信息的动态相关变量选择方法

与随机变量不同, 多元动态输出变量与时间有关, 其相关性分析与变量选择需从时间序列整体的角度进行分析.一些传统的随机变量相关性分析方法对于多元动态变量同样适用, 例如Pearson系数、Kendall系数、Spearman系数等, 但无法用于动态输出变量具有多个样本(时间序列)的情况.此外, 一些统计学分析方法, 如Granger因果关系分析^[21]、典型相关分析^[22]、Copula分析^[23]、灰色关联分析^[24]以及互信息分析^[25]等同样能够用于多变量的相关性分析. Granger因果关系分析只能定性地分析变量间的因果关系, 而无法得到定量的结果; 典型相关分析对观测值的顺序不会做出响应, 因此无法解决时间序列问题; Copula分析需要建立在对边缘分布的合理假设之上, 使其应用受到限制; 灰色关联分析仅从形状相关性的角度对时间序列进行分析, 其相关性度量不够全面.

基于互信息的相关性分析方法能够度量动态输出变量间任意类型的关系, 互信息以信息熵为理论基础, 它能够度量变量取值的不确定性程度, 进而描述变量的信息含量大小^[26], 通常用于多种类型时间序列的特征提取和结构化预测^[27].然而, 互信息同样存在不能完整刻画变量集相关性结构的缺点, 因此本文引入类可分性和变量可分性提取多元动态输出的相关变量子集^[28].假设$ {\pmb{Y}}\!_{{sD}}^{{}} $、$ {\pmb{Y}}\!_{{rD}}^{{}} $为仿真和参考多元动态输出变量集, 同样以$ {\pmb{Y}}\!_{{sD}}^{{}} $为例, 给出变量选择步骤如下.

步骤1. 计算$ {\pmb{Y}}\!_{{sD}}^{{}} $的$ m\times m $维互信息矩阵, 具体算法见文献[26];

步骤2. 分别计算每一维变量的类间离散度$ {{\Omega }_{{b}i}} $和类内离散度$ {{\Omega }_{{w}i}} $:

式中, $ {{C}_{{sam}}} $为样本类别总数, $ {{q}_{i}} $为属于第i类的样本个数, $ \mu = ({{1}}/{{{n}_{{s}}}})\;\sum\nolimits_{i = 1}^{{{n}_{{s}}}}{y_{{sD}}^{i}} $, $ {{\mu }_{i}} = ({{1}}/{{{q}_{i}}})\;\sum\nolimits_{i = 1}^{{{q}_{i}}}{y_{{sD}}^{i}} $.按照每个变量的类可分离性大小, 进行变量排序:

步骤3.  取$ {{J}_{i}} $值最大的变量为变量子集$ {\pmb{G}}_{{sD}}^{i} $的第一个变量;

步骤4.  选择使下式最大的变量为$ {\pmb{G}}_{{sD}}^{i} $的下一个变量:

式中, $ \left| {\pmb{G}}_{{sD}}^{i} \right| $为子集$ {\pmb{G}}_{{sD}}^{i} $的变量个数, $ {{\mu }_{ko}} $为子集$ {\pmb{G}}_{{sD}}^{i} $中属于第$ o $类的第$ k $个变量的均值;

步骤5. 当$ \left| {\pmb{G}}_{{sD}}^{i} \right| = \varepsilon $, 则算法终止, 其中, $ \varepsilon $为预设值, 否则转至步骤4.通过上述步骤得到相关变量子集$ {\pmb{G}}_{{sD}}^{i} $、$ {\pmb{G}}_{{rD}}^{j} $如下:

式中, $ i = 1, 2, \cdots , {{\beta }_{{s}}} $, $ j = 1, 2, \cdots , {{\beta }_{{r}}} $. $ {{\beta }_{{s}}} $、$ {{\beta }_{{r}}} $的含义与式(3)相同.

3.   多元输出仿真结果验证方法

考虑不确定的影响, 若对每一时刻的多元动态输出变量进行分析势必导致维数爆炸.为此提出基于特征的验证方法, 首先提取用户关注的输出数据特征, 而后计算每个特征下多变量联合概率分布的差异, 并将其转化为可信度结果, 最后综合多个验证结果得到模型可信度.

3.1   多元输出数据特征提取方法

针对于静态输出, 选取数据本身作为变量特征.假设$ {{{\pmb{Y}}}\!_{{s}}}\in {{\bf {R}}^{{{n}_{{s}}}\times m}} $与$ {{{\pmb{Y}}}\!_{{r}}}\in {{\bf {R}}^{{{n}_{{r}}}\times m}} $为多元仿真和参考静态输出变量, 其数据特征描述为

对于动态输出$ {{y}_{{s}ij}}(t) $、$ {{y}_{{r}ik}}(t) $, $ i = 1, 2, \cdots , m $, $ j = 1, 2, \cdots , {{n}_{{s}}} $, $ k = 1, 2, \cdots , {{n}_{{r}}} $则选取$ {{n}_{{s}}} $、$ {{n}_{{r}}} $次系统运行得到的输出均值曲线$ {{\bar{y}}_{{s}i}}(t) $、$ {{\bar{y}}_{{r}i}}(t) $作为基准, 与每次实验得到的输出曲线进行对比, 求取相应的均值曲线的第$ l $个特征$ e_{{s}ij}^{l} $、$ e_{{r}ik}^{l} $:

式中, $ {{L}_{i}} $为第i个输出的特征数, $ {{\Phi }_{l}}\left( \cdot \right) $为第l个特征度量模型.

提取动态输出特征前, 需要先对动态输出进行归类^[1].以第j个动态输出的第i次实现$ {{y}_{ij}}\left( t \right) $为例, 其对应的时间变化序列为$ \left[ {{t}_{1}}, {{t}_{2}}, \cdots , {{t}_{N}} \right] $.则定义$ {{y}_{ij}}\left( t \right) $随时间变化的频率为:

式中, $ {{F}_{ij}}\ge 0 $为$ {{y}_{ij}}\left( t \right) $的变化频率; $ \Delta {{y}_{ij}}\left( {{t}_{k}} \right) = {{y}_{ij}}\left( {{t}_{k+1}} \right)-{{y}_{ij}}\left( {{t}_{k}} \right) $; $ \Delta {{t}_{k}} = {{t}_{k+1}}-{{t}_{k}} $; $ \left| {{{\bar{y}}}_{ij}} \right| = {\sum\nolimits_{k = 1}^{N}{\left| {{y}_{ij}}\left( {{t}_{k}} \right) \right|}}/{N}\;\ne 0 $.给定$ {{F}_{0}} $为判断时间序列变化快慢的临界值, 若$ {{F}_{ij}}\ge {{F}_{0}} $, 则认为$ {{y}_{ij}}\left( t \right) $为速变数据, 否则为缓变数据.

为刻画不确定性对系统输出的影响, 从距离和形状两方面提取缓变数据的特征.在前期工作^[29]的基础上, 给出第j个仿真输出的第i次实现$ {{y}_{{s}ij}}\left( t\right) $与其均值曲线$ {{\bar{y}}_{{s}j}}(t) $的距离和形状差异$ e_{{sd}}^{ij} $、$ e_{{sc}}^{ij} $的度量公式如下.

式中, $ z\left( t \right) = y_{{s}ij}^{{}}\left( t \right)-\bar{y}_{{s}j}^{{}}\left( t \right) $, $ t = {{t}_{1}}, {{t}_{2}}, \cdots , {{t}_{N}} $, $ \bar{z} = {\sum\nolimits_{t = {{t}_{1}}}^{{{t}_{N}}}{z\left( t \right)}}/{N}\; $.另外, 选取谱密度特征度量速变数据$ {{y}_{{s}ij}}\left( t \right) $与相应均值曲线$ {{\bar{y}}_{{s}j}}(t) $的差异$ e_{{sh}}^{ij} $, 定义如下.

式中, $ e_{{sh}}^{ij} $表示速变数据$ {{{\pmb{y}}}_{{s}ij}} $与$ {{\bar{{\pmb{y}}}}}_{{s}j} $的谱密度差异; M表示$ {{y}_{{s}ij}}\left( t \right) $和$ {{\bar{y}}_{{s}j}}(t) $转换至频域中的点数; $ \gamma $表示通过相容性检验的点数.根据式(2)$ \sim $(7)得到$ {{{\pmb{Y}}}\!_{{s}}} $的第i个相关变量子集$ {{{\pmb{G}}}_{{s}i}} $、$ {{{\pmb{Y}}}\!_{{r}}} $的第j个相关变量子集$ {{{\pmb{G}}}_{{r}j}} $关于第l个特征的差异度量矩阵分别为

式中, $ i = 1, 2, \cdots , \alpha $, $ j = 1, 2, \cdots , \beta $.若$ {{{\pmb{G}}}_{{s}i}} $和$ {{{\pmb{G}}}_{{r}j}} $均为静态输出变量子集, 则$ {{L}_{i}} = 1 $; 若$ {{{\pmb{G}}}_{{s}i}} $和$ {{{\pmb{G}}}_{{r}j}} $均为缓变输出变量子集, 则$ {{L}_{i}} = 2 $; 若$ {{{\pmb{G}}}_{{s}i}} $和$ {{{\pmb{G}}}_{{r}j}} $均为速变输出变量子集, 则$ {{L}_{i}} = 1 $.需要说明的是, 在某些特殊仿真应用中, 除了上述特征外, 通常还需关注数据本身的一些特征, 例如控制系统阶跃响应中的上升时间、超调量以及稳态误差, 位置数据中的变化趋势, 测量数据中的噪声等.在进行实际验证中, 特征矩阵包含两部分内容, 一部分为上文给出的数据特征, 另一部分为根据具体领域知识确定的数据特征.

3.2   基于概率分布差异的特征一致性分析方法

以$ {\pmb{E}}_{{s}i}^{l} $为例进行分析, 用$ \upsilon $维随机变量$ {{x}_{1}}, {{x}_{2}}, \cdots , {{x}_{\upsilon }} $代替其列向量$ \left[ {\pmb{e}}_{{s}i1}^{l}, {\pmb{e}}_{{s}i2}^{l}, \cdots , {\pmb{e}}_{{s}i{{m}_{{s}i}}}^{l} \right] $, $ \upsilon = {{m}_{{s}i}} $, 采用多维随机变量概率分布定义$ {{x}_{1}}, {{x}_{2}}, \cdots , {{x}_{\upsilon }} $的联合CDF:

将$ \upsilon $维空间划分为等尺寸的$ {{\rho }^{\upsilon }} $个区域, 遍历$ \upsilon $维变量的$ \rho $个取值区间, 若$ {{x}_{1}}<X_{1}^{0} $, $ {{x}_{2}}<X_{2}^{0}, \cdots, {{x}_{\upsilon }}<X_{\upsilon }^{0} $, 则$ F\left( {{x}_{1}}, {{x}_{2}}, \cdots, {{x}_{\upsilon }} \right) = 0 $; 若$ {{x}_{1}}<X_{1}^{k} $, $ {{x}_{2}}<X_{2}^{k}, \cdots, {{x}_{\upsilon }}<X_{\upsilon }^{k} $, 则$ F\left({{x}_{1}}, {{x}_{2}}, \cdots , {{x}_{\upsilon }} \right) = {k}/{{{\rho }^{\upsilon}}} $等.如果变量集$ {{x}_{1}}, {{x}_{2}}, \cdots , {{x}_{\upsilon }} $在第k个区间内的样本量为1, 则F在$ {{x}_{1}}, {{x}_{2}}, \cdots , {{x}_{\upsilon }} $点的跳跃度为$ {1}/{{{\rho }^{\upsilon }}} $, 如果变量集$ {{x}_{1}}, {{x}_{2}}, \cdots , {{x}_{\upsilon }} $在第k个区间内有$ \varepsilon $个样本, 则F在$ {{x}_{1}}, {{x}_{2}}, \cdots , {{x}_{\upsilon }} $点的跳跃度是$ {\varepsilon }/{{{\rho }^{\upsilon }}} $.给出$ {\pmb{E}}_{{s}i}^{l} $和$ {\pmb{E}}_{{r}j}^{l} $联合CDF间的差异如下.

为计算联合CDF的差异$ D\left( {{F}_{{r}}}, {{F}_{{s}}} \right) $, 可将上式改写为下面积分之差的形式:

式中$ {\pmb{x}} = \left[ {{x}_{1}}, {{x}_{2}}, \cdots , {{x}_{\upsilon }} \right] $.假设分别用$ {{\hat{I}}_{{s}}} $和$ {{\hat{I}}_{{r}}} $估计$ {{I}_{{s}}} $和$ {{I}_{{r}}} $, 并用$ \hat{D} = {{\hat{I}}_{{s}}}-{{\hat{I}}_{{r}}} $估计$ \hat{D} $, 则$ \hat{D} $的方差为:

显然, 在$ {Var}\left( {{{\hat{I}}}_{{s}}} \right) $和$ {Var}\left( {{{\hat{I}}}_{{r}}} \right) $确定后, $ {{\hat{I}}_{{s}}} $和$ {{\hat{I}}_{{r}}} $的正相关度越高, 则$ \hat{D} $的方差越小.本文采用重要性抽样法估计$ {{I}_{{s}}} $和$ {{I}_{{r}}} $, 即改写D为

其中, $ {\pmb{x}} = {{x}_{1}}, {{x}_{2}}, \cdots , {{x}_{\upsilon }} $, $ {{g}_{{s}}}\left( {\pmb{x}} \right) $、$ {{g}_{{r}}}\left( {\pmb{x}} \right) $是两个密度函数, $ {{H}_{{s}}}\left( {\pmb{x}} \right) = {{{F}_{{s}}}\left( {\pmb{x}} \right)}/{{{g}_{{s}}}\left( {\pmb{x}} \right)} $, $ {{H}_{{r}}}\left( {\pmb{x}} \right) = {{{F}_{{r}}}\left( {\pmb{x}} \right)}/{{{g}_{{r}}}\left( {\pmb{x}} \right)} $.首先, 由$ {{g}_{{s}}}\left( {\pmb{x}} \right) $、$ {{g}_{{r}}}\left( {\pmb{x}} \right) $各产生P个相互独立的$ \upsilon $维随机数$ {{{\pmb{T}}}_{{s}1}}, \cdots , {{{\pmb{T}}}_{{s}P}} $和$ {{{\pmb{T}}}_{{r}1}}, \cdots , {{{\pmb{T}}}_{{r}P}} $, 并计算

采用逆变换方法由同一个$ \upsilon $维联合均匀分布$ U\left( 0, 1 \right) $产生$ {{{\pmb{T}}}_{{s}1}}, \cdots , {{{\pmb{T}}}_{{s}P}} $和$ {{{\pmb{T}}}_{{r}1}}, \cdots , {{{\pmb{T}}}_{{r}P}} $, 能够保证两组随机数具有较高的正相关程度, 进而使$ {Var}( {\hat{D}} ) $较小, 对$ \hat{D} $的估计值趋于稳定.

需要说明的是, $ D\left( {{F}_{{s}}}, {{F}_{{r}}} \right)\in \left[ 0, \infty \right) $仅是仿真和参考输出特征的联合CDF的差异(如图 1所示), 其取值范围为$ \left[ 0, \infty \right) $, 此时无法给出仿真和参考输出的一致性程度(即取值为$ \left[ 0, 1 \right] $的相对值).因此, 提出将$ D\left( {{F}_{{s}}}, {{F}_{{r}}} \right) $向可信度$ C\left( {{F}_{{s}}}, {{F}_{{r}}} \right) $转化的公式如下.

图 1  参考与仿真输出的CDF对比

Fig. 1  Comparing CDF curves of reference and simulation output

下载: 全尺寸图片 幻灯片

式中, $ \prod\nolimits_{i = 1}^{\upsilon }{\left( X_{i}^{\max }-X_{i}^{\min } \right)} $表示$ \upsilon $维样本空间所占区域的大小; $ X_{i}^{\min } = \min \left( X_{{s}i}^{\min }, X_{{r}i}^{\min } \right) $, $ X_{i}^{\max } = \max \left( X_{{s}i}^{\max }, X_{{r}i}^{\max } \right) $表示第i维变量的样本极值.显然, $ C\left( {{F}_{{s}}}, {{F}_{{r}}} \right) $满足如下性质^[13], 进而能够用于度量仿真模型可信度.

性质1. 非负性: $ C\left( {{F}_{{s}}}, {{F}_{{r}}} \right)\ge 0 $;

性质2. 交换性: $ C\left( {{F}_{{s}}}, {{F}_{{r}}} \right) = C\left( {{F}_{{r}}}, {{F}_{{s}}} \right) $;

性质3. 有界性: $ C\left( {{F}_{{s}}}, {{F}_{{r}}} \right)\in \left[ 0, 1 \right] $;

性质4. 同一性: $ C\left( {{F}_{{s}}}, {{F}_{{r}}} \right) = 1 $, 当且仅当$ {{F}_{{s}}} = {{F}_{{r}}} $.

3.3   多元输出仿真结果验证方法流程

基于前文所述方法, 给出考虑相关性的多元输出仿真结果验证流程如图 2所示.

图 2  考虑相关性的多元输出仿真结果验证方法流程

Fig. 2  Procedures of multivariate simulation result validation under correlation

下载: 全尺寸图片 幻灯片

1) 考虑不确定因素的影响, 分别进行n次仿真运行和实际试验, 获取多元仿真和参考输出$ {{{\pmb{Y}}}\!_{{s}}} = \left\{ {{{\pmb{Y}}}\!_{{s}1}}, {{{\pmb{Y}}}\!_{{s}2}}, \cdots, {{\bf{Y}}_{{s}m}} \right\} $, $ {{{\pmb{Y}}}\!_{{r}}} = \left\{ {{{\pmb{Y}}}\!_{{r}1}}, {{{\pmb{Y}}}\!_{{r}2}}, \cdots, {{{\pmb{Y}}}\!_{{r}m}} \right\} $;

2) 利用多元输出变量选择方法提取$ {{{\pmb{Y}}}\!_{{s}}} $、$ {{{\pmb{Y}}}\!_{{r}}} $的相关变量子集$ {{{\pmb{G}}}_{{s}i}}, i = 1, \cdots, {{\beta }_{{s}}} $, $ {{{\pmb{G}}}_{{r}j}}, j = 1, \cdots , {{\beta }_{{r}}} $;

3) 若$ {{\beta }_{{s}}} = {{\beta }_{{r}}} $且$ {{{\pmb{G}}}_{{s}i}} = {{{\pmb{G}}}_{{r}j}} $, 则依据式(8)$ \sim $(13)提取$ {{{\pmb{G}}}_{{s}i}} $、$ {{{\pmb{G}}}_{{r}j}} $中各变量的数据特征$ {\pmb{e}}_{{s}ik}^{l} $、$ {\pmb{e}}_{{r}ik}^{l} $; 反之, 若$ {{\beta }_{{s}}}\ne {{\beta }_{{r}}} $或$ {{{\pmb{G}}}_{{s}i}}\ne {{{\pmb{G}}}_{{r}j}} $的相关变量子集, 则认为该仿真模型不可信, 即C = 0, 算法结束;

4) 依据式(14)分别计算数据特征变量集$ e_{{s}i1}^{l}, e_{{s}i2}^{l}, \cdots , e_{{s}i{{m}_{{s}i}}}^{l} $和$ e_{{r}j1}^{l}, e_{{r}j2}^{l}, \cdots , e_{{r}j{{m}_{{r}j}}}^{l} $的联合CDF: $ {{F}_{{s}il}} $、$ {{F}_{{r}jl}} $;

5) 依据式(15)$ \sim $(19)计算特征变量集的联合CDF: $ {{F}_{{s}il}} $、$ {{F}_{{r}jl}} $的差异$ D_{i}^{l}\left( {{F}_{{s}il}}, {{F}_{{r}jl}} \right) $;

6) 依据式(20)将$ D_{i}^{l}\left( {{F}_{{s}il}}, {{F}_{{r}jl}} \right) $转化为可信度结果$ C_{i}^{l}\left( {{F}_{{s}il}}, {{F}_{{r}jl}} \right) $;

7) 通过2)可知, $ \alpha $个相关变量子集之间是相互独立的, 且用户关注的多个数据特征(包括位置、形状、频谱)间也可认为是独立的, 进而可采用加权方法综合多个可信度结果$ C_{i}^{l}\left( {{F}_{{s}il}}, {{F}_{{r}jl}} \right) $, $ l = 1, \cdots , {{L}_{i}} $; $ i = 1, \cdots , {{\beta }_{{s}}} $; $ j = 1, \cdots , {{\beta }_{{r}}} $. 图 2中“Integrate($ \cdot $)”表示加权综合算子.同时第$ \sigma $个动态输出的均值曲线$ {{\bar{y}}_{{s}\sigma }} $、$ {{\bar{y}}_{{r}\sigma }} $可认为是对系统输出的一次抽样, 不考虑不确定性影响时的多元输出数据是近似独立的, 进而综合得到最终验证结果如下所示.

其中, $ \sigma = 1, \cdots , {{m}_{{dynamic}}} $表示第$ \sigma $个动态输出变量, $ \zeta = 1, \cdots , {{L}_{\sigma }} $表示动态输出均值曲线的第$ \zeta $个特征, $ {{w}_{l}} $、$ {{w}_{\zeta }} $代表第l、$ \zeta $个数据特征的可信度结果权重, $ {{w}_{i}} $代表第i个相关变量子集的一致性分析结果权重, $ w_{1}^{-} $、$ w_{2}^{-} $代表相关变量子集和动态输出均值曲线一致性的权重.

4.   应用实例与分析

为验证本文方法的有效性, 针对文献[2]中给出的某飞行器纵向平面内末制导阶段的仿真模型进行结果验证.该模型包括飞行器制导模型和目标运动模型. 图 3给出纵向平面内弹目相对运动几何关系.目标以恒定速度$ {{v}_{{T}}} $沿$ x $轴向右行驶.假设飞行器无动力飞行且航向已对准目标, 忽略地球自转, 给出以时间为自变量的飞行器纵向质心运动方程

图 3  纵向平面内弹目相对运动几何关系

Fig. 3  Geometrical relationship of relative missile-target movement in longitudinal plane

下载: 全尺寸图片 幻灯片

式中, v为速度, $ \theta $为弹道倾角, h为高度, d为水平距离.阻力$ D = 0.5\rho {{v}^{2}}S{{C}_{{D}}}\left( Ma, \alpha \right) $, 升力$ L = 0.5\rho {{v}^{2}}S{{C}_{{L}}}\left( Ma, \alpha \right) $. $ {{C}_{{D}}} $与$ {{C}_{{L}}} $分别为阻力系数与升力系数. $ \alpha $为攻角, 马赫数$ Ma = v/{{v}_{{s}}} $, S为参考面积, M为质量, $ {{\alpha }_{{M}}} $为法向加速度, $ \lambda $为视线角.声速$ {{v}_{s}} $与大气密度$ \rho $根据1976年美国标准大气计算.相应的制导律设计可见文献[30].根据上述信息建立该飞行器纵向平面内末制导仿真模型.

利用此仿真模型精确研究此飞行器在末制导阶段的特性, 需要考虑其受到的不确定性因素.飞行器升力与阻力均存在不确定性, 引入升力系数扰动$ {{C}_{{LC}}} $, 与阻力系数扰动$ {{C}_{{DC}}} $模拟阻力系数与升力系数的不确定性, 因此有升力$ L = 0.5\rho {{v}^{2}}S{{C}_{{LC}}}{{C}_{{L}}}\left( Ma, \alpha \right) $, 阻力$ D = 0.5\rho {{v}^{2}}S{{C}_{{DC}}}{{C}_{{D}}}\left( Ma, \alpha \right) $, 分别采用不同分布对$ {{C}_{{LC}}} $和$ {{C}_{{DC}}} $进行描述.同时, 大气密度会影响升力和阻力, 因每次飞行环境不同需考虑其不确定性的影响, 采用大气密度系数$ {{C}_{{ }\!\!\rho\!\!{ }}} $表示.此外, 飞行器进入末制导阶段时的初始视线角$ {{\lambda }_{0}} $与弹道倾角$ {{\theta }_{0}} $亦具有不确定性.选取仿真模型和参考系统的不确定参数如表 2所示.

表 2  飞行器末制导过程的不确定参数取值

Table 2  Uncertainty parameters values in the terminal guidance process of flight vehicle

变量名	仿真模型参数分布	参考系统参数分布
大气密度系数${{C}_{{ }\!\!\rho\!\!{ }}}$	$N\left( 0, 0.033 \right)$	$N\left( 0, 0.033 \right)$
升力系数${{C}_{{D}}}$	$N(0, 0.05)$	$N(0.02, 0.07)$
阻力系数${{C}_{{L}}}$	$N(0, 0.033)$	$N(0.02, 0.033)$
初始弹道倾角${{{\theta }_{0}}}~/{{\rm rad}}$	$N\left( 0.17, 0.09 \right)$	$N\left( 0.26, 0.07 \right)$
初始视线角${{{\lambda }_{0}}}~/{{\rm rad}}$	$N\left( 0.17, 0.09 \right)$	$N\left( 0.17, 0.09 \right)$

下载: 导出CSV 

| 显示表格

选取用户关注的多元输出变量如表 3所示.选取静态输出变量有飞行器的最终落点位置坐标$ \left( {{x}_{{f}}}, {{z}_{{f}}} \right) $和目标终点位置坐标$ \left( {{x}_{{Tf}}}, {{z}_{{Tf}}} \right) $, 同时选取待验证的动态输出变量有弹道倾角$ \theta $、攻角$ \alpha $、视线角$ \lambda $、弹目相对距离$ {{D}_{{MT}}} $、目标速度$ {{v}_{{T}}} $.利用拉丁超立方抽样法, 对模型不确定性参数进行抽样, 给定初始样本数为1 000, 运行仿真模型共得到1 000组输出.改变飞行器模型参数(见表 2), 采用拉丁超立方抽样获得的1 000组数据作为参考输出.系统输出数据的包络线如图 4$ \sim $10所示.其中目标速度$ {{v}_{{T}}} $、目标终点位置$ {{z}_{{Tf}}} $为恒定值, 未在图中标出.

表 3  待验证的模型输出

Table 3  Model outputs to be validated

变量类型	变量名
动态	弹道倾角${\theta }~$(rad)
动态	攻角${\alpha }$ (rad)
动态	视线角${\lambda }~$(rad)
动态	弹目相对距离${{{D}_{{MT}}}}~$(m)
动态	目标速度${{v}_{{T}}}~$(m/s)
静态	飞行器落点X坐标${{{x}_{{f}}}}~$(m)
静态	飞行器落点Z坐标${{{z}_{{f}}}}~$(m)
静态	目标终点位置X坐标${{{x}_{{Tf}}}}~$(m)
静态	目标终点位置Z坐标${{{z}_{{Tf}}}}~$(m)

下载: 导出CSV 

| 显示表格

图 4  弹道倾角输出包络线

Fig. 4  Envelope lines of launch angle

下载: 全尺寸图片 幻灯片

图 5  攻角输出包络线

Fig. 5  Envelope lines of angle of attack

下载: 全尺寸图片 幻灯片

图 6  视线角输出包络线

Fig. 6  Line-of-sight angle envelopes

下载: 全尺寸图片 幻灯片

图 7  弹目相对距离输出包络线

Fig. 7  Envelope lines of the missile-target relative distance

下载: 全尺寸图片 幻灯片

图 8  飞行器落点X坐标输出散点图

Fig. 8  Scatter diagram of X-direction drop point coordinates of the flight vehicle

下载: 全尺寸图片 幻灯片

图 9  目标终点位置X坐标输出散点图

Fig. 9  Scatter diagram of X-direction terminal point coordinates of the target vehicle

下载: 全尺寸图片 幻灯片

图 10  飞行器落点Z坐标输出散点图

Fig. 10  Scatter diagram of the terminal point of the target in the Z direction

下载: 全尺寸图片 幻灯片

利用本文方法对该飞行器末制导仿真输出进行验证.首先利用多元输出变量选择方法分别对仿真和参考的动静态输出变量进行相关性分析及变量选择, 得到相关性分析结果如表 4所示.通过分析可知, 仿真和参考输出变量具有相同的变量子集, 动态输出变量$ \theta $、$ \alpha $、$ \lambda $、$ {{D}_{{MT}}} $具有相关性, 故将其归为一类. $ {{v}_{{T}}} $为定值(不随时间改变)并与变量子集Ⅰ相互独立; 静态输出$ {{x}_{{f}}} $、$ {{x}_{{Tf}}} $具有相关性, 通过验证可知两者满足线性关系(如图 11所示), 同时$ {{z}_{{f}}} $与$ {{x}_{{f}}} $相互独立(如图 12所示), 进而可得$ {{z}_{{f}}} $与$ {{x}_{{f}}} $相互独立, $ {{z}_{{Tf}}} $为定值0形成了变量子集Ⅲ.由上述分析结果验证了变量选择方法的有效性.

表 4  多元输出变量选择结果

Table 4  Variables selection results of multiple outputs

输出类型	变量子集Ⅰ	变量子集Ⅱ	变量子集Ⅲ
动态	$\theta $, $\alpha$, $\lambda$, ${{D}_{{MT}}}$	${{v}_{{T}}}$	-
静态	${{x}_{{f}}}$, ${{x}_{{Tf}}}$	${{z}_{{f}}}$	${{z}_{{Tf}}}$

下载: 导出CSV 

| 显示表格

图 11  飞行器落点X坐标与目标终点位置X坐标间的关系

Fig. 11  Relationship of X-direction coordinates between drop point of flight vehicle and terminal point of target

下载: 全尺寸图片 幻灯片

图 12  飞行器落点X方向坐标与Z方向坐标间的关系

Fig. 12  Relationship between X-direction and Z-direction coordinates of the drop point of flight vehicle

下载: 全尺寸图片 幻灯片

根据表 4的变量选择结果求取各变量子集关于某特征的联合CDF, 选取动态输出的位置和形状特征, 分别求取变量子集Ⅰ的联合CDF, 变量子集Ⅱ为恒定值在验证过程中直接采用相对误差方法进行一致性分析即可; 对于静态输出变量子集Ⅰ关于数据本身的联合CDF如图 13所示, 变量子集Ⅱ的CDF曲线如图 14所示.进而得到动态输出均值曲线的一致性结果(见表 5)以及多个变量组关于多个特征的CDF差异和可信度结果(见表 6).依据式(21)综合多个可信度结果得到最终验证结果为0.82, 由于仿真和参考输出变量$ {{v}_{{T}}} $、$ {{z}_{{Tf}}} $均相等且恒为0, 故在计算模型可信度时不予考虑, 为方便计算采用均权的方式进行综合.

图 13  仿真和参考静态输出变量子集Ⅰ的联合CDF对比

Fig. 13  JCDF comparison of variable subset I between static simulation and reference output

下载: 全尺寸图片 幻灯片

图 14  仿真和参考静态输出变量子集Ⅱ的CDF对比

Fig. 14  Comparison of variable subset Ⅱ between static simulation and reference output

下载: 全尺寸图片 幻灯片

表 5  动态输出均值曲线的一致性分析结果

Table 5  Consistency analysis results of the mean curves of dynamic outputs

变量名	位置特征一致性	形状特征一致性
$\theta $	0.92	0.74
$\alpha$	0.63	0.60
$\lambda$	0.98	0.74
${{D}_{{MT}}}$	0.97	0.61

下载: 导出CSV 

| 显示表格

表 6  仿真和参考输出变量子集的一致性分析结果

Table 6  Consistency analysis results of the variables subset of the simulation and reference outputs

输出变量类型	变量子集标号	累积概率分布差异	可信度结果
动态	变量子集Ⅰ	位置差异: $8.92\times {{10}^{{-8}}}$	位置特征: 0.99
动态	变量子集Ⅰ	形状差异: $1.1\times {{10}^{{-3}}}$	形状特征: 0.94
动态	变量子集Ⅱ	0	1
静态	变量子集Ⅰ	$1.6\times {{10}^{5}}$	0.84
静态	变量子集Ⅱ	0.5	0.9
静态	变量子集Ⅲ	0	1

下载: 导出CSV 

| 显示表格

此外, 为进一步验证本文方法对参数不确定性度量的有效性, 针对上述应用实例分别设计两组验证实验(不确定性参数取值见表 7).固定仿真模型和参考系统的不确定性参数大气密度系数$ C_ \rho $、升力系数$ {{C}_{{D}}} $、阻力系数$ {{C}_{{L}}} $和初始视线角$ {{\lambda }_{0}} $的取值.分别调节仿真模型初始弹道倾角$ {{\theta }_{0}} $的均值和方差, 得到最终验证结果如图 15$ \sim $16所示.通过实验可得, 该方法能够度量仿真模型不确定参数取值的离散程度对验证结果的影响, 证明过大或过小的参数不确定度均会降低模型的可信度; 同时该方法能够度量不确定性参数的均值差异对验证结果的影响.综上所述, 所提方法能够用于解决带有相关性的多元输出仿真结果验证问题.

表 7  验证实验的不确定参数取值

Table 7  Uncertainty parameters values for validation experiments

试验编号	参考系统${{\theta }_{0}}$取值	实验组Ⅰ ${{\theta }_{0}}$取值	实验组Ⅱ ${{\theta }_{0}}$取值
1	$N\left( 0.26, 0.07 \right)$	0.26	$N\left( 0.15, 0.07 \right)$
2	$N\left( 0.26, 0.07 \right)$	$N\left( 0.26, 0.04 \right)$	$N\left( 0.21, 0.07 \right)$
3	$N\left( 0.26, 0.07 \right)$	$N\left( 0.26, 0.07 \right)$	$N\left( 0.26, 0.07 \right)$
4	$N\left( 0.26, 0.07 \right)$	$N\left( 0.26, 0.1 \right)$	$N\left( 0.31, 0.07 \right)$
5	$N\left( 0.26, 0.07 \right)$	$N\left( 0.26, 0.13 \right)$	$N\left( 0.37, 0.07 \right)$

下载: 导出CSV 

| 显示表格

图 15  实验组Ⅰ验证结果

Fig. 15  Validation result of experiment Ⅰ

下载: 全尺寸图片 幻灯片

图 16  实验组Ⅱ验证结果

Fig. 16  Validation result of experiment Ⅱ

下载: 全尺寸图片 幻灯片

5.   结论

针对带有相关性的多元输出仿真模型验证问题, 提出了考虑不确定性的联合验证方法.首先对多变量输出提取相关变量子集, 并对各输出变量提取数据特征, 利用联合CDF差异法度量各相关变量子集的一致性程度, 进而综合得到模型可信度.利用单变量验证方法进行多变量验证时需要满足输出变量相互独立的条件.本文方法考虑了多变量间的相关关系, 基于相关变量子集进行联合验证, 较单变量验证方法应用更合理; 同时在验证前引入了变量相关性分析, 使其能够适应输出变量之间关系未知的情况, 使验证结果更准确, 但也增加了计算开销.此外, 该方法能够度量不确定性因素对模型可信度的影响.

需要说明的是, 本文仅考虑同一类型输出(动态或静态)存在相关性的情况, 涉及的变量选择方法本质上属于数据挖掘方法, 为确保方法的准确性要求具备足够的样本容量, 对于参考数据缺乏的情况, 可采用专家给出参考数据的大致分布, 或可利用已有的历史数据、可信度较高且类似的半实物/纯数字仿真模型所产生的数据代替.此外, 刻画动态输出的数据特征不限于距离、形状及频谱, 可依据具体应用需求而定(例如, 超调量、相位误差等).后续将对动态、静态输出间的相关性分析及变量选择方法进行研究; 同时针对参考数据缺乏以及存在认知和固有混合不确定性时的多元输出仿真结果验证问题展开进一步研究.