随着生产实践和科学技术的发展，压力传感器广泛应用于航空航天、武器装备、工业生产等领域^[1-4]。目前，对于压力传感器静态特性的研究已经比较深入和全面，然而，在实际应用中大部分被测压力的幅值随时间不断变化，呈现一定的动态特性，如果用静态特性指标来评价压力传感器在动态测量中的性能必然会带来误差。因此，有必要对压力传感器进行动态校准，获取其时域和频域动态特性参数在测量时的变化规律^[5]。

激波管作为一种能产生近似理想阶跃压力信号的装置，广泛应用于压力传感器的动态校准^[6-8]。通常，需要校准的压力传感器时域动态特性参数包括上升时间、调节时间和超调量；频域动态特性参数包括谐振频率和工作频带^[9]。在实际校准中，由于受到激波二次冲击、温度剧烈变化、外界噪声等影响，每次动态校准得到的压力传感器动态特性参数值不同，导致重复校准结果中存在着差异较大的不确定度。时频域动态特性参数的不确定度幅值直接影响压力传感器的测量范围和应用场合。因此，对压力传感器的动态特性参数进行合理的不确定度评定，不仅可以得出传感器的动态测量性能，还可以成为动态校准装置不确定度溯源的基础，具有十分重要的理论意义和应用前景。

近年来，对压力传感器动态校准不确定度的评定研究已经取得了一定的进展。Crowder和Moyer^[10]提出两步蒙特卡罗法评定非线性系统的测量不确定度，为压力传感器校准的不确定度提供了一种新思路。文献[11]提出从动态校准系统着手，分析影响校准结果的主要因素，建立各因素的不确定度传递模型，通过合成得出动态校准总不确定度。该方法通过建立不确定度传递链，为后续不确定度的溯源提供基础，但是由于不确定度评定结果的可靠性过于依赖各因素，在实际校准过程中，不同校准实验的不确定度来源往往存在着很大的差异，许多因素的不确定度传递模型又很难获得，导致计算过程的实施非常困难，通用性比较差。文献[12]采用正弦压力发生装置对压力传感器进行动态校准，在假设各校准影响因素互不相关的前提下，通过合成得到幅值灵敏度和相移不确定度。然而该方法对动态校准过程中不确定度来源的分析不够全面，得到结果的准确度存在着很大的质疑。文献[13]提出了压力传感器动态校准不确定度的灰色评定模型，采用正弦压力对传感器进行重复校准，根据传感器的输出信号，采用灰色方法计算特定频率点的不确定度，通过最小二乘法拟合得到工作频带中所有频率点的不确定度。在实际动态校准中，由于噪声的干扰，传感器的输出信号往往不能精确表征传感器的特性，因此不确定度评定精度受校准装置噪声的影响很大。文献[14]提出一种模型辨识的压力传感器动态校准不确定度评定方法，采用广义最小二乘法对压力传感器进行建模，根据模型估计传感器工作频带值，采用贝塞尔方法计算工作频带的动态不确定度。该方法通过传感器的数学模型获取工作频带值，得到的不确定度受噪声干扰比较小，但是直接采用传感器的响应数据进行建模，使模型的精度难以得到有效的提高，加之没有考虑传感器的其他时频域动态特性参数，且在重复校准数据量较少的情况下使用统计方法计算不确定度，降低了不确定度评定结果的可靠性。

现有方法虽然在一定程度上初步实现了压力传感器动态校准不确定度的评定，但是仍然存在许多问题。例如传感器数学模型的准确性普遍比较差、动态特性参数的分析过于单一、动态校准不确定度的计算方法不合理等。本文提出一种压力传感器动态特性参数不确定度评定方法，全面分析了校准过程中噪声对校准结果的影响，充分考虑了传感器的时频域动态特性参数。首先采用基于经验模态分解(Empirical Mode Decomposition, EMD)的传感器输出信号预处理方法消除噪声对压力传感器动态校准结果的影响；然后采用自适应最小二乘法估计模型的最优阶数和参数，提高传感器数学模型的可靠性；最后针对重复校准得到的动态特性参数序列样本量小的特点，提出一种基于自助采样的不确定度评定方法，得到了被校准压力传感器动态特性参数在不同置信概率下的扩展不确定度和相对不确定度。

1 动态特性参数不确定度评定方法 1.1 评定原理

压力传感器动态特性参数的不确定度评定过程包括以下5个步骤：

1) 利用激波管系统对给定的压力传感器进行动态校准，得到传感器的输出信号。

2) 采用经验模态分解对传感器输出信号进行预处理，减小噪声对输出信号的影响。

3) 根据传感器输入输出信号，采用自适应最小二乘法建立传感器的数学模型，得到传递函数。

4) 根据传递函数模型获取压力传感器的阶跃响应曲线和幅频曲线，计算时频域动态特性参数。

5) 多次重复动态校准实验，得到压力传感器动态特性参数序列，采用自助法计算参数自助样本的概率密度直方图，估计不同置信水平下的扩展不确定度和相对不确定度。

评定原理如图 1所示。

1.2 信号预处理

在采用激波管对压力传感器动态校准过程中，由于噪声的干扰，被校准传感器的输出通常出现失真现象，如果直接进行建模就会得到不可靠的模型，影响动态特性参数的估计精度，因此需要在建模之前对原始输出信号进行预处理。由于传感器的理想输出未知，采用传统滤波器难以估计出其截止频率。本文采用经验模态分解方法将原始输出信号分解为不同频带的子序列；根据各子序列的频谱与传感器输出信号频谱之间的关系，剔除其中的噪声子序列，减小传感器输出信号的失真。其中分解得到的子序列称为本征模态函数(Intrinsic Mode Function, IMF)，满足2个条件^[15]：

1) 整个数据长度中极值点和过零点的个数必须相等或者至多相差一个。

2) 任意点处的极大值包络线和极小值包络线的均值为零。

假设传感器的原始输出信号为y₀(t)，具体分解步骤如下：

步骤1  提取y₀(t)的全部局部最大值和局部最小值。采用三次样条曲线分别连接所有局部最大值和局部最小值，得到的曲线分别称为上包络线u(t)和下包络线l(t)。

步骤2  计算上包络线和下包络线的均值为

步骤3  y₀(t)与m₁(t)之间的差值为

如果h₁⁽¹⁾(t)满足IMF的2个条件，则其为y₀(t)的第1个IMF分量；否则令y₀(t)=h₁⁽¹⁾(t)，将步骤1~步骤3重复k次，直到h₁^(k)(t)满足IMF的2个条件^[16]为止。此时的第1个IMF分量记为

显然，c₁(t)包含y₀(t)中的最高频成分。

步骤4  从原始信号y₀(t)中减去c₁(t)，得到对应的残余分量为

将r₁(t)当作y₀(t)重复以上步骤i次，则第i个IMF分量c_i(t)被提取出来，表示为

继续执行上述分解步骤，直到最终的残余分量r_m(t)成为单调函数或者只有一个极值点，此时不能再从中提取出更多的IMF分量。

将式(4)与式(5)相加，则传感器原始输出信号可以表示为多个IMF分量和最终残余分量之和

被分解出来的IMF分量的频带从高到低，如果某些IMF分量的频带远离原始信号y₀(t)的振铃频率，则认为这些IMF为噪声分量，予与剔除。将剩下的IMF分量之和进行重构，得到预处理后的传感器输出信号y(t)。

1.3 压力传感器模型辨识

得到预处理的压力传感器输出信号之后，采用自适应最小二乘法估计传感器数学模型。压力传感器一般可以描述为单输入、单输出的时不变线性系统^[17]，其差分方程为

式中：x(k)、y(k)(k=1, 2, …, N)为传感器的输入、输出序列，N为序列长度；ε(k)为随机噪声；a_i、b_i(i=1, 2, …, n)为模型参数，n为模型阶数。

激波管产生的阶跃压力信号通常可以认为是理想阶跃信号^[18]，因此输入信号中的随机噪声可以忽略，则式(7)中随机噪声为输出序列中的噪声，记为ε(k)=e_y(k)。基于自适应最小二乘的模型参数辨识方法的基本步骤如下：

步骤1  根据压力传感器的输入序列x(k)和输出序列y(k)，将式(7)表示为

式中：

当k＜0时，x(k)=y(k)=0。

采用最小二乘法进行迭代，模型参数向量θ_n由式(9)进行估计：

式中：，a⁽⁰⁾和b⁽⁰⁾为模型参数初值。

迭代目标函数为

步骤2  在进行第l次迭代时，模型参数为a^(l)=[a₁^(l), a₂^(l), …, a_n^(l)]和b^(l)=[b₁^(l), b₂^(l), …, b_n^(l)]，对应的和可由式(11)计算：

当k＜0时，。

步骤3  执行第l+1次迭代，采用最小二乘法对式(12)分别估计出a^(l+1)和b^(l+1)：

式中：，d^-1为移位算子；, b₀为第一个模型参数。

步骤4  令l=l+1，重复步骤2和步骤3，直到满足|(J^(l+1)-J^(l))/J^(l)|＜δ或者l=L，其中δ和L分别为收敛指标和最大迭代次数。一旦迭代结束，即可得到最优模型参数。

需要注意的是，当噪声的频带靠近或者覆盖传感器输出信号的振铃频率时，模型参数的估计精度将会降低甚至于不能得到收敛的结果。原因是经验模态分解不能有效地消除这种类型噪声的影响，导致目标函数J的初值太大，在有限的迭代次数下不能满足迭代条件|(J^(l+1)-J^(l))/J^(l)|＜δ。

另一方面，在描述压力传感器的动态特性时，通常认为它是一个二阶线性模型^[19]。然而实际动态校准的结果很容易被一些不可控的因素影响，采用二阶线性模型表示其动态特性必然会增大动态特性参数的估计误差。为了解决这个问题，本文采用残余方差准则来估计模型的最优阶数。

定义残余方差为

式中：y为压力传感器的输出信号；为模型参数向量的估计值。

随着模型阶数的增加，残余方差的值逐渐减小。当残余方差趋近于稳定时，就可以得到模型的最优阶数。

1.4 动态特性参数的不确定度估计

得到压力传感器的数学模型之后，根据模型时域输出曲线和幅频特性曲线，分别计算出时频域动态特性参数^[5]。其中时域参数主要包括上升时间t_r、调节时间t_s和超调量σ；频域参数主要包括谐振频率ω_r和幅值误差为±10%的工作频带ω。时频域动态特性参数的定义如图 2所示。

采用激波管装置对压力传感器进行单次校准实验的时间一般都比较长，校准的成本也比较大，实际校准的重复性实验次数普遍比较少。假设进行c次重复性校准实验，根据本文提出的建模方法，分别得到时频域动态特性参数序列，序列的长度都是c。在参数序列样本量比较小的情况下，如果沿用传统的贝塞尔法计算每个参数的不确定度，很可能得到不可靠的评定结果，并且无法得到不同置信水平下的不确定度^[20]。

下面以工作频带为例，采用自助法对每个参数序列分别进行不确定度评定。

设c次重复校准得到的工作频带序列W为

从W中按照等概率可放回的方式抽样c次，得到自助样本W_b为

式中：ω_b(i)为W_b的第i个数据，i=1, 2, …, c。

计算自助样本序列W_b的均值为

重复自助采样B次，从原始工作频带序列中得到B个自助样本，构成矩阵：

式中：W_b为第b个自助样本，b=1, 2, …, B。

由B个自助样本得到相应的均值序列：

B通常是一个很大的数，因此ω_mean是一个大样本序列，称之为工作频带均值自助样本序列。将ω_mean从小到大排序，并且等分为Q组；采用直方图法建立ω_mean序列的自助概率函数为F_ω。

工作频带的估计真值ω_ev用加权均值表示为

式中：ω_mq为第q组的中值；F_ωq为第q组的自助概率，q=1, 2, …, Q。

在置信概率P下，工作频带真值估计区间为

式中：α=1-P为显著性水平，α∈[0, 1]；ω_α/2和ω_1-α/2分别为对应自助累加概率为和1-的工作频带值；ω_L和ω_U分别为工作频带估计区间的下界和上界^[21]。

定义工作频带的动态校准扩展不确定度为

相对不确定度表示为

2 实例分析

本实验使用北京长城计量测试技术研究所建立的激波管动态校准系统，如图 3所示。该激波管包含高压室和低压室2部分，长度分别为4 m和7 m，直径都为0.1 m，中间用铝膜片隔开。在进行动态校准实验时向高压室不断充入空气，当高低压室的差压到达铝膜片承受的最大压力时，膜片破裂。产生的激波向低压室的端面方向传播，在端面处形成反射激波阶跃压力。安装在端面上的被校准压力传感器受该压力的激励产生阶跃响应信号。其中被校准压力传感器由不锈钢的圆形载体支撑，且其安装位置为该载体的中心，对应于激波管横截面的中心位置。在动态校准过程中，传感器的输出会收到激波管端面的振动、温度的剧烈变化和二次激波冲击等的综合影响，如果仍然采用二阶线性模型表示其动态特性，必然会得到不准确的结果。

被校准压力传感器型号为PCB M102A02，灵敏度为2.67 mV/kPa，最大测量范围为690 kPa，铝膜片厚度为0.07 mm，数据采集系统的采样频率为5 MHz。得到被校准压力传感器输出如图 4所示。其中幅频谱由阶跃响应曲线的傅里叶变换得到，为了减小栅栏效应，幅频谱长度L_f和阶跃响应曲线长度L_t满足：2^i-1＜L_t＜L_f=2ⁱ，其中i为正整数。

从图 4(a)可以看出，压力传感器输出的振荡幅值随时间并不逐渐减小。其原因在于采用激波管系统对传感器进行动态校准时，受传感器敏感面与激波管端面不齐平、端面振动、破膜不全面导致激波幅值不稳定等因素的干扰，导致不规则的阶跃响应信号，如在大约2.1 ms处曲线出现明显的下降段，在2.6 ms之后仍然存在波动。

从图 4(b)的幅频特性曲线可以看出，在大约10⁰~10² kHz频率范围内的频谱幅值较大，并且存在明显的波动，其原因在于校准过程中受到低频噪声的干扰，如传感器安装不齐平导致的管腔效应等。这些现象会降低压力传感器的建模精度。在幅频特性曲线曲线最大峰值处对应的频率，称之为振铃频率。通常，用于动态测量场合的压力传感器的阻尼比很小，振铃频率值与传感器的谐振频频几乎相等，因此在压力传感器的动态校准时振动频率附近的频段对其谐振频率的校准结果影响较大，在进行信号预处理时应尽量保留该频段信号。

对压力传感器输出信号进行经验模态分解处理的结果如图 5所示。可以看出，原始输出信号分解成8个IMF分量(IMF1~IMF8)和一个残余分量。为了分离原始输出信号中的噪声分量，分别计算各IMF和原始输出信号的相关系数以及在振铃频率处的幅值比，如表 1所示。可以看出IMF1与原始信号的相关系数为0.942，远大于其他IMF的相关系数，说明IMF1包含原始信号中的大部分有用信息。从表 1的第3列可以看出，IMF1、IMF2和IMF3包含原始输出信号在振铃频率处99.653%的频谱幅值，其他IMF的幅值比都接近于0，判断其余IMF分量不包含有用信息的噪声分量。因此选择IMF1、IMF2和IMF3为重构传感器的输出信号，其余IMF分量予以剔除。

将IMF1、IMF2和IMF3相加，得到传感器输出的重构信号如图 6所示。可以看出经过经验模态分解预处理之后，原始输出信号中的低频噪声被剔除，特别是在2.1 ms处曲线的下降段被消除，在2.6 ms之后的曲线也趋于稳定。

得到重构信号之后，采用自适应最小二乘法对传感器进行建模。首先确定模型的最优阶数，分别计算不同阶数的模型残余方差和运行时间。其中残余方差由式(13)计算得到，结果如图 7所示。可以看出，残余方差值随模型阶数的变大而减小。当阶数大于4时，残余方差趋于稳定。另一方面，运行时间随着阶数的变大明显增长。综合2个参数的变化判断出模型的最优阶数为4。

确定了传感器数学模型的最优阶数后，模型参数a_i、b_i(i=1, 2, 3, 4)可以通过自适应最小二乘估计获得，进而得到传感器的差分方程模型为

对式(23)进行Z变换，并令z=[1+s/(2f)]/[1-s/(2f)]，其中f=5 MHz为采样频率。则传感器的连续传递函数为

根据求得的压力传感器传递函数，可以得到模型的归一化时域输出曲线与幅频特性曲线分别如图 8(a)、(b)所示。

对压力传感器重复校准6次，每次的实验结果都按照上述过程进行预处理和建模，得到时频域动态特性参数如表 2所示。在不同实验中得到的结果有所不同。需要对校准结果的差异进行分析，评定出各动态特性参数的不确定度。

以工作频带为例，采用自助法计算时频域动态特性参数不确定度。设置采样次数B=3 000，组数Q=30，得到的自助概率直方图以及自助累加概率与各组工作频带均值的关系分别如图 9和图 10所示。

从图 10可以计算出不同置信概率下工作频带的扩展不确定度。类似地，采用自助法对其他4个动态特性参数序列进行采样，计算自助累加概率和各组参数均值之间的关系，得到估计真值和不同置信概率下的扩展不确定度。在评定传感器的动态校准不确定度时，可以根据实际要求设置相应的置信水平，得到合适的不确定度数值。

采用贝塞尔法和蒙特卡罗法^[22]估计动态特性参数的估计真值及其扩展不确定度。其中贝塞尔法取包含因子为2；执行蒙特卡罗法时，根据重复实验得到传感器数学模型参数的估计均值和标准差，针对每个模型参数，采用MATLAB中的randn函数生成服从其估计均值和标准差的正态分布随机数，个数为100 000，得到模型参数的样本值，进而可以求得其估计真值和扩展不确定度。3种方法对于传感器动态特性参数真值和扩展不确定度的估计结果分别如表 3和表 4所示。从表 3可知，3种方法得到动态特性参数估计真值之间的相对误差小于1%，说明本文方法可以有效地估计时频域动态特性参数的真值。

对比表 4中的结果可以看出：一方面，采用传统贝塞尔法得到的扩展不确定度，比置信概率在100%下自助法得到的扩展不确定度大得多。说明在样本量比较少的情况下，如果仍然沿用贝塞尔法，得到的不确定度偏大。这是因为在采集数据样本量比较小的情况下，统计标准差无法准确地表征传感器动态特性的真实分散性。另一方面，对比本文方法与蒙特卡罗法，当置信概率与包含概率取值相同时，2种方法得到的结果比较接近，对应相对误差小于10%，说明了本文方法评定动态校准不确定度的准确性和有效性。

压力传感器动态特性参数的相对不确定度如图 11所示。对于频域动态特性参数，在不同置信概率下谐振频率的相对不确定度都最小，说明传感器重复校准对谐振频率的影响很小；工作频带的相对不确定度均大于2%，在置信概率为P=100%时相对不确定度高达7.12%，说明激波管在动态校准过程中存在一定的噪声，并且对传感器的工作频带影响较大。因此在压力传感器的动态校准过程中，应当尽量减小低频噪声的干扰，确保工作频带校准结果的可靠性。在时域动态特性参数方面，上升时间和调节时间的相对不确定度相差很小，都小于2%；超调量的相对不确定度略大，在P=100%时相对不确定度为3.14%。这很可能是由于每次实验中膜片的破裂情况不同，产生的激波到达低压室端面时的状态有所差异，导致传感器产生了不同超调量的响应信号。

3 结论

就激波管系统对压力传感器动态校准时的原理和特点，提出了一种压力传感器动态特性参数不确定度的评定方法。通过动态校准实验得到：

1) 采用基于经验模态分解和自适应最小二乘的传感器模型辨识方法，信号预处理有效地改善了输出信号的低频波动和曲线的不规则性。迭代使用最小二乘法估计得到压力传感器数学模型的最优参数，并根据残余方差准则估计出模型的最优阶数，得到了传感器的时频域动态特性参数。

2) 采用自助法对压力传感器的时频域动态特性参数进行不确定度评定，与贝塞尔法和蒙特卡罗法对比可知:一方面，本文方法弥补了贝塞尔法在小样本量情况下难以准确评定动态校准不确定度的缺点；另一方面，本文方法与蒙特卡罗法得到的扩展不确定度相对误差小于10%，验证了本文方法的准确性和有效性。

3) 在不同置信概率下，通过比较时频域动态特性参数的相对不确定度，得到压力传感器的工作频带和超调量受噪声干扰均比较大，在置信概率为100%时的相对不确定度分别为7.12%和3.14%。在压力传感器的动态校准过程中，主要应当考虑如何有效地消除低频噪声的干扰，保证激波状态的一致性，从而提高压力传感器动态校准结果的可靠程度。