CN104268414A - 一种故障自修复的鲁棒的煤矿巷道剪应力采集方法 - Google Patents

Abstract

一种故障自修复的鲁棒的煤矿巷道剪应力采集方法。煤矿现场通常用杆体来检测巷道的应变程度，主要测量杆体的剪应力分布，一般采用在杆体上安装应变片的方法，但是应变片数量多且容易损坏，而应变片数据采集都是自动采集，无人员实时监控，所以应变片发生故障损坏时不能被及时发现，得到错误的弯矩拟合函数和剪应力分布，这对于煤矿巷道观察被测杆体应变特性是非常不利的；其次，在煤矿现场安装杆体需要用特种水泥将其固定在巷道里，应变片发生故障损坏时将无法更换新的应变片，如不进行故障修复，少数应变片故障将造成杆体的浪费。

Description

一种故障自修复的鲁棒的煤矿巷道剪应力采集方法

技术领域

本发明专利涉及一种故障自修复的鲁棒的煤矿巷道剪应力采集方法，属于测控和自动化技术领域。

背景技术

剪应力是指物体由于外因(受力、湿度变化等)而变形时，在物体内各部分之间产生相互作用的内力，来抵抗外因的作用，并力图使物体从形变后的位置回复到形变前的位置。对于杆体，准确掌握整个长度范围内的剪应力分布，是研究整个杆体结构特性的重点。

煤矿巷道中，对于一些大型结构(桥梁等)，通常测量内置杆体的整体形变来反映外部结构的形变，准确掌握大型结构的形变，对于及时做好预防措施、减少安全隐患具有重要意义。因此，准确检测杆体的形变，特别是准确掌握杆体剪应力分布至关重要。

测量杆体的剪应力分布，通常采用在杆体上安装应变片的方法，由于杆体一般用于测量大型结构，所以杆体比较长，所需要的应变片数量很大。由于应变片体积小且薄，容易在煤矿巷道的环境中损坏，而应变片数据采集都是自动采集，无人员实时监控，所以当应变片发生故障损坏时不能被及时发现，故障应变片的错误采集值仍参与计算，得到错误的弯矩函数拟合和剪应力分布，因此，工作人员检查结构时看到错误的杆体性质，会相应的采取错误的措施，影响工作效率，并且杆体一旦浇筑进大型结构中无法拿出，如果出现故障不能解决，那么杆体检测失去了意义。本发明提出一种故障自修复的鲁棒的煤矿巷道剪应力采集方法，通过实时监测所有应变片是否有故障，当应变片发生故障时对其进行相应的采集数据剔除或者进行误差模型补偿，使得被测杆体在无人监控的环境下实现故障自修复，保证弯矩拟合函数的准确性，得到准确的剪应力的函数，从而获得整个被测杆体上的剪应力分布。

发明内容

本发明解决的问题是：为了克服工业中杆体上应变片发生故障时检测出错误的杆体剪应力分布，影响杆体结构性质判断的缺点，结合以上背景和需求，本发明提供一种故障自修复的鲁棒的煤矿巷道剪应力采集方法，通过检测出有故障的应变片，对其采集数值进行误差补偿，用补偿后的采样数值进行拉伸应变、弯曲应变、弯矩的计算，然后对每个测量点的弯矩进行函数拟合获得整个被测杆体上的弯矩函数，由弯矩拟合函数得到沿测点长度范围剪应力的函数，根据这个函数可以很容易得到杆体上每个测点处的剪应力。

本发明的技术解决方案是：

1、一种故障自修复的鲁棒的煤矿巷道剪应力采集方法，其特征在于，包括

以下步骤：

步骤1：设煤矿巷道剪应力测试杆体长度L，在该杆体上贴N对应变片，每对应变片之间的间隔d＝L/(N-1)，每对应变片正对贴在杆体上，设上面一层应变片的编号为X_i1，正对面一层应变片编号为X_i2(i＝0、1……N-1，i为第i对应变片)，将该杆体安装在对其进行弯曲、拉伸的相关设备上；

步骤2：误差模型的输出值的评估：

步骤2.1：利用相关设备对装有2N个应变片的煤矿巷道剪应力测试杆体以时间间隔t₁分别进行C次拉伸应变和C次弯曲应变，采集并记录弯曲应变、拉伸应变数值，记为ε_i1 ^拉伸，ε_i2 ^拉伸，ε_i1 ^弯曲，ε_i2 ^弯曲(i＝0、1……N-1)，此处时间间隔t₁时间间隔较短，模拟的是煤矿巷道应变突变的情况下的数据采集；

步骤2.2：分别计算煤矿巷道剪应力测试杆体模拟煤矿巷道应变突变的情况下每一对应变片的拉伸应变和弯曲应变数值误差δ_ij ^拉伸＝ε_i1 ^拉伸-ε_i2 ^拉伸，δ_ij ^弯曲＝ε_i1 ^弯曲+ε_i2 ^弯曲(i＝0、1……N-1，j＝0、1……C-1)；

步骤2.3：利用相关设备对装有2N个应变片的煤矿巷道剪应力测试杆体以时间间隔t₂分别进行C′次拉伸应变和C′弯曲应变，采集并记录弯曲应变、拉伸应变数值，记为(i＝0、1……N-1)，此处时间间隔t₂时间间隔较长，模拟的是煤矿巷道应变变化缓慢的情况下的数据采集；

步骤2.4：分别计算设煤矿巷道剪应力测试杆体模拟煤矿巷道应变变化缓慢的情况下每一对应变片的拉伸应变和弯曲应变数值误差(i＝0、1……N-1，j′＝0、1……C′-1)；

步骤2.5：分别计算设煤矿巷道剪应力测试杆体每一对应变片拉伸应变和弯曲应变应的数值综合工况下误差平均值，(i＝0、1……N-1)，该平均值即为煤矿巷道综合工况下的应变误差模型输出值的估算值，该估算值即为应变偏差补偿的数据源；

步骤3：将剪应力测试杆体安装在煤矿巷道里，进行煤矿巷道综合工况下的数据采集。

步骤4：根据采集到的数据判断应变片是否故障：

步骤4.1：每隔时间t对杆体上的应变片采集数值进行读取；

步骤4.2：读取杆体上2N个应变片的采集数值ε_i1，ε_i2(i＝0、1……N-1)；

步骤4.3：如果连续三次ε_i1或者ε_i2值为0，则认为应变片发生故障，进入步骤5；

步骤4.4：如果ε_i1，ε_i2(i＝0、1……N-1)都不为0，即应变片未出现故障，则此时有效应变片对数N′＝N，进入步骤7；

步骤5：应变片故障类型判断：

步骤5.1：如果应变片X_i1出现故障，应变片X_i2正常工作或者应变片X_i2出现故障，应变片X_i1正常工作，进入步骤6.1；

步骤5.2：如果应变片X_i1和X_i2都出现故障，则记录应变片X_i1和X_i2都出现故障的应变片对数为P，进入步骤6.2；

步骤6：根据故障类型进入相应的应变片故障处理流程：

步骤6.1：对于故障应变片的采集数值进行误差补偿，如果X_i1出现故障，则计算拉伸应变时对X_i1应变片进行补偿，令则计算弯曲应变和弯矩时令则可以得到第i个应变片的弯矩：M_i ^弯矩＝π·E·R³·(ε′_i1-ε_i2)/8，其中π是圆周率，E为弹性系数，R为杆体的截面半径；如果X_i2出现故障，则计算拉伸应变时对X_i2进行补偿，令则计算弯曲应变和弯矩时令则可以得到第i个应变片的弯矩：M_i ^弯矩＝π·E·R³·(ε_i1-ε′_i2)/8，进入步骤8；

步骤6.2：如果第i对两个应变片都是有故障的，则剔除这一对应变片的采样数值，将剩下的正常工作的应变片作为有效测点，则修正后有效应变片对数N′＝N-P进入步骤7；

步骤7：被测杆体应变参数计算，E为弹性系数，R为杆体的截面半径，A为杆体横截面积，i′为故障处理之后有效应变片的序号，ε_i′1,ε_i′2为步骤6中进行误差补偿之后的应变值：

步骤7.1：计算被测杆体上每个应变片测点拉伸应变参数，拉伸应变为拉伸应力为(i′＝0、1……N′-1)；

步骤7.2：；计算被测杆体上每个应变片测点弯曲应变参数，弯曲应变为弯曲应力为(i′＝0、1……N′-1)；

步骤7.3：；计算被测杆体每个有效应变片点弯矩：

M_i′＝π·E·R³·(ε_i′1-ε_i′2)/8，(i′＝0、1……N′-1)    (1)

步骤8：对被测杆体的N′对有效应变片的弯矩M_i′进行函数拟合，在被测杆体由应变片分割的N′-1个区间段上使用三次样条插值拟合获取被测杆体上的弯矩函数，每个区间段上的拟合函数的一阶导数和二阶导数均为连续函数，则获得的拟合弯矩函数曲线是光滑曲线，具体拟合步骤如下：

步骤8.1：被测杆体的N′对有效应变片，相邻两个应变片之间是一个区间，所以有(N′-1)个区间，设每个区间段[z_k,z_k+1](k＝0,1,...,N′-2)，k表示(N′-1)个区间的序号。设第k个区间[z_k,z_k+1]的三次样条插值拟合函数为S_k(x)，S_k(x)表示距离被测杆体端点为x处的弯矩，即：

S_k(x)＝a_k+b_k(x-z_k)+c_k(x-z_k)²+d_k(x-z_k)³    (2)

其中x是距离被测杆端点的距离，z_k＝k×d(k＝0,1,...,N′-2)是当应变片i＝k序号时，该应变片距离被测杆端点的距离，a_k,b_k,c_k,d_k为三次样条插值拟合函数需要求解的参数。

根据三次样条插值拟合函数的要求，其插值条件为：

当x＝z_k时，S_k(z_k)＝M_k    (3)

其中，M_k是当i′＝k时M_i′的值，M_i′的值从公式(1)得到。

三次样条插值拟合函数的连续性条件为：

$\underset{x\→z_k}{\lim }{S}_k\left(x\right)=S_k\left(z_k\right)---\left(4\right)$

三次样条插值拟合函数的一阶导数连续条件为：

$\underset{x\→z_k}{\lim }{S_k}^{\′}\left(x\right)={S_k}^{\′}\left(z_k\right)---\left(5\right)$

其中，S_k′(x)是S_k(x)的一次导数。

三次样条插值拟合函数的二阶导数连续条件为：

$\underset{x\→z_{i^{\′}}}{\lim }{S_k}^{\′\′}\left(x\right)={S_k}^{\′\′}\left(z_k\right)---\left(6\right)$

其中，S_k″(x)是S_k(x)的二次导数。

以上插值条件和连续性条件总共能够得到(4N′-2)条等式。而三次样条插值在每一个区间上有4个未知参数，因此总共有4N′个待求解未知参数。被测杆体两端的端点位置没有受到让它们弯曲的力，因此端点条件为自由边界，所以边界端点约束提供了两个等式，为：

通过以上4N′条等式，能够求解出S_k(x)中4N′个未知参数a_k,b_k,c_k,d_k(k＝0、1……N′-2)。

步骤8.2：求解三次样条插值拟合函数为S_k(x)。

步骤8.2.1：将拟合函数S_k(x)的二阶导数作为待定参数，令T_k＝S_k″(x)；

步骤8.2.2：因为S_k(x)是三次多项式，所以其二阶导数是一次多式，

令 $S_k^{\′\′}\left(x\right)=\frac{T_{k+1}-T_k}{h_k}(x-z_k)+T_k,x\∈{[z}_k,z_{k+1}]---\left(8\right)$

其中，h_k＝z_k+1-z_k。

对二阶导数积分，得到一阶导数：

$S_k^{\′}\left(x\right)=\frac{T_{k+1}-T_k}{{2h}_k}{(x-z_k)}^2+T_k(x-z_k)+p_k---\left(9\right)$

其中，p_k为常量参数。

对一阶导数积分，得到函数S_k(x)：

$S_k\left(x\right)=\frac{T_{k+1}-T_k}{{6h}_k}{(x-z_k)}^3+\frac{T_k}{2}{(x-z_k)}^2+p_k(x-z_k)+q_k---\left(10\right)$

其中，q_k为常量参数。

将x＝z_k带入公式(9)中，可以得到S_k(z_k)＝q_k，q_k即当x＝z_k时S_k(x)的值，根据公式(3)可知S_k(z_k)等于序号为k的应变片的弯矩值，即M_k，因此q_k是已知量。

因为拟合曲线平滑时端点处函数值相等，所以

S_k(z_k+1)＝S_k+1(z_k+1)＝q_k+1(11)将公式(11)代入公式(10)中，可以得到：

$S_k\left(z_{k+1}\right)=\frac{T_{k+1}-T_k}{{6h}_k}{(z_{k+1}-z_k)}^3+\frac{T_k}{2}{(z_{k+1}-z_k)}^2+p_k(z_{k+1}-z_k)+q_k=q_{k+1}---\left(12\right)$

由公式(12)得到

$p_k=f[z_k,z_{k+1}]-\frac{T_{k+1}+{2T}_k}{6}{h}_k---\left(13\right)$

其中， $f[z_k,z_{k+1}]=\frac{q_{k+1}-q_k}{z_{k+1}-z_k}---\left(14\right)$

将公式(13)带入公式(10)中可以得到三次样条插值拟合函数为：

$S_k\left(x\right)=\frac{T_{k+1}-T_k}{{6h}_k}{(x-z_k)}^3+\frac{T_k}{2}{(x-z_k)}^2+\left(f\right[z_k,z_{k+1}]-\frac{T_{k+1}+{2T}_k}{6}{h}_k)(x-z_k)+q_k---\left(15\right)$

公式(15)中，T_k,T_k+1是未知量，需要继续求解。

步骤8.2.3：确定二阶导数待定参数T_k。

因为拟合曲线一次导数平滑所以端点处相邻两个函数值相等：

S′_k-1(z_k)＝S′_k(z_k)    (16)

将公式(16)代入公式(9)中，即：

$\begin{array}{c}\frac{T_k-T_{k-1}}{{2h}_{k-1}}{(z_k-z_{k-1})}^2+T_k(z_k-z_{k-1})+f[z_{k-1},z_k]-\frac{T_k+{2T}_{k-1}}{6}{h}_{k-1}\\ =f[z_k,z_{k+1}]-\frac{T_{k+1}+{2T}_k}{6}{h}_k\end{array}---\left(17\right)$

化简公式(17)，可以得到：

$\frac{h_{k-1}}{h_{k-1}+h_k}{T}_{k-1}+{2T}_k+(1-\frac{h_{k-1}}{h_{k-1}+h_k})T_{k+1}=6\frac{f[z_k,z_{k+1}]-f[z_{k-1},z_k]}{h_{k-1}+h_k}---\left(18\right)$

令：

${\mathrm{\μ}}_k=\frac{h_{k-1}}{h_{k-1}+h_k},{\mathrm{\λ}}_k=1-{\mathrm{\μ}}_k,$

$d_k=\frac{f[z_k,z_{k+1}]-f[z_{k-1},z_k]}{h_{k-1}+h_k}(k=1......N^{\′}-2),d_0=f[z_0,z_1],d_{N^{\′}-1}=f[z_{N^{\′}-2},z_{N^{\′}-1}]$

则公式(18)可以化简为：

μ_kT_k-1+2T_k+λ_kT_k+1＝6d_k    (19)

又因为边界条件公式(7)可以得到

$\left\{\begin{array}{c}{S_0}^{\′\′}\left(z_0\right)=0\\ {S^{\′\′}}_{N^{\′}-1}\left(z_{N^{\′}-1}\right)=0\end{array}\right.$ 即 $\left\{\begin{array}{c}{T}_0=0\\ T_{N^{\′}-1}=0\end{array}\right.---\left(20\right)$

可以得到求解T_k(k＝1,2......N′-2)的矩阵方程，(N′-2)个方程可以求解出(N′-2)个T_k的值：

步骤8.3：将求得的T_k(k＝0,1......N′-1)代入三次样条插值拟合函数S_k(x)公式(15)中，获得每一个区间段[z_k,z_k+1](k＝0,1,...,N′-2)对应的弯矩函数：

$S_k\left(x\right)=\frac{T_{k+1}-T_k}{{6h}_k}{(x-z_k)}^3+\frac{T_k}{2}{(x-z_k)}^2+\left(f\right[z_k,z_{k+1}]-\frac{T_{k+1}+{2T}_k}{6}{h}_k)(x-z_k)+q_k---\left(22\right)$

步骤9：沿测点长度范围剪应力分布函数：

步骤9.1：对公式(22)的函数S_k(x)进行求导运算得到S′_k(x)，根据拟合的弯矩函数的导数S′_k(x)得到每个区间段剪应力分布函数：

τ_k(x)＝4·S′_k(x)/(3π·R²)    (23)

其中x是距离被测杆端点的距离，R是杆体横截面的半径；

步骤9.2：将被测杆体上任意一点距离杆体端点的距离z代入公式(23)中，可以得到该点的剪应力。

本发明的有益效果包括：

(1)本发明及时对应变片进行故障检测以及相应的故障处理，该发明将故障分为一对应变片中有一片发生故障及一对应变片中两片均发生故障，根据不同的故障类型采用不同的故障处理流程，保证了整个杆体剪应力分布的准确性，避免了因为应变片故障导致整个杆体应变检测失效的问题，延长了杆体检测的使用寿命；

(2)本发明采用自适应的故障处理策略，当判断故障类型为一对应变片中只有一个应变片发生故障的情况时将采用本发明中提出的煤矿巷道应变误差模型偏差补偿法对故障进行修复，若判断为一对应变片都发生故障的情况采用剔除该应变片数据，将剩下的应变片数据作为煤矿巷道剪应力计算的数据源，保证在分析整个杆体剪应力分布时不会受到错误数据干扰，并对部分情况下错误数据进行合理修复，提高煤矿巷道剪应力采集的鲁棒性，同时提高杆体的使用寿命；

(3)本发明对于一对应变片中只有一个应变片发生故障的情况采用了一种误差模型补偿的故障修复策略，该误差模型是基于模拟煤矿巷道应力突变及模拟煤矿巷道应力平缓综合数据采集的情况下建立的，保证了该误差模型更加接近煤矿巷道剪应力检测杆体在煤矿综合工况下的误差真实模型，保证了故障修复的鲁棒性及准确性；

(4)本发明对被测杆体的弯矩进行三次样条插值拟合，插值函数节点处的波动只对节点两侧的分段有影响，离该点较远的分段影响逐渐变小，有较好的稳定性。同时，由于个别端点的测量值异常不会引起被测杆体整体数据的异常，提高精度。

附图说明

图1为本发明的工作原理流程图；

具体实施方法

一种故障自修复的鲁棒的煤矿巷道剪应力采集方法，其特征在于，包括以下

步骤：

步骤1：设煤矿巷道剪应力测试杆体长度L，在该杆体上贴N对应变片，每对应变片之间的间隔d＝L/(N-1)，每对应变片正对贴在杆体上，设上面一层应变片的编号为X_i1，正对面一层应变片编号为X_i2(i＝0、1……N-1，i为第i对应变片)，将该杆体安装在对其进行弯曲、拉伸的相关设备上；

步骤2：误差模型的输出值的评估：

步骤2.1：利用相关设备对装有2N个应变片的煤矿巷道剪应力测试杆体以时间间隔t₁分别进行C次拉伸应变和C次弯曲应变，采集并记录弯曲应变、拉伸应变数值，记为ε_i1 ^拉伸，ε_i2 ^拉伸，ε_i1 ^弯曲，ε_i2 ^弯曲(i＝0、1……N-1)，此处时间间隔t₁时间间隔较短，模拟的是煤矿巷道应变突变的情况下的数据采集；

步骤2.2：分别计算煤矿巷道剪应力测试杆体模拟煤矿巷道应变突变的情况下每一对应变片的拉伸应变和弯曲应变数值误差δ_ij ^拉伸＝ε_i1 ^拉伸-ε_i2 ^拉伸，δ_ij ^弯曲＝ε_i1 ^弯曲+ε_i2 ^弯曲(i＝0、1……N-1，j＝0、1……C-1)；

步骤2.3：利用相关设备对装有2N个应变片的煤矿巷道剪应力测试杆体以时间间隔t₂分别进行C′次拉伸应变和C′弯曲应变，采集并记录弯曲应变、拉伸应变数值，记为(i＝0、1……N-1)，此处时间间隔t₂时间间隔较长，模拟的是煤矿巷道应变变化缓慢的情况下的数据采集；

步骤2.4：分别计算设煤矿巷道剪应力测试杆体模拟煤矿巷道应变变化缓慢的情况下每一对应变片的拉伸应变和弯曲应变数值误差(i＝0、1……N-1，j′＝0、1……C′-1)；

步骤2.5：分别计算设煤矿巷道剪应力测试杆体每一对应变片拉伸应变和弯曲应变应的数值综合工况下误差平均值，(i＝0、1……N-1)，该平均值即为煤矿巷道综合工况下的应变误差模型输出值的估算值，该估算值即为应变偏差补偿的数据源；

步骤3：将剪应力测试杆体安装在煤矿巷道里，进行煤矿巷道综合工况下的数据采集。

步骤4：根据采集到的数据判断应变片是否故障：

步骤4.1：每隔时间t对杆体上的应变片采集数值进行读取；

步骤4.2：读取杆体上2N个应变片的采集数值ε_i1，ε_i2(i＝0、1……N-1)；

步骤4.3：如果连续三次ε_i1或者ε_i2值为0，则认为应变片发生故障，进入步骤5；

步骤4.4：如果ε_i1，ε_i2(i＝0、1……N-1)都不为0，即应变片未出现故障，则此时有效应变片对数N′＝N，进入步骤7；

步骤5：应变片故障类型判断：

步骤5.1：如果应变片X_i1出现故障，应变片X_i2正常工作或者应变片X_i2出现故障，应变片X_i1正常工作，进入步骤6.1；

步骤5.2：如果应变片X_i1和X_i2都出现故障，则记录应变片X_i1和X_i2都出现故障的应变片对数为P，进入步骤6.2；

步骤6：根据故障类型进入相应的应变片故障处理流程：

步骤6.1：对于故障应变片的采集数值进行误差补偿，如果X_i1出现故障，则计算拉伸应变时对X_i1应变片进行补偿，令则计算弯曲应变和弯矩时令则可以得到第i个应变片的弯矩：M_i ^弯矩＝π·E·R³·(ε′_i1-ε_i2)/8，其中π是圆周率，E为弹性系数，R为杆体的截面半径；如果X_i2出现故障，则计算拉伸应变时对X_i2进行补偿，令则计算弯曲应变和弯矩时令则可以得到第i个应变片的弯矩：M_i ^弯矩＝π·E·R³·(ε_i1-ε′_i2)/8，进入步骤8；

步骤6.2：如果第i对两个应变片都是有故障的，则剔除这一对应变片的采样数值，将剩下的正常工作的应变片作为有效测点，则修正后有效应变片对数N′＝N-P进入步骤7；

步骤7：被测杆体应变参数计算，E为弹性系数，R为杆体的截面半径，A为杆体横截面积，i′为故障处理之后有效应变片的序号，ε_i′1,ε_i′2为步骤6中进行误差补偿之后的应变值：

步骤7.1：计算被测杆体上每个应变片测点拉伸应变参数，拉伸应变为拉伸应力为(i′＝0、1……N′-1)；

步骤7.2：；计算被测杆体上每个应变片测点弯曲应变参数，弯曲应变为弯曲应力为(i′＝0、1……N′-1)；

步骤7.3：；计算被测杆体每个有效应变片点弯矩：

M_i′＝π·E·R³·(ε_i′1-ε_i′2)/8，(i′＝0、1……N′-1)    (1)

步骤8：对被测杆体的N′对有效应变片的弯矩M_i′进行函数拟合，在被测杆体由应变片分割的N′-1个区间段上使用三次样条插值拟合获取被测杆体上的弯矩函数，每个区间段上的拟合函数的一阶导数和二阶导数均为连续函数，则获得的拟合弯矩函数曲线是光滑曲线，具体拟合步骤如下：

步骤8.1：被测杆体的N′对有效应变片，相邻两个应变片之间是一个区间，所以有(N′-1)个区间，设每个区间段[z_k,z_k+1](k＝0,1,...,N′-2)，k表示(N′-1)个区间的序号。设第k个区间[z_k,z_k+1]的三次样条插值拟合函数为S_k(x)，S_k(x)表示距离被测杆体端点为x处的弯矩，即：

S_k(x)＝a_k+b_k(x-z_k)+c_k(x-z_k)²+d_k(x-z_k)³    (2)

其中x是距离被测杆端点的距离，z_k＝k×d(k＝0,1,...,N′-2)是当应变片i＝k序号时，该应变片距离被测杆端点的距离，a_k,b_k,c_k,d_k为三次样条插值拟合函数需要求解的参数。

根据三次样条插值拟合函数的要求，其插值条件为：

当x＝z_k时，S_k(z_k)＝M_k    (3)

其中，M_k是当i′＝k时M_i′的值，M_i′的值从公式(1)得到。

三次样条插值拟合函数的连续性条件为：

$\underset{x\→z_k}{\lim }{S}_k\left(x\right)=S_k\left(z_k\right)---\left(4\right)$

三次样条插值拟合函数的一阶导数连续条件为：

$\underset{x\→z_k}{\lim }{S_k}^{\′}\left(x\right)={S_k}^{\′}\left(z_k\right)---\left(5\right)$

其中，S_k′(x)是S_k(x)的一次导数。

三次样条插值拟合函数的二阶导数连续条件为：

$\underset{x\→z_{i^{\′}}}{\lim }{S_k}^{\′\′}\left(x\right)={S_k}^{\′\′}\left(z_k\right)---\left(6\right)$

其中，S_k″(x)是S_k(x)的二次导数。

以上插值条件和连续性条件总共能够得到(4N′-2)条等式。而三次样条插值在每一个区间上有4个未知参数，因此总共有4N′个待求解未知参数。被测杆体两端的端点位置没有受到让它们弯曲的力，因此端点条件为自由边界，所以边界端点约束提供了两个等式，为：

通过以上4N′条等式，能够求解出S_k(x)中4N′个未知参数a_k,b_k,c_k,d_k(k＝0、1……N′-2)。

步骤8.2：求解三次样条插值拟合函数为S_k(x)。

步骤8.2.1：将拟合函数S_k(x)的二阶导数作为待定参数，令

T_k＝S_k″(x)；

步骤8.2.2：因为S_k(x)是三次多项式，所以其二阶导数是一次多式，

令 $S_k^{\′\′}\left(x\right)=\frac{T_{k+1}-T_k}{h_k}(x-z_k)+T_k,x\∈{[z}_k,z_{k+1}]---\left(8\right)$

其中，h_k＝z_k+1-z_k。

对二阶导数积分，得到一阶导数：

$S_k^{\′}\left(x\right)=\frac{T_{k+1}-T_k}{{2h}_k}{(x-z_k)}^2+T_k(x-z_k)+p_k---\left(9\right)$

其中，p_k为常量参数。

对一阶导数积分，得到函数S_k(x)：

$S_k\left(x\right)=\frac{T_{k+1}-T_k}{{6h}_k}{(x-z_k)}^3+\frac{T_k}{2}{(x-z_k)}^2+p_k(x-z_k)+q_k---\left(10\right)$

其中，q_k为常量参数。

将x＝z_k带入公式(9)中，可以得到S_k(z_k)＝q_k，q_k即当x＝z_k时S_k(x)的值，根据公式(3)可知S_k(z_k)等于序号为k的应变片的弯矩值，即M_k，因此q_k是已知量。

因为拟合曲线平滑时端点处函数值相等，所以

S_k(z_k+1)＝S_k+1(z_k+1)＝q_k+1    (11)

将公式(11)代入公式(10)中，可以得到：

$S_k\left(z_{k+1}\right)=\frac{T_{k+1}-T_k}{{6h}_k}{(z_{k+1}-z_k)}^3+\frac{T_k}{2}{(z_{k+1}-z_k)}^2+p_k(z_{k+1}-z_k)+q_k=q_{k+1}---\left(12\right)$

由公式(12)得到

$p_k=f[z_k,z_{k+1}]-\frac{T_{k+1}+{2T}_k}{6}{h}_k---\left(13\right)$

其中， $f[z_k,z_{k+1}]=\frac{q_{k+1}-q_k}{z_{k+1}-z_k}---\left(14\right)$

将公式(13)带入公式(10)中可以得到三次样条插值拟合函数为：

$S_k\left(x\right)=\frac{T_{k+1}-T_k}{{6h}_k}{(x-z_k)}^3+\frac{T_k}{2}{(x-z_k)}^2+\left(f\right[z_k,z_{k+1}]-\frac{T_{k+1}+{2T}_k}{6}{h}_k)(x-z_k)+q_k---\left(15\right)$

公式(15)中，T_k,T_k+1是未知量，需要继续求解。

步骤8.2.3：确定二阶导数待定参数T_k。

因为拟合曲线一次导数平滑所以端点处相邻两个函数值相等：

S′_k-1(z_k)＝S′_k(z_k)    (16)

将公式(16)代入公式(9)中，即：

$\begin{array}{c}\frac{T_k-T_{k-1}}{{2h}_{k-1}}{(z_k-z_{k-1})}^2+T_k(z_k-z_{k-1})+f[z_{k-1},z_k]-\frac{T_k+{2T}_{k-1}}{6}{h}_{k-1}\\ =f[z_k,z_{k+1}]-\frac{T_{k+1}+{2T}_k}{6}{h}_k\end{array}---\left(17\right)$

化简公式(17)，可以得到：

$\frac{h_{k-1}}{h_{k-1}+h_k}{T}_{k-1}+{2T}_k+(1-\frac{h_{k-1}}{h_{k-1}+h_k})T_{k+1}=6\frac{f[z_k,z_{k+1}]-f[z_{k-1},z_k]}{h_{k-1}+h_k}---\left(18\right)$

令：

${\mathrm{\μ}}_k=\frac{h_{k-1}}{h_{k-1}+h_k},{\mathrm{\λ}}_k=1-{\mathrm{\μ}}_k,$

$d_k=\frac{f[z_k,z_{k+1}]-f[z_{k-1},z_k]}{h_{k-1}+h_k}(k=1......N^{\′}-2),d_0=f[z_0,z_1],d_{N^{\′}-1}=f[z_{N^{\′}-2},z_{N^{\′}-1}]$

则公式(18)可以化简为：

μ_kT_k-1+2T_k+λ_kT_k+1＝6d_k    (19)

又因为边界条件公式(7)可以得到

$\left\{\begin{array}{c}{S_0}^{\′\′}\left(z_0\right)=0\\ {S^{\′\′}}_{N^{\′}-1}\left(z_{N^{\′}-1}\right)=0\end{array}\right.$ 即 $\left\{\begin{array}{c}{T}_0=0\\ T_{N^{\′}-1}=0\end{array}\right.---\left(20\right)$

可以得到求解T_k(k＝1,2......N′-2)的矩阵方程，(N′-2)个方程可以求解出(N′-2)个T_k的值：

步骤8.3：将求得的T_k(k＝0,1......N′-1)代入三次样条插值拟合函数S_k(x)公式(15)中，获得每一个区间段[z_k,z_k+1](k＝0,1,...,N′-2)对应的弯矩函数：

$S_k\left(x\right)=\frac{T_{k+1}-T_k}{{6h}_k}{(x-z_k)}^3+\frac{T_k}{2}{(x-z_k)}^2+\left(f\right[z_k,z_{k+1}]-\frac{T_{k+1}+{2T}_k}{6}{h}_k)(x-z_k)+q_k---\left(22\right)$

步骤9：沿测点长度范围剪应力分布函数：

步骤9.1：对公式(22)的函数S_k(x)进行求导运算得到S′_k(x)，根据拟合的弯矩函数的导数S′_k(x)得到每个区间段剪应力分布函数：

τ_k(x)＝4·S′_k(x)/(3π·R²)    (23)

其中x是距离被测杆端点的距离，R是杆体横截面的半径；

步骤9.2：将被测杆体上任意一点距离杆体端点的距离z代入公式(23)中，可以得到该点的剪应力。

以上的实施方法只是已实现的有效的具体实施方式之一，本领域的技术人员在本发明技术方案范围内进行的通常变化和替换都应包含在本发明的保护范围内。

Claims (1)

1.一种故障自修复的鲁棒的煤矿巷道剪应力采集方法，其特征在于，包括以下步骤：

步骤1：设煤矿巷道剪应力测试杆体长度L，在该杆体上贴N对应变片，每对应变片之间的间隔d＝L/(N-1)，每对应变片正对贴在杆体上，设上面一层应变片的编号为X_i1，正对面一层应变片编号为X_i2(i＝0、1……N-1，i为第i对应变片)，将该杆体安装在对其进行弯曲、拉伸的相关设备上；

步骤2：误差模型的输出值的评估：

步骤2.1：利用相关设备对装有2N个应变片的煤矿巷道剪应力测试杆体以时间间隔t₁分别进行C次拉伸应变和C次弯曲应变，采集并记录弯曲应变、拉伸应变数值，记为ε_i1 ^拉伸，ε_i2 ^拉伸，ε_i1 ^弯曲，ε_i2 ^弯曲(i＝0、1……N-1)，此处时间间隔t₁时间间隔较短，模拟的是煤矿巷道应变突变的情况下的数据采集；

步骤2.2：分别计算煤矿巷道剪应力测试杆体模拟煤矿巷道应变突变的情况下每一对应变片的拉伸应变和弯曲应变数值误差δ_ij ^拉伸＝ε_i1 ^拉伸-ε_i2 ^拉伸，δ_ij ^弯曲＝ε_i1 ^弯曲+ε_i2 ^弯曲(i＝0、1……N-1，j＝0、1……C-1)；

步骤2.3：利用相关设备对装有2N个应变片的煤矿巷道剪应力测试杆体以时间间隔t₂分别进行C′次拉伸应变和C′弯曲应变，采集并记录弯曲应变、拉伸应变数值，记为ε_i1 ^′拉伸，ε_i2 ^′拉伸，ε_i1 ^′弯曲，ε_i2 ^′弯曲(i＝0、1……N-1)，此处时间间隔t₂时间间隔较长，模拟的是煤矿巷道应变变化缓慢的情况下的数据采集；

步骤2.4：分别计算设煤矿巷道剪应力测试杆体模拟煤矿巷道应变变化缓慢的情况下每一对应变片的拉伸应变和弯曲应变数值误差δ_ij ^′拉伸＝ε_i1 ^′拉伸-ε_i2 ^′拉伸，δ_ij ^′弯曲＝ε_i1 ^′弯曲+ε_i2 ^′弯曲(i＝0、1……N-1，j′＝0、1……C′-1)；

步骤2.5：分别计算设煤矿巷道剪应力测试杆体每一对应变片拉伸应变和弯曲应变应的数值综合工况下误差平均值， 该平均值即为煤矿巷道综合工况下的应变误差模型输出值的估算值，该估算值即为应变偏差补偿的数据源；

步骤3：将剪应力测试杆体安装在煤矿巷道里，进行煤矿巷道综合工况下的数据采集。

步骤4：根据采集到的数据判断应变片是否故障：

步骤4.1：每隔时间t对杆体上的应变片采集数值进行读取；

步骤4.2：读取杆体上2N个应变片的采集数值ε_i1，ε_i2(i＝0、1……N-1)；

步骤4.3：如果连续三次ε_i1或者ε_i2值为0，则认为应变片发生故障，进入步骤5；

步骤4.4：如果ε_i1，ε_i2(i＝0、1……N-1)都不为0，即应变片未出现故障，则此时有效应变片对数N′＝N，进入步骤7；

步骤5：应变片故障类型判断：

步骤5.1：如果应变片X_i1出现故障，应变片X_i2正常工作或者应变片X_i2出现故障，应变片X_i1正常工作，进入步骤6.1；

步骤5.2：如果应变片X_i1和X_i2都出现故障，则记录应变片X_i1和X_i2都出现故障的应变片对数为P，进入步骤6.2；

步骤6：根据故障类型进入相应的应变片故障处理流程：

步骤6.1：对于故障应变片的采集数值进行误差补偿，如果X_i1出现故障，则计算拉伸应变时对X_i1应变片进行补偿，令则计算弯曲应变和弯矩时令则可以得到第i个应变片的弯矩：M_i ^弯矩＝π·E·R³·(ε′_i1-ε_i2)/8，其中π是圆周率，E为弹性系数，R为杆体的截面半径；如果X_i2出现故障，则计算拉伸应变时对X_i2进行补偿，令则计算弯曲应变和弯矩时令则可以得到第i个应变片的弯矩：M_i ^弯矩＝π·E·R³·(ε_i1-ε′_i2)/8，进入步骤8；

步骤6.2：如果第i对两个应变片都是有故障的，则剔除这一对应变片的采样数值，将剩下的正常工作的应变片作为有效测点，则修正后有效应变片对数N′＝N-P进入步骤7；

步骤7：被测杆体应变参数计算，i′为故障处理之后有效应变片的序号，ε_i′1,ε_i′2为步骤6中进行误差补偿之后的应变值：

步骤7.1：计算被测杆体上每个应变片测点拉伸应变参数，拉伸应变为拉伸应力为

步骤7.2：；计算被测杆体上每个应变片测点弯曲应变参数，弯曲应变为弯曲应力为

步骤7.3：；计算被测杆体每个有效应变片点弯矩：

M_i′＝π·E·R³·(ε_i′1-ε_i′2)/8，(i′＝0、1……N′-1)    (1)

步骤8：对被测杆体的N′对有效应变片的弯矩M_i′进行函数拟合，在被测杆体由应变片分割的N′-1个区间段上使用三次样条插值拟合获取被测杆体上的弯矩函数，每个区间段上的拟合函数的一阶导数和二阶导数均为连续函数，则获得的拟合弯矩函数曲线是光滑曲线，具体拟合步骤如下：

步骤8.1：被测杆体的N′对有效应变片，相邻两个应变片之间是一个区间，所以有(N′-1)个区间，设每个区间段[z_k,z_k+1](k＝0,1,...,N′-2)，k表示(N′-1)个区间的序号。设第k个区间[z_k,z_k+1]的三次样条插值拟合函数为S_k(x)，S_k(x)表示距离被测杆体端点为x处的弯矩，即：

S_k(x)＝a_k+b_k(x-z_k)+c_k(x-z_k)²+d_k(x-z_k)³   (2)

其中x是距离被测杆端点的距离，z_k＝k×d(k＝0,1,...,N′-2)是当应变片i＝k序号时，该应变片距离被测杆端点的距离，z_k+1＝(k+1)×d(k＝0,1,...,N′-2)是当应变片i＝k+1序号时，该应变片距离被测杆端点的距离，a_k,b_k,c_k,d_k为三次样条插值拟合函数需要求解的参数。

根据三次样条插值拟合函数的要求，其插值条件为：

当x＝z_k时，S_k(z_k)＝M_k    (3)

其中，M_k是当i′＝k时M_i′的值，M_i′的值从公式(1)得到。

三次样条插值拟合函数的连续性条件为：

$\underset{x\→z_k}{\lim }{S}_k\left(x\right)=S_k\left(z_k\right)---\left(4\right)$

三次样条插值拟合函数的一阶导数连续条件为：

$\underset{x\→z_k}{\lim }{S_k}^{\′}\left(x\right)={S_k}^{\′}\left(z_k\right)---\left(5\right)$

其中，S_k'(x)是S_k(x)的一次导数。

三次样条插值拟合函数的二阶导数连续条件为：

$\underset{x\→z_k}{\lim }{S_k}^{\′\′}\left(x\right)={S_k}^{\′\′}\left(z_k\right)---\left(6\right)$

其中，S_k″(x)是S_k(x)的二次导数。

以上插值条件和连续性条件总共能够得到(4N′-2)条等式。而三次样条插值在每一个区间上有4个未知参数，因此总共有4N′个待求解未知参数。被测杆体两端的端点位置没有受到让它们弯曲的力，因此端点条件为自由边界，所以边界端点约束提供了两个等式，为：

其中，z₀＝0×d＝0表示当应变片i＝0序号时，该应变片距离被测杆端点的距离为0，z_N′-1＝(N′-1)×d是当应变片i＝N′-1序号时，该应变片距离被测杆端点的距离，通过以上4N′条等式，能够求解出S_k(x)中4N′个未知参数a_k,b_k,c_k,d_k(k＝0、1……N′-2)。

步骤8.2：求解三次样条插值拟合函数为S_k(x)。

步骤8.2.1：将拟合函数S_k(x)的二阶导数作为待定参数，令T_k＝S_k″(x)；

步骤8.2.2：因为S_k(x)是三次多项式，所以其二阶导数是一次多式，

令 $\begin{array}{cc}{S}_k^{\′\′}\left(x\right)=\frac{T_{k+1}-T_k}{h_k}(x-z_k)+T_k& x\∈[z_k,z_{k+1}]---\left(8\right)\end{array}$

其中，h_k＝z_k+1-z_k。

对二阶导数积分，得到一阶导数：

$S_k^{\′}\left(x\right)=\frac{T_{k+1}-T_k}{{2h}_k}{(x-z_k)}^2+T_k(x-z_k)+p_k---\left(9\right)$

其中，p_k为常量参数。

对一阶导数积分，得到函数S_k(x)：

$S_k\left(x\right)=\frac{T_{k+1}-T_k}{{6h}_k}{(x-z_k)}^3+\frac{T_k}{2}{\left({x-z}_k\right)}^2+p_k(x-z_k)+q_k---\left(10\right)$

其中，q_k为常量参数。

将x＝z_k带入公式(9)中，可以得到S_k(z_k)＝q_k，q_k即当x＝z_k时S_k(x)的值，根据公式(3)可知S_k(z_k)等于序号为k的应变片的弯矩值，即M_k，因此q_k是已知量。

因为拟合曲线平滑时端点处函数值相等，所以

S_k(z_k+1)＝S_k+1(z_k+1)＝q_k+1   (11)

将公式(11)代入公式(10)中，可以得到：

$S_k\left(z_{k+1}\right)=\frac{T_{k+1}-T_k}{{6h}_k}{(z_{k+1}-z_k)}^3+\frac{T_k}{2}{(z_{k+1}-z_k)}^2+p_k(z_{k+1}-z_k)+q_k=q_{k+1}---\left(12\right)$

由公式(12)得到

$p_k=f[z_k,z_{k+1}]-\frac{T_{k+1}+{2T}_k}{6}{h}_k---\left(13\right)$

其中， $f[z_k,z_{k+1}]=\frac{q_{k+1}-q_k}{z_{k+1}-z_k}---\left(14\right)$

将公式(13)带入公式(10)中可以得到三次样条插值拟合函数为：

$S_k\left(x\right)=\frac{T_{k+1}-T_k}{{6h}_k}{(x-z_k)}^3+\frac{T_k}{2}{(x-z_k)}^2+\left(f\right[z_k,z_{k+1}]-\frac{T_{k+1}+{2T}_k}{6}{h}_k)(x-z_k)+q_k---\left(15\right)$

公式(15)中，T_k,T_k+1是未知量，需要继续求解。

步骤8.2.3：确定二阶导数待定参数T_k。

因为拟合曲线一次导数平滑所以端点处相邻两个函数值相等：

S′_k-1(z_k)＝S′_k(z_k)         (16)

将公式(16)代入公式(9)中，即：

$\begin{array}{c}\frac{T_k-T_{k-1}}{{2h}_{k-1}}{(z_k-z_{k-1})}^2+T_k(z_k-z_{k-1})+f[z_{k-1},z_k]-\frac{T_k+{2T}_{k-1}}{6}{h}_{k-1}\\ =f[z_k,z_{k+1}]-\frac{T_{k+1}+{2T}_k}{6}{h}_k\end{array}---\left(17\right)$

化简公式(17)，可以得到：

$\frac{h_{k-1}}{h_{k-1}+h_k}{T}_{k-1}+{2T}_k+(1-\frac{h_{k-1}}{h_{k-1}+h_k})T_{k+1}=6\frac{f[z_k,z_{k+1}]-f[z_{k-1},z_k]}{h_{k-1}+h_k}---\left(18\right)$

令：

$\begin{array}{c}{\mathrm{\μ}}_k=\frac{h_{k-1}}{h_{k-1}+h_k},{\mathrm{\λ}}_k=1-{\mathrm{\μ}}_k,\\ d_k=\frac{f[z_k,z_{k+1}]-f[z_{k-1},z_k]}{h_{k-1}+h_k}(k=1......N^{\′}-2),d_0=f[z_0,z_1],d_{N^{\′}-1}=f[z_{N^{\′}-2},z_{N^{\′}-1}]\end{array}$

则公式(18)可以化简为：

μ_kT_k-1+2T_k+λ_kT_k+1＝6d_k(19)

又因为边界条件公式(7)可以得到

$\left\{\begin{array}{c}{S_0}^{\′\′}\left(z_0\right)=0\\ {S^{\′\′}}_{N^{\′}-1}\left(Z_{N^{\′}-1}\right)=0\end{array}\right.$ 即 $\left\{\begin{array}{c}{T}_0=0\\ T_{N^{\′}-1}\end{array}\right.---\left(20\right)$

可以得到求解T_k(k＝1,2......N′-2)的矩阵方程，(N′-2)个方程可以求解出(N′-2)个T_k的值：

步骤8.3：将求得的T_k(k＝0,1......N′-1)代入三次样条插值拟合函数S_k(x)公式(15)中，获得每一个区间段[z_k,z_k+1](k＝0,1,...,N′-2)对应的弯矩函数：

$S_k\left(x\right)=\frac{T_{k+1}-T_k}{{6h}_k}{(x-z_k)}^3+\frac{T_k}{2}{(x-z_k)}^2+\left(f\right[z_k,z_{k+1}]-\frac{T_{k+1}+{2T}_k}{6}{h}_k)(x-z_k)+q_k---\left(22\right)$

步骤9：沿测点长度范围剪应力分布函数：

步骤9.1：对公式(22)的函数S_k(x)进行求导运算得到S′_k(k)，根据拟合的弯矩函数的导数S′_k(x)得到每个区间段剪应力分布函数：

τ_k(x)＝4·S′_k(x)/(3π·R²)     (23)

其中x是距离被测杆端点的距离，R是杆体横截面的半径；

步骤9.2：将被测杆体上任意一点距离杆体端点的距离z代入公式(23)中，可以得到该点的剪应力。