-
开放科学（资源服务）标识码（OSID）：
-
A/B基准值是表征航空复合材料层合板、高强铝合金、钛合金等关键材料性能统计分布下限的核心指标，是飞机机翼、机身等主承力结构强度设计与适航认证的主要依据。其中，A基准值表示95%的置信度下99%的性能数值群的最小值；B基准值表示95%的置信度下90%的性能数值群的最小值。在金属材料设计中，根据结构的具体需求选用相应的基准值。在需要冗余设计的情况下，倾向于采用B基准值[1-4]。例如地梁，即使单个构件发生故障，也能安全地将负荷重新分配到其他构件，且不会超过结构的最大承载极限。而对于依赖单一载荷路径的结构，例如机翼和机身，通常建议使用更严格的A基准值，因为这些结构的故障可能导致飞机遭受灾难性的后果。因此，与A基准值相比，B基准值提供了更大的设计容许范围。研究A/B基准值的统计方法为结构强度设计提供了可靠依据，极大提升了飞机在复杂载荷环境下的安全裕度。同时，精确的许用值评估有助于充分挖掘材料性能潜力，实现结构轻量化，显著改善燃油经济性与运载效益。此外，在面对高昂试验成本与小样本挑战时，健全的统计方法能够优化试验设计，支撑适航符合性。
众多学者对B基准值的确定方法进行了深入探讨。B基准值的确定方法主要有蒙特卡罗模拟法与统计分析法，这2种方法都基于材料和几何特性及其统计分布，但在确定B基准值的策略上有所区别。Cardoso等[5]基于仿真数据训练人工神经网络代理模型，并借助该模型求解典型载荷工况下的A/B基准值。Yan等[6]提出一种有限拉伸试验数据结合数字孪生模型分析的混合方法，可以高效获取SiC/SiC复合材料拉伸的B基准值。Furtado等[7]将机器学习技术应用于预测复合材料层压板的统计设计容许值，通过分析生成的数据，使用4种机器学习算法预测复合材料层压板的缺口强度及其统计分布，以及相关材料属性和几何特征的不确定性。结果表明除了随机森林算法外，所有机器学习算法均表现出色，特别是高斯过程算法在数据点较少时表现突出，而人工神经网络在训练集较大时更为有效。本文将重点研究统计分析法，对于蒙特卡罗模拟法不作过多展开。
国内学者对基于统计分析法的A/B基准值估计进行了大量研究[8-10]。傅惠民[11]提出了一种新的单侧容限系数法，这种方法特别适用于正态分布和对数正态分布，并且在样本量较小的情况下，相比国际上常用的单侧容限系数法更精确地近似真实值。此外，傅惠民[12]还提出二维单侧容限系数法，该方法可以利用过往数据，在维持同等精度的同时，减少一半的试验量。在三参数威布尔分布领域，傅惠民等[13]提出了一种基于秩的非参数化方法，并开发了多种回归分析工具。马小兵等[14]为小样本提供了一种新的Fiducial推断方法，用于确定分位数的置信区间。赵新攀等[15]采用Bootstrap技术，极大提高了小样本概率加权矩法的准确性。章洛[16]基于二分迭代算法给出了三参数威布尔分布参数估计A/B基准值的解析式，开发了融合相关系数法、概率加权矩估计与Bootstrap重采样技术的通用计算程序，可处理任意设定的置信度水平与性能数群规模的最小值计算问题。
除了以上经典的统计方法，贝叶斯理论在可靠性分析领域也备受关注。为了获得估计参数的后验分布，需要将历史数据等多种信息源有效地纳入现场数据进行推理。Hannig等[17]提出了广义先验推断(GFI)方法，由故障数据提供其先验信息，为贝叶斯方法提供了一种客观的替代方法。研究表明，广义基准推断区间在保持所述覆盖率的同时，其平均区间长度比其他方法更短。Yan等[18]使用GFI方法分别估计了广义指数分布和洛马克斯分布的参数区间。结果表明，GFI方法的均方根误差小于传统的参数估计方法。Yang等[19]根据GFI方法对三参数威布尔分布进行点估计和区间估计，与无先验信息的贝叶斯估计相比更加稳健、准确。总体来说，基于统计分析法计算三参数威布尔分布A/B基准值的方法较少，尤其是在小样本下的估计结果均不太理想。
本文针对服从三参数威布尔分布的非结构数据(单批样本数据)，构建了4种母体置信限曲线方法(非参数估计与相关系数法结合的NP-CCM方法、非参数估计与最小偏差法结合的NP-MDM方法、非参数估计与误差变量法结合的NP-EIV方法和非参数估计与最小二乘法结合的NP-LSM方法)进行A/B基准值的统计分析。在不同的样本量和不同的威布尔分布下，将4种母体置信限曲线方法与MMPDS-16手册中的方法[20]进行对比分析，以验证本文方法的适用性与准确性。
HTML
-
根据MIL-HDBK-17F标准，在计算基准值时首先推荐使用三参数威布尔分布，其次是正态分布和对数正态分布。尽管美国联邦航空管理局(FAA)指出使用威布尔分布可能会得出较为保守的结果，但众多文献均表明，工程材料的性能通常遵循三参数威布尔分布。因此，本文优先采用三参数威布尔分布进行数据拟合分析。
三参数威布尔分布W(β，η，γ)的概率密度函数和累积分布函数分别为：
其中：β、η和γ分别为形状参数、尺度参数和位置参数。
-
异常值通常指在一组数据中显著低于或高于其他观测值的数值。在进行统计分析之前，需要对数据集进行异常值的检测。本文使用最大范数残差法识别数据集中可能存在的异常值。
采用最大范数残差法，假设除异常值之外，样本数据遵循正态分布。如果一个样本观测值与样本均值的绝对偏差显著超出样本标准差的阈值，则该观测值将被视为异常值。最大范数残差法每次仅能识别出一个异常值，这使其更适合对单个数据点进行评估。
最大范数残差统计量的定义为：
式中：xi为样本观测值，i=1，2，…，n；x为样本均值；s为样本标准差。
样本标准差为：
最大范数的临界值C定义为：
式中：t为自由度n-2的t分布的
$1-\frac{\alpha}{2 n}$ 分位数；α为显著性水平，本文取0.05。当最大范数残差统计量MNR低于临界值C时，以1-α的置信水平认为样本中没有异常值；如果MNR超过临界值C，则以相同的置信水平认为与MNR相关的观测值可能是异常值。
-
MMPDS-16手册中针对三参数威布尔分布提出改进的Anderson-Darling检验方法，用于确定给定数据集的曲线是否可以使用三参数威布尔曲线来拟合。该检验的实质是对观测数据的累积分布函数与整个测量范围内拟合威布尔曲线的累积分布函数进行数值比较，该检验与原始的Anderson-Darling检验不同，它强调低尾部。该检验方法可用于完全数据或删减数据的分析处理。
MMPDS-16手册介绍了一种估计三参数威布尔分布位置参数的方法，可用于拟合优度检验和计算T99或T90值(A/B基准值)。令K小于等于
$\min \left\{\frac{4 n}{15}, (1-p) \frac{n}{3}\right\}$ 的最大整数，p为右尾部被删减的比例(p=0、0.2、0.5)，定义函数W(γ)为：式中：
使用MMPDS-16手册中的表 9.10.8确定参数r。表 9.10.8中的第1列用于估计与此处所述的Anderson-Darling拟合优度检验相关的位置参数γ50，第2列和第3列用于估算γ99和γ90，可确定T99和T90。方程W(γ)=r的解是位置参数γ的估计值。
当得到位置参数估计值后，将样本值减去位置参数估计值得到新的样本值，此时相当于两参数威布尔分布，再使用极大似然估计对形状参数和尺度参数进行估计，得到3个参数的估计值。估计出参数后，可根据以下步骤计算Anderson-Darling统计量。
当i=1，…，r时，则：
令fn+1=1，则：
Anderson-Darling统计量定义为：
如果存在：
可以得出抽取样本不是三参数威布尔分布的结论(误差风险为5%)。否则，不能拒绝样本服从三参数威布尔分布的假设。
-
为了计算三参数威布尔分布总体的置信下限值，必须具备以下条件：①位置参数的估计；②形状和尺度参数的估计；③三参数威布尔分布的单侧容限系数表。若x1，…，xn表示任意顺序的样本观测值，x(1)，…，x(n)表示从小到大顺序排列的样本观测值，使用1.2节中描述的方法分别计算T99和T90的位置参数估计值，并分别用γ99和γ90表示。此外，形状参数估计值用β99和β90表示，尺度参数估计值用η99和η90表示。
当样本观测值减去估计的位置参数，则形状参数和尺度参数可使用两参数威布尔分布的极大似然法进行估计。根据样本{x(i)-γ99：i=1，…，r}计算估计值β99和η99；同样，根据样本{x(i)-γ90：i=1，…，r}计算估计值β90和η90。
A/B基准值计算公式分别为：
其中：
$Q_{99}=\eta_{99}(0.010\;05)^{\frac{1}{\beta_{99}}}$ ；$Q_{90}=\eta_{90}(0.010\;05)^{\frac{1}{\beta_{90}}}$ 。V99和V90的通用计算公式为：
式中：
$g(p)=0.45+0.779\;7 \ln [-\ln (1-p)]$ ；计算V99时p=0.01，计算V90时p=0.1；d=0.779 696 8；c=1.645；$k_n=\sqrt{\frac{n}{n-1}}$ ；无截尾数据时，a00=0.607 9、a01=-0.474 0、a11=0.977 5。
1.1. 异常值的检测
1.2. 分布参数检验
1.3. 基准值的计算
-
相关系数法(Correlation Coefficient Method，CCM)是当样本的相关系数取得最大值时得到位置参数估计值，然后使用最小二乘法求得形状参数和尺度参数的方法。
样本的相关系数为：
式中：
$x_i=\ln \left(t_i-\gamma\right)$ ；$\bar{x}=\sum\limits_{i=1}^n \frac{x_i}{n}$ ；$y=\ln \left[-\ln \left(1-F_i\right)\right]$ ；$\bar{y}=\sum\limits_{i=1}^n \frac{y_i}{n}$ 。F为累积分布函数，可以使用中位秩来估计，即：
式中：i为按升序排列的寿命样本序号，i=1，2，…，n；n为样本量。
获得位置参数估计后，将原来的样本减去位置参数得到新的样本，使用最小二乘方程可得到形状参数和尺度参数，即：
-
最小偏差法(Minimum Discrepancy Method，MDM)由Xie等[21]提出，通过将尺度参数之间的差异最小化，可估计出正确的形状参数和位置参数。
由式(2)转换为：
式中：t为寿命数据。累积分布函数可采用式(21)估计。
形状参数和位置参数的计算公式为：
威布尔尺度参数的估计为：
-
误差变量法[22](Errors-in-Variables，EIV)是当样本观测值与根据累积分布函数转换得到的寿命平方差最小时，估计出形状、尺度和位置参数的方法，即：
式中：ti为失效数据。
tsi的表达式如下：
式中：累积分布函数可采用式(21)进行估计。
-
最小二乘法(Least Squares Method，LSM)是累积分布函数的偏差平方和最小时，估计出形状、尺度和位置参数的方法，即：
式中：累积分布函数可采用式(21)进行估计。
-
样本t1，t2，…，tn来自三参数威布尔分布累积分布函数F(t；β，η，γ)总体，样本次序统计量为t(1)，t(2)，…，t(n)。Pi为第i个秩统计量，表示总体中个体小于第i个观测值的百分率，Pi=F(t(1)；β，η，γ)，i=1，2，…，n，其概率密度函数为：
Pui和Pli定义如下：
Pui和Pli表示置信度为1-α的非参数置信上限和置信下限，即：
由于秩分布是一个贝塔分布，可将上式变换得到：
由式(34)求得置信度1-α的概率P=F(tP；β，η，γ)的总体百分位值tP的置信下限：
同样，总体百分位值tP的置信上限可由下式求得：
-
由于R=1-P，使用式(34)和式(35)可计算出可靠寿命tR的置信上下限。利用式(34)和式(35)可代替三参数威布尔分布相关系数法的中位秩式(21)估计出βl，ηl，γl和βu，ηu，γu。通过蒙特卡罗抽样从威布尔分布W(2，1 00 0，1 00 0)中随机产生20个样本，分别为1 197.10、1 305.94、1 374.27、1 375.18、1 449.16、1 472.29、1 504.98、1 738.40、1 773.36、1 916.89、1 957.51、2 154.36、2 160.85、2 305.60、2 388.28、2 389.72、2 412.30、2 496.05、2 572.18和3 318.25。当置信水平为90%时，使用不同的经验分布函数式(21)、式(34)和式(35)估计该样本的3个参数，估计结果如表 1所示。
根据表 1的参数估计结果得到以下3个累积分布函数：
将式(36)-式(38)的累积分布函数分别在图 1中表示。由图 1可以看出，深蓝色曲线既可以表示累积分布函数的置信上限，也可以表示可靠寿命的置信下限；浅蓝色曲线既可以表示累积分布函数的置信下限，也可以表示可靠寿命的置信上限；红色虚线表示点估计的累积分布。
本节采用4种母体置信限曲线方法(NP-CCM、NP-MDM、NP-EIV和NP-LSM)对A/B基准值进行估计。
NP-CCM：利用式(34)代替三参数威布尔分布相关系数法的中位秩式(21)估计βl、ηl、γl，根据可靠寿命公式得到非参数估计和相关系数法结合的A/B基准值评估方法，即：
NP-MDM：利用式(34)代替三参数威布尔分布最小偏差法的中位秩式(21)估计βl、ηl、γl，根据可靠寿命公式得到非参数估计和最小偏差法结合的A/B基准值评估方法，计算公式同式(39)和式(40)。
NP-EIV：利用式(34)代替三参数威布尔分布误差变量法的中位秩式(21)估计βl、ηl、γl，根据可靠寿命公式得到非参数估计和误差变量法结合的A/B基准值评估方法，计算公式同式(39)和式(40)。
NP-LSM：利用式(34)代替三参数威布尔分布最小二乘法的中位秩式(21)估计出βl、ηl、γl，根据可靠寿命公式得到非参数估计和最小二乘法结合的A/B基准值评估方法，计算公式同式(39)和式(40)。
2.1. 相关系数法
2.2. 最小偏差法
2.3. 误差变量法
2.4. 最小二乘法
2.5. 非参数估计
2.6. 4种母体置信限曲线法
-
为比较4种母体置信限曲线方法(NP-CCM、NP-MDM、NP-EIV和NP-LSM)与MMPDS-16标准方法的估计结果，本文进行了大量的蒙特卡罗仿真计算。对不同威布尔分布W(2，1 000，1 000)、W(2，1 000，2 000)、W(3，1 000，1 000)和W(3，1 000，2 000)进行蒙特卡罗抽样，样本量选用n=30、50、100和200，重复抽取1 000次，使用以上5种方法进行A/B基准值估计并进行对比。随机选取100组样本，采用5种方法估计W(2，1 000，1 000)和W(3，1 000，1 000)的A/B基准值(n=30)，并表示在图 2-图 3中；计算出1 000次估计结果的偏差(Std)和均方根误差(RMSE)，并表示在图 4-图 7中。
1) 针对A基准值，以上5种方法分别在样本量为n=30、50、100和200，威布尔分布为W(2，1 000，1 000)、W(2，1 000，2 000)、W(3，1 000，1 000)和W(3，1 000，2 000)的不同情况下进行比较。在样本量较小时，4种母体置信限曲线法估计的A基准值结果比MMPDS-16标准方法估计的结果更加稳健、更加准确。当样本量增加到100和200时，标准方法的计算结果更加冒进，因此不推荐MMPDS-16标准方法。
2) 针对B基准值，以上5种方法分别在样本量为n=30、50、100和200，威布尔分布为W(2，1 000，1 000)、W(2，1 000，2 000)、W(3，1 000，1 000)和W(3，1 000，2 000)的不同情况下进行比较。在样本量较小时，4种母体置信限曲线方法估计的B基准值结果比MMPDS-16标准方法估计的结果更加稳健，MMPDS-16标准方法估计B基准值的结果比4种母体置信限曲线法估计的结果更加准确，因此从准确性角度推荐MMPDS-16标准方法。
-
选用文献[20]的一组航空材料单批样本(数据见表 2)进行实例分析。首先根据异常值的检测，计算出最大范数残差统计量MNR低于临界值C(MNR=2.896，C=3.384，α=0.05)，以1-α的置信水平认为样本中没有异常值，并且通过了Anderson-Darling检验(AD＜0.395 1+4.186×10-5n，AD=0.196)。再分别采用4种母体置信限曲线方法(NP-CCM、NP-MDM、NP-EIV和NP-LSM)以及MMPDS-16标准方法进行A/B基准值的计算，计算结果如表 3所示。
表 3列出了以上5种方法计算实例样本A/B基准值的结果，构建的4种母体置信限曲线方法估计的A/B基准值都与MMPDS-16标准方法的计算结果十分接近，说明这4种母体置信限曲线方法是可行的。从上述100个样本中随机抽取30个样本，分别使用5种方法进行A/B基准值的计算。将样本量为100时MMPDS-16标准方法的计算结果作为真实值，同时将样本量为30时的5种方法计算结果进行对比，计算相对误差的绝对值列于表 4中。由表 4可知，在样本量为30时4种母体置信限曲线方法计算A/B基准值的结果比MMPDS-16标准方法的计算结果更加准确。
同样，从上述100个样本中随机抽取50个样本，分别使用5种方法进行A/B基准值的计算。将样本量为100时MMPDS-16标准方法的计算结果作为真实值，同时将样本量为50时的5种方法计算结果进行对比，计算相对误差的绝对值列于表 5中。由表 5可知，在样本量为50时4种母体置信限曲线方法计算A/B基准值的结果比MMPDS-16标准方法的计算结果更加准确。
-
1) 针对符合三参数威布尔分布的非结构数据进行统计分析，构建了4种母体置信限曲线方法(NP-CCM、NP-MDM、NP-EIV和NP-LSM)估计A/B基准值。
2) 在不同的样本量和不同的威布尔分布下，将4种母体置信限曲线方法与MMPDS-16标准方法进行对比，在样本量较小时4种母体置信限曲线方法估计的A基准值结果比MMPDS-16标准方法估计的结果更加稳健、准确；当样本量增加到100和200时，MMPDS-16标准方法的计算结果更加冒进，因此不推荐MMPDS-16标准方法。
3) 4种母体置信限曲线方法估计的B基准值结果比MMPDS-16标准方法估计的结果更加稳健，MMPDS-16标准方法估计B基准值的结果比4种母体置信限曲线方法估计的结果更加准确，从准确性角度推荐MMPDS-16标准方法。
4) 使用4种母体置信限曲线方法和MMPDS-16标准方法对一组航空材料单批样本进行了实例分析，验证了4种方法的实用性和准确性。本文构建的4种母体置信限曲线法为小样本条件下航空材料A/B基准值的高精度估计提供了有效且可靠的解决方案。
基于本文的研究基础与现有局限，未来可以结合贝叶斯理论引入先验信息，构建贝叶斯母体置信限，完成小样本下A/B基准值的评估。
DownLoad: