动态力计量建模与分析方法综述

江文松，吴娅辉，李新良，王中宇，张力*

（1.中国计量大学计量测试工程学院，浙江 杭州 310018；2.航空工业北京长城计量测试技术研究所计量与校准技术重点实验室，北京 100095；3.北京航空航天大学仪器科学与光电工程学院，北京 100191）

0 引言

动态计量数据作为一种典型的科研数据，是国家质量基础的核心内容之一［1-2］。在计量领域，测量数据的正确建模与分析是处理复杂系统量值溯源与传递的关键。

计量需要选取合适的动态发生信号作为标准量的激励源，用以激发被校准传感器、测量装置、实验装置等工业仪器仪表的动态特性［3-4］。例如：在动态力校准中，既可选择电磁激振器、机械激振器等产生标准正弦激励信号来激发被校准仪表的动态频域特性，也可以选择落锤式冲击力校准装置产生标准半正弦脉冲激励信号来激发被校准仪表的瞬态响应特性。

对被校准仪表建立数学模型既能帮助计量检定人员了解其工作原理，选择合适的计量方式，也可为生产厂家开展仪表参数的辨识和校准提供量值溯源的可能［5-6］。例如：在建立好压力传感器的二阶参数化模型之后，采用激波校准法实现其参数的辨识和二阶模型响应特性的量化［7］。

动态力校准装置通过施加标准激励和可溯源的量值测量，就能通过辨识的数学模型进行动态时域特性和频域特性评价。例如：通过时域上的上升时间、峰值时间或频域上的穿越频率可以评价被校准仪表的响应速度；通过时域上的调节时间、超调量、或频域上的相位裕量可以评价被校准仪表的稳定性；通过稳态误差或动态降落可以判断被校准仪表的抗干扰特性等［8-12］。

由于对计量活动认识的局限性，测量过程不可避免地存在测量误差、建模误差、数据分析误差等，产生具有一定量值水平的测量不确定度［13］。通过这些评价指标，实现动态量值的立体计量。例如：Heisenberg在1927年针对量子物理学实验首次提出了测不准原理及概念，1993年国际不确定度工作组联合各个专业国际组织和各国计量部门，制定了《测量不确定度表示指南》（《Guide to the Expres⁃sion of Uncertainty in Measurement》，简称《GUM》），被世界各国及国际组织沿用至今，并被广泛应用于航空航天、电力电子等现代高端制造行业［14-15］。

因此，激励源设计、被校准仪表建模、动态时频域评价、测量不确定度评定等数据处理方法是开展动态校准、量值溯源与传递等的关键环节。本文以动态力校准为例，分别通过讨论激励源的分类形式、工作原理及模型，分析动态计量参数的建模方法，建立其计量数据的时频域评价体系，介绍动态校准数据的测量不确定度评定方法，系统地概述动态力计量的建模和分析理论，为动态力计量数据评价机制的建立提供理论基础。

1 动态力激励信号的量值表征

在动态计量领域，通常采用阶跃信号、脉冲信号、正弦信号、斜坡信号、抛物线信号等特征数据作为典型的激励源。在动态力校准实践中，前三种类型的信号较为常见，属于常规激励信号，通常利用这些激励开展不同形式动态量值的测量、溯源、校准及评价等计量活动，如图1所示。

图1 动态量值计量过程的建模与分析方法Fig.1 Modeling and analysis method of traceable measurement for the dynamic quantity

1.1 常规激励源

1.1.1 正弦激励信号

正弦激励信号r（t）属于周期变化的信号，在交流供电电源、动态力、动态压力、振动、动态温度等现场校准过程中较为常见，经常在动态校准过程中，用于信号的频谱分析，其数学表达式为

式中：A为信号的幅值；ω 为角频率；t为时间，s；在频域表达上，r（t）需要进行拉普拉斯变换，L（r（t））为对r（t）进行拉普拉斯变换，得到其频域的表征形式值R（s）；s为频率，Hz。

在动态校准过程中，将正弦信号作为激励，通过频率特性法计算被校准仪表的频率响应特性，实现其幅频特性、稳定性等动态过程指标的分析。例如：正弦激振器为正弦力激励源，输出周期频率可变的动态力，如图2所示。

图2 正弦激励源及信号发生原理Fig.2 Sinusoidal excitation and its signal generation principle

由于常规激励源具有激发动态特性的良好能力，研究人员近年来针对正弦激励源设计了动态力计量装置，并做了大量的研究。通过被校准力传感器、负载质量块和参考传感器来实现动态载荷的精确计量与校准［16］，动态力校准的矢量受力方程为

式中：m1为标准质量块质量，kg；m2为连接件的质量，kg。

正弦力校准装置最早由德国联邦物理技术研究院（Physikalisch-Technische Bundesanstalt，PTB）提出，利用波形发生器、功率放大器和高精度振动台产生特定频率和幅值的正弦力。2002年，T.Bruns研究员和R.Kumme教授在此基础上对引入气浮技术、减震技术和加速度传感器网络等，实现横向偏摆、惯性质量、加速度分布不均匀等影响因素的干扰消除，在10 Hz～2 kHz 的频率范围内实现了100，800 N和10 kN的正弦力，测量不确定度达0.4%［19］。法国国家实验室（Laboratoire national de métrologie et d'essais，LNE）和西班牙计量中心（Centro Español de Metrología，CEM）及北京长城计量测试技术研究所（Changcheng Institute of Me⁃trology&Measurement，CIMM）等也基于该原理设计了系列正弦力校准装置（如图3 所示），并实现国际比对［20］。

图3 不同型号正弦力计量装置Fig.3 Different sinusoidal force measuring devices

1.1.2 冲击激励信号

冲击信号持续时间极短，其数学表达式为

式中：H为冲击幅值；ε为脉宽。通常当H=1 且脉宽逼近时，即为常见的单位理想冲击函数δ(t)，其面积为脉冲强度。据此可知δ(t)信号的拉普拉斯变换为1，具有瞬间无穷大的能量，因此这种信号无法在实际工程中获得，如图4所示。

图4 冲击信号与单位理想冲击函数Fig.4 Impulse signal and the unit ideal impulse function

在计量实践中，通常把ε远小于被测传感器响应时间的窄脉冲半正弦信号近似为δ(t)函数。例如：落锤式发生装置可以近似为冲击力激励源，能够输出不同脉冲宽度和脉冲幅值的动态力，如图5 所示。

图5 冲击激励源及信号发生原理图Fig.5 Impulse excitation and its signal generation principle

冲击力校准装置由激励源和测量单元组成。各国也基于该激励原理实现冲击力计量装置的研制，如图6 所示。例如：中国计量科学研究院（National Institute of Metrology，China，NIM）20 世纪80年代就开展了500 N～200 kN 冲击力校准典型装置研制，其在0.6 ms 的脉冲宽度的测量不确定度为2%，此时自由度取k=3［21］；印度国家物理实验室（National Physical Laboratory India，NPLI）在2012年联合德国PTB 研制了1 kN～1 MN 的冲击力校准装置，实现了100，200，500，1000 kN 等冲击力的溯源，重复性误差接近1%［22］；CIMM 研制的冲击力校准装置，能够实现0.5 ms脉宽、180 kN冲击力的溯源，重复性测量误差不大于1%［23］。

图6 不同型号冲击力计量装置Fig.6 Different impact force measuring devices

1.1.3 阶跃激励信号

阶跃信号为通过瞬间突变后保持的信号，其数学模型可表示为

式中：R0为常量，当R0=1 时称为单位阶跃信号ε(t)，在进行动态特性分析时，可将其进行拉普拉斯变换为R（s）。

阶跃信号的发生比较特殊，既需要具备瞬态突变特征，又需要有量值保持能力。例如：激波管的膜片破裂、脆性梁的断裂、电源突然通断、负载突变等均近似地认为会输出阶跃信号。另外，阶跃信号是传感器或测量系统时域性能分析的常用激励，因此被广泛应用于动态性能分析。在动态力校准中，脆性梁的瞬间断裂产生的负载通常为阶跃力，用来产生标准阶跃激励，如图7所示。

图7 阶跃激励源及信号发生原理Fig.7 Step excitation and its signal generation principle

阶跃力校准主要以卸荷的方式产生，例如，德国PTB 实现了20 kN 和250 kN 的阶跃力发生，北京航空航天大学黄俊钦教授实现了196 kN 阶跃力校准［24］；北京理工大学也利用该方法实现了100 kN 负阶跃力校准［25］；北京长城计量测试技术研究所采用快速卸载方式在10～100 μs 的上升时间实现了1.2 MN 负阶跃力校准等［18］。但相比于正弦力和冲击力校准方式，由于阶跃信号的保持难度较大，因此在计量实践中并不多见。

1.2 其他激励源

1.2.1 斜坡激励信号

斜坡信号是随时间t线性增长的信号，其数字表达式为

式中：v0为常量，当v0=1时描述的激励称为单位斜坡信号。在进行动态特性分析时，可将其进行拉普拉斯变换为R（s）。

1.2.2 抛物线信号

抛物线信号随时间t以等加速度增长，因此也称为等加速度信号，其数学表达式及拉普拉斯变换形式可以表示为

式中：a0为常量，当a0=1 时，称为单位等加速度激励信号。

由于斜坡信号和抛物线信号作为激励源不仅不容易准确实现信号发生与放大，还很难激发被校准传感器或仪器的动态性能，因此在动态校准中并不常见。

2 动态力校准参数的建模与分析方法

2.1 动态系统建模

当确定了激励源之后，建立被校准仪表的模型是动态计量的关键环节。为此需要利用被校准仪表上的激励和响应信号，对其模型进行辨识。

针对原理已知的白箱系统，通常对其各个环节所遵循的基本物理规律建立微分方程组，通过消除中间变量和归一化处理，建立传感器输入与输出之间的微分方程（由输出量的各阶导数和输入量的各阶导数及系统结构参数构成），从而得到动态系统的数学模型。

令n阶线性定常系统，其数学模型可用n阶线性定常微分方程描述为

图8 其它激励源Fig.8 Other excitations

式中：c（t）为输出量；r（t）为输入量，即激励源信号；ai和bj分别为系统结构参数，i=1…n，j=1…m。可见，动态系统的数学模型能揭示不同类型传感器物理系统的共性内在特征，给动态校准过程带来便捷。

为了求得被校准仪表的传递函数G（s），在零初始条件下，对式（7）进行拉普拉斯变换，即

针对黑箱模型，通常还需要在多次重复实验的条件下，采用参数辨识的方式对激励和响应信号的卷积、归一化、聚类、模态分解或正则化等进行信号处理，实现黑箱模型参数及微分关系的辨识，从而建立被校准仪表的动态模型。

以动态力计量建模为例，在力传感器动态建模过程中，为了溯源更为准确的动态力，通常还会在式（2）的基础上增加修正因子。例如，采用多传感器融合的方式修正加速度分布不均引起的测量误差，建立的时域修正模型为

式中：m3为传感器的惯性质量，kg；ω为振动台的激振频率，Hz；l为负载质量块的长度，mm；ρ为负载质量块的密度；E 为杨氏模量；n为质量块上表面参考加速度计的数量。

为了便于动态特性分析，力传感器通常可以近似为二阶模型。其中，压电材料具有弹性和阻尼特征，弹性结构的劲度系数为k，阻尼结构的阻尼系数为f，顶部质量块的质量为m1，基部质量块的质量为m2，在动态力F（t）的作用下，建立质量块m1和m2及其在空间上产生的相对时变位移x1（t）和x2（t）之间的受力平衡关系，从而得到力传感器的开环传递函数G（s）为

根据传递函数，获得频率特性表达式，即

通过式（11）可得到G（s）的幅频特性A(ω)和相频特性φ(ω)。因此，只要确定了被校准力传感器的频率特性，就能辨识其关键参数水平。力传感器动态特性由其传递函数确定，采用上述的正弦、阶跃或脉冲等常规激励作用，即可求解传递函数模型的各参数，实现动态力测量过程的动态建模。

2.2 时域评价方法

在动态计量中，被校准传感器的时域动态指标通常采用阶跃激励r（t）作用，评价被校准仪表的跟随性能及稳定性能。典型的被校准传感器的阶跃响应曲线如图9所示，其中，反映响应特性的指标包括上升时间tr、峰值时间tp、超调量σ%、调节时间ts、振荡次数、动态降落、恢复时间及稳态误差等。

图9 阶跃响应曲线及性能指标Fig.9 Step response curve and performance index

设某被校准传感器为二阶线性定常系统，即n=2，在阶跃激励r（t）作用下的输出c（t）的拉普拉斯变换为

式中：ξ为被校准传感器的阻尼系数，Φ(s)为被校准传感器的传递函数，根据其分母s2+ 2ξωns+即可求得两个特征根，分别为

通常根据ξ在[0，+ ∞)上取值的不同，特征值s1，2有可能为实数根、复数根或重根。被校准传感器动态时域性能评价指标分别如下：

1）上升时间

上升时间tr是系统输出响应从稳态值的10%上升至稳态值的90%所需时间。因此，根据式（8）可求得c（t）为

由定义即可求出上升时间tr为

2）峰值时间

峰值时间tp是指传感器或系统的输出响应自零时刻开始到首次达到峰值所经历的时间。对于二阶系统，此时的c（t）在tp时刻的一阶导数为0，即

于是可求得峰值时间为

3）超调量

超调量σ%为传感器或系统的输出响应超出稳态值的最大偏离量占稳态值的百分比，即

由式（18）可知，σ%越小，被校准仪表的动态响应越平稳。

4）调节时间

调节时间ts是从零时刻开始，被校准仪表的输出响应进入并保持在某一误差带Δ之内所需的时间，通常Δ取±5%（或±2%）。

显然，误差带越小，对被校准仪表的计量特性要求就越高；调节时间越小，表明被校准仪表的动态响应过程越短，即响应越快速。

5）稳态误差

从定性的角度看，被校准仪表的输出量围绕稳态值摆动的次数越少、摆动的幅度越小，表明系统稳定性越好。这种摆动次数和摆动幅度通常用稳态误差ess来定量表示。稳态误差是指系统期望值r0（t）与实际稳态值r（t）之间的差值，即

由式（19）可知，ess越小，系统稳态精度越高。

6）动态降落值

动态降落值ΔCmax用于衡量被校准仪表的抗扰动性能，通过对稳态过程施加一个标准负阶跃激励，就能测量出被校准仪表响应的最大降落值ΔCmax，从扰动激励到最大降落值之间的过渡时间称为降落时间tm。

当采用脉冲信号、正弦信号等其它激励源时，还需要通过被校准传感器的输出计量数据的峰谷值、脉冲宽度、平均值等特征指标来评价其动态计量性能。

2.3 频域评价方法

以典型二阶系统为例，无论是时域评价方法还是频域评价方法，关注的重点均为被校准仪表输出响应的快速性和稳定性。区别于时域评价方法，常用频率特性中的穿越频率ωc来反映系统的快速性，用相位裕量γ来衡量系统的相对稳定性。

令s=jω，根据式（8）可知，被校准仪表的幅频特性A(ω)可表示为

相频特性φ(ω)为

1）相位裕量

将穿越频率ωc代入式（16）并令A（ωc）=1，即可得到相位裕度γ为

对于典型二阶系统，相位裕量γ只与系统的阻尼比ξ有关，具体地，当ξ越大则γ越大，被校准仪表的平稳性及相对稳定性越高，如图10所示。

图10 典型二阶系统幅频特性曲线Fig.10 Amplitude-frequency characteristic curves of typical second-order systems

2）穿越频率

将时域分析中的调节时间ts，令误差带Δ=±5%，则可得到ts与穿越频率ωc之间满足

由此可知，调节时间ts，ωc与γ之间满足函数关系。具体地，在相位裕度γ不变时，穿越频率ωc越大，调节时间ts越短，被校准传感器及系统的响应越快。

另外，在进行典型高阶测试系统的频域特性评价时，谐振频率和阻尼比系数也常常被用于判断其系统的跟随性能及稳定性能。由于二者均可与相位裕量、穿越频率建立算法转换联系，因此在功能原理上类似，在此不再赘述。

3 测量不确定度评定

测量不确定度是一个体现测量值分散性的参数，它由基于概率论和统计理论的A类评定和基于先验知识的B类评定得到，是两类传统评定方法。

3.1 两类传统评定方法

3.1.1 A类评定方法

假设某被测量x，在相同测量条件下独立重复测量n次，观测值为xi（i=1，2，…，n），任意一次测量值xk的A 类标准不确定度用其实验标准差s（xk）表示，因此，采用统计分析的方法求得的标准差即为A 类评定方法，以贝塞尔公式法为例，标准差为

除此之外，标准差还可采用极差法、最大残差法、最大方差法、彼特斯法等进行求解。

3.1.2 B类评定方法

该方法是基于信任度的主观概率方法。当输入量Xi的估计值xi不由重复观测得到时，估计方差或标准不确定度可根据Xi的可能变化的全部有关信息的判断来评定，这些信息来源可能包含：以前的测量数据、对有关材料和仪器特性的了解、生产厂提供的技术说明书、校准证书或其他证书提供的数据、手册给出的参考数据的不确定度、规定实验方法的国家标准或类似技术文件中给出的重复性限r或复现性限R。

在动态力校准过程中，误差主要来源于各个观测值。其中，冲击力校准装置主要测量误差源包括有效碰撞质量的测量误差、激光光轴与运动方向的夹角误差、激光干涉仪的测量误差、落锤的加速度分布不均匀引入的误差、落锤的横向偏摆引入的误差、二次仪表误差和重复性测量误差等；正弦力校准装置的测量误差主要来源于等效质量的测量误差、激光光轴与运动方向的夹角误差、激光干涉仪的测量误差、质量块的加速度分布不均匀引入的误差、激振器的横向偏摆引入的误差、激振器波形失真和重复性测量误差等。其中，重复性测量误差可采用A 类评定方法求解标准差，其他误差来源可采用B类评定方法求解标准差。

3.2 蒙特卡洛法

传统评定方法仅适用于高斯分布模型，当测量模型具有偏导、积分等复杂函数时，测量不确定度的传统评定方法很难解决。为此，ISO/IEC Guide 98-3系列标准规定了蒙特卡洛法评定测量不确定度的实施步骤。

已知测量模型为Y=f（X），抽样样本X=（x1，x2，…，xN），其概率分布函数为px(ξ)，对应输出模型样本为Y=（y1，y2，…，yN），对应的概率分布函数为py(ξ)。蒙特卡洛法的具体步骤如下：

1）选择蒙特卡洛仿真试验次数

蒙特卡洛仿真试验次数M决定输出量的样本容量。M值越大，样本容量越大，越接近于输出量的真实情况。通常M的取值使包含区间达到95%及以上为宜。

2）确定输入样本的概率密度函数

蒙特卡洛法的评定结果依赖于测量模型，因此各输入量的概率密度函数是评定结果的约束。由于各输入量的概率密度函数依赖于数据本身规律，因此可以通过直接法、贝叶斯法、最大熵原理等统计分析方法来设定X的概率密度函数。

3）输入样本的概率密度函数抽样

根据上述概率密度函数，对输入样本X进行M次随机抽样，得到M×N维输入样本的伪随机数矩阵VM×N。

4）输出样本值的计算

将上述矩阵VM×N代入测量模型Y=f（X）中，即可得到对应的M组输出量WM×N。

5）输出评定结果

根据输出量Y的抽样WM×N，计算输出量Y的数学期望和标准差，其中数学期望即为Y的估计值，标准差即为标准测量不确定度；通过分布函数py(ξ)得到Y的包含区间和置信概率。

3.3 其他评定方法

由于传统评定方法依据大数定律和中心极限定理，当先验信息缺乏或测量数据为小样本量时，GUM 因缺乏理论基础会出现评估偏差和置信度下降的问题。为此，依据先验知识的贝叶斯理论、灰模型法及自助法等非统计分析方法被用于测量不确定度的评定上。

3.3.1 贝叶斯理论方法

设测量模型的参数为θ，测量数据样本集为X=（x1，x2，…，xN），该样本的似然函数为L（x|θ），采用贝叶斯理论可计算后验分布，即

由此即可确定后验分布的概率密度函数，测量结果取n次测量数据的平均值μ，将后验分布密度函数的数学期望表示最佳估计值μ̂，其标准差表示标准不确定度u，即

由此即可实现测量不确定度的评定。

3.3.2 灰模型法

设被测量重复测量n次，测量数据为X=(x(k)|k= 1，2，…，n)，首先，将数据序列X按升序的方式排列成序列X（0），利用灰色模型将X（0）一次累加生成（1-AGO）得到新序列X（1），即

同时，定义

对式（28）取最大值Δmax，表示测量序列累加曲线与X（0）极值连线之间的偏差最大值，即

若用单次测量值xi作为测量结果，则标准不确定度可以表示为

式中：c为灰度系数，由灰色GM（0，N）模型求得。

由此可以看出，贝叶斯法侧重于建立后验分布规律，灰模型法侧重于小样本数据的挖掘，这两种方法都旨在解决数据规律未知情况下的不确定度评价。但鉴于这类新方法在适用场合上的特殊性和对数据的二次加工处理，通常只会作为两类传统评定方法和MCM方法的补充。

3.4 不确定度合成与表示方法

根据测量模型Y=f（X），测量的估计值ŷ为

根据上述方法求得的标准不确定度，被测量估计值的合成标准不确定度uc（y）表示为

式中：xi为X的估计值；u（xi）为估计值xi的标准不确定度；∂f/∂xi为y与xi之间的相关系数。当确定了各标准不确定度分量u（xi）的自由度vi，即可利用式（33）计算合成标准不确定度的自由度v，即

为了保证扩展不确定度是被测量可能值包含区间的半宽度，扩展不确定度U由合成标准不确定度uc乘以包含因子k得到，即

测量不确定度的表示方式有很多，最为常见的有以下两种情况：

1）当不确定度用合成标准不确定度表示时，表示形式为

2）当不确定度用扩展不确定度表示时，表示形式为

通过上述过程，即可实现测量不确定度的评价。

4 结论

系统性地综述了动态力计量的建模和分析方法，探讨了激励源、参数建模、数据评价、测量不确定度评定等计量基础理论在实际计量工程案例中的应用。针对动态力校准过程中标准激励发生方式，探讨了各类激励信号源的原理、数学表示形式，简要介绍了不同量值溯源中的激励装置。针对被校准仪表的黑白箱模型，以动态力校准为例，分别介绍了基于物理环节的微分-传递函数建模方法和基于卷积、归一化、聚类、模态分解或正则化等信号处理的参数辨识方法。建立动态力校准模型的时域评价体系和频域评价体系，实现动态力校准模型的动态性、跟随性和稳定性的综合评价；系统地建立了动态力校准数据的测量不确定度评价方法体系，实现测量数据的可溯源。