一种横向滑坡埋地管道工况参数智能辨识方法

技术领域

本发明涉及工程参数反分析领域，具体涉及一种横向滑坡埋地管道工况参数智能辨识方法。

背景技术

准确获取横向滑坡作用下埋地管道的工况参数是通过力学分析方法进行管道安全预警的重要基础。传统的通过现场试验测定工况参数的方式，效率低下，实用性较低。随着横向滑坡长输油气管道安全事故的频发，管道安全监测工作不断开展，从而积累了大量的监测信息，利用管道的监测信息辨识管道工况参数成为一种技术趋势。

近些年，智能参数辨识方法得到了广泛的研究。这些辨识技术的区别主要体现在目标函数、是否采用及采用何种近似方法代替数值计算、参数优化算法等方面。在目标函数方面，主要包括误差平方和方法与概率方法；在计算方面，主要包括商业软件模拟计算和机器学习方法，主要使用到的机器学习方法有支持向量机和神经网络等；在优化算法方面，主要包括最速梯度下降、粒子群算法、模拟退火算法、遗传算法以及各种改进算法等。这些辨识技术在利用静态监测信息进行参数辨识方面取得了一定的进展，但不能充分利用连续的时间序列监测信息，参数辨识结果可靠性较低。

发明内容

本发明的目的在于提供一种横向滑坡埋地管道工况参数智能辨识方法，以解决现有辨识技术因不能利用连续监测信息导致参数识别结果可靠性较低的技术问题。

为解决上述技术问题，本发明具体提供下述技术方案：

一种横向滑坡埋地管道工况参数智能辨识方法，包括如下步骤：

S100、通过大量的模板工况及各个模板工况下管道横截面的轴向最大应变数据建立模板数据库；

S200、监测实际工况中管道各监测截面的轴向最大应变，形成监测应变时间序列；

S300、利用多维动态时间规整算法，计算模板数据库和监测应变时间序列的相似度，通过比较相似度大小辨识实际工况参数；

其中，模板工况越多，模板数据库越大，辨识结果越准确。

作为本发明的一种优选方案，步骤S100具体包括；

S101、估算实际工况各参数的取值范围；

S102、将各参数在各自的取值范围内进行离散，并组合离散参数以形成多个模板工况；

S103、利用力学模型和遗传算法计算各模板工况下管道横截面的轴向最大应变，形成模板数据库；

作为本发明的一种优选方案，步骤S101中的参数具体包括滑坡中心位置d'、滑坡推力q的变化范围、滑坡宽度L、以及管下地基土体抗力系数k。

其中，d'的取值范围为[d

q的变化范围无需估计；

L的取值范围[L

k的取值范围[k

其中，d'为以监测点分布中心为坐标原点、以管道轴向为坐标轴x'时，滑坡中心的坐标。

作为本发明的一种优选方案，步骤S102具体包括：

在参数取值范围内，将滑坡宽度L离散为L

将管下地基土体抗力系数k离散为k

记M＝ab，将离散的L和k组成M个模板工况，即(L

其中，a和b分别为L和k的离散个数。

作为本发明的一种优选方案，步骤S103具体包括：

S1031、对于各模板工况，按一定的载荷增量增大滑坡推力，直至管道最大的轴向拉应力大于管道的屈服强度，从而管道发生破坏；

S1032、对于每一载荷步，选择相应的力学模型计算管道横截面的轴向最大应变，得到相关的方程组；

S1033、通过遗传算法求解上述力学模型的方程组，最终得到各模板工况所有载荷步的管道横截面的轴向最大应变。

作为本发明的一种优选方案，步骤S1032中的力学模型具体包括全弹性、塑性区小于滑坡区、以及塑性区大于滑坡区三种模型；

其中，管土相互作用采用理想弹塑性本构关系，全弹性表示整个管道的管土之间均处于弹性相互作用，塑性区小于滑坡区表示管土之间的塑性相互作用区域没有超出滑坡区，塑性区大于滑坡区表示管土之间的塑性相互作用区域超出滑坡区；

根据管道挠曲微分方程可以求得三种力学模型下管道各段的挠度通解表达式，结合管道的边连条件、临界条件和伸长量条件，可以得到相应的待求方程组；

然后，由于各模型关于滑坡中心对称，取模型右半部分为研究对象，将滑坡区内的管道称为第一和第二段管道，将滑坡区外的管道称为第三段管道；

对于全弹性模型，管道各段的挠曲微分方程、挠度通解和定解条件为：

当0≤x≤L/2时，管道第一段挠曲微分方程为

EIω

挠度通解为

当L/2≤x<∞时，管道第三段挠曲微分方程为

EIω

挠度通解为

ω

管道的定解条件为

边连条件

伸长量条件

记作

对于塑性区小于滑坡区模型，管道各段的挠曲微分方程、挠度通解和定解

条件为：

当0≤x≤Δ时，管道第一段挠曲微分方程为

EIω

挠度通解为

当Δ≤x≤L/2时，管道第二段挠曲微分方程为

EIω

挠度通解为

当L/2≤x<∞时，管道第三段挠曲微分方程为

EIω

挠度通解为

ω

管道的定解条件为

边连条件

伸长量条件

记作

对于塑性区大于滑坡区模型，管道各段的挠曲微分方程、挠度通解和定解条件

当0≤x≤L/2时，管道第一段挠曲微分方程为

EIω

挠度通解为

当L/2≤x≤Δ时，管道第二段挠曲微分方程为

EIω

挠度通解为

当L/2≤x<∞时，管道第三段挠曲微分方程为

EIω

挠度通解为

ω

管道的定解条件为

边连条件

临界条件和伸长量条件

记作

其中，w

另外，c

作为本发明的一种优选方案，步骤S1033的遗传算法为：

根据塑性区半宽度Δ和管道轴力T进行迭代，迭代过程中，通过边连条件计算方程组中剩余的未知数，直至收敛。

作为本发明的一种优选方案，在步骤S1033中的遗传算法分别求解三种力学模型的具体操作步骤为：

S10331、对于全弹性模型，随机初始化T对应的种群，种群个体为T

S10332、通过边连条件计算剩余未知数，然后计算父代种群各个体的适应度，对各个体进行选择、交叉和变异三步操作，得到子代种群；

对于全弹性模型，适应度计算公式为

对于塑性区小于滑坡区模型和塑性区大于滑坡区模型，适应度计算公式为

S10333、通过边连条件计算剩余未知数，然后计算子代种群各个体的适应度，记录适应度最大的个体数值；

S10334、若S10333中最大适应度小于设定的阈值，则以子代为父代，重复S10331～S10333，直至最大适应度大于设定的阈值，记录此适应度对应的个体，该个体转换为十进制后即为Δ和T的最优解，算法终止。

S10335、通过遗传算法求得方程组的未知数后，将各未知数的值代入管道的挠度通解表达式，根据挠度与横截面的最大轴向应变的关系

可以求得管道横截面的轴向最大应变通解表达式；

S10336、在模板工况中，以滑坡中心为坐标原点，以管道轴向为坐标轴x，将整个管道间隔相同距离，取N个管道横截面的轴向最大应变，将其称作离散点应变时间序列，最终求得M个模板工况的离散点应变时间序列Y，具体表示为

Y＝[Y

其中，Y

作为本发明的一种优选方案，步骤S200中时间序列形成的具体操作为：

在实际工况中，通过检测管道的各监测截面的轴向最大应变，形成监测应变时间序列X，具体表示为

其中，p为最大时刻，n为监测截面的个数，x

作为本发明的一种优选方案，步骤S300的具体步骤包括：

S301、设定移动时窗用来从每个模板工况离散点应变时间序列Y

W为t×n的矩阵；

其中，t为移动时窗所含载荷步数，即时窗长度，其变化范围为[t

S302、设模板工况滑坡中心在坐标轴x'的坐标为d，d的变化范围为[d

x＝x′-d

计算各监测截面在模板工况坐标轴x中的坐标，将这些坐标与离散点的x坐标进行比较，获得与这些坐标最近的离散点；

S303、开始时，计算d

S304、移动时窗每到达一个位置，即得到一个二维矩阵W，下面对移动时窗当前位置的W与X进行动态时间规整，定义距离d

遍历W和X里的各元素，即可得到距离矩阵

对于矩阵d，以d

S305、定义累计距离矩阵

其中Dij为d11为到dij的最小累计距离，计算公式为

到此可知，当前W和X的最小累计距离为D

S306、随着移动窗口在Yi中移动，可得相似度矩阵C

当t在[t

S307、对第i个模板工况的相似度矩阵，比较矩阵中各相似度的大小，求得该模板工况对应的最大相似度，以及相应的滑坡中心位置d，滑坡推力q变化范围；

S308、对M个模板工况S

S309、比较每个模板工况的最大相似度，可以得到所有模板工况下的模板数据库里与监测应变时间序列X最相似的序列片段，则该序列片段对应的工况参数即为辨识出的与实际工况最接近的工况参数。

本发明与现有技术相比较具有如下有益效果：

本发明采用遗传算法求解横向滑坡作用下埋地管道的力学模型，并生成模板数据库，采用多维动态时间规整算法，充分利用实际工况的连续时空监测信息，准确辨识实际工况参数，避免现场测试耗费大量的人力物力。

附图说明

为了更清楚地说明本发明的实施方式或现有技术中的技术方案，下面将对实施方式或现有技术描述中所需要使用的附图作简单地介绍。显而易见地，下面描述中的附图仅仅是示例性的，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图引申获得其它的实施附图。

图1为本发明实施例中的流程示意图；

图2是本发明实施例中力学模型的全弹性模型概化图；

图3是本发明实施例中力学模型的塑性区小于滑坡区模型概化图；

图4是本发明实施例中力学模型的塑性区大于滑坡区模型概化图；

图5是本发明实施例中管土作用本构关系示意图；

图6是本发明实施例中力学模型选择流程图；

图7是本发明实施例中遗传算法流程图；

图8是本发明实施例中多维动态时间规整算法示意图。

图9是本发明实施例中工况参数辨识模型示意图。

图10是本发明实施例中相似度计算方式(时窗长度45)示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，本发明提供了一种横向滑坡埋地管道工况参数智能辨识方法，包括如下步骤：

S100、通过大量的模板工况及各个模板工况下管道横截面的轴向最大应变数据建立模板数据库；

S200、监测实际工况中管道各监测截面的轴向最大应变，形成监测应变时间序列；

S300、利用多维动态时间规整算法，计算模板数据库和监测应变时间序列的相似度，通过比较相似度大小辨识实际工况参数；

其中，模板工况越多，模板数据库越大，辨识结果越准确。

需要特别说明的是，埋地管道采用X70管材，弹性模量E＝210GPa，外径D＝1.016m，壁厚t＝0.0175m，截面惯性矩I＝6.84×10

对于横向滑坡作用下埋地管道工程而言，需要识别的实际工况参数包括滑坡变形体中心位置d’，管道监测应变时间序列对应的滑坡横向推力q的变化范围，滑坡区宽度L，地基土体抗力系数k。具体识别流程如下：

(1)估计实际工况各参数的取值范围，在该取值范围内设计各模板工况的参数组合。本例对滑坡宽度L和地基土体抗力系数k分别取三个水平，形成9种模板工况。为了对辨识结果进行对比验证，同时也设定了实际工况参数。各工况参数如表1所示。

(2)利用力学模型和遗传算法计算各模板工况下管道横截面的轴向最大应变，形成模板数据库。对于每个模板工况，分100步按固定增量增加滑坡推力，直至管道发生破坏，计算每一步的管道横截面的轴向最大应变；在整个管道轴向均匀设置离散点，各模板工况离散点处横截面的轴向最大应变组成模板数据库。

(3)获得监测应变时间序列。本实施例为说明辨识方法有效性，把力学模型和遗传算法计算的实际工况下管道横截面的轴向最大应变作为监测应变时间序列。在实际工况中，同样分100步按固定增量增加滑坡推力，直至管道发生破坏；在管道上设置监测点，取后50步载荷步下所有监测点的应变组成监测应变时间序列。

(4)计算相似度，辨识实际工况参数。采用多维动态时间规整算法，计算模板数据库与监测应变时间序列的相似度，相似度最大的模板工况参数即为辨识出的实际工况参数。

其中，力学模型用来获得埋地管道横截面的轴向最大应变通解和相应的定解条件，遗传算法用来求解通解里的未知量，最终得到埋地管道横截面的轴向最大应变的特解。

图2-4为横向滑坡作用下埋地管道的模型概化图，分别是全弹性模型，塑性区小于滑坡区模型，塑性区大于滑坡区模型。图5是力学模型里用到的管土作用本构关系。对于某一工况，在某一滑坡推力下，管道的应变具体应采用何种力学模型进行计算，图6给出了模型选择流程。

首先进行全弹性模型计算，若计算结果显示管道的滑坡中心处的挠度w0小于临界wc，则管道为全弹性模型，计算结果可用；若计算结果显示管道的滑坡中心处的挠度w0大于临界wc，则进行塑性区小于滑坡区模型计算，若计算结果显示塑性区小于滑坡区，则管道为塑性区小于滑坡区模型，计算结果可用；若计算结果显示塑性区大于滑坡区，则进行塑性区大于滑坡区模型计算。其中wc＝0.1D。

图7是遗传算法流程图，用来求解力学模型的方程组，具体算法流程见本发明的发明内容；图8是多维动态时间规整算法示意图，用来计算模板数据库和监测应变时间序列的相似度，具体计算方法见本发明的发明内容。

工况设计如表1所示，R为设定的实际工况，S1～S9为模板工况，模板数据库和监测应变时间序列均通过力学模型和遗传算法计算得到。实际工况的滑坡推力分50个载荷步，从1.5351×10

表1

如图9所示，实际工况中，在管道轴向设置11个监测点，建立以监测点分布中心为坐标原点的坐标系x'o'y'，各监测点x'方向坐标为-20m，-16m，-12m，-8m，-4m，0m，4m，8m，12m，16m，20m，滑坡变形体中心x'方向坐标为d'＝2.5m。模板工况中，建立以滑坡变形体中心为坐标原点的坐标系xoy，在整个管道轴向均匀设置离散点，离散点间距为0.2m，离散点总数为501个。滑坡变形体中心o在x'o'y'坐标系下的坐标为d，d

x＝x′-d

d以0.2m增量从d

实际工况有50个载荷步，每个载荷步有11个监测点的应变，监测应变时间序列是50×11的二维矩阵X；每个模板工况有100个载荷步，每个载荷步有501个离散点的应变，则每个模板工况的离散点应变时间序列是100×501的二维矩阵Yi；对于某个d，从Yi取出与11个监测点对应的11个离散点的应变时间序列，这11个离散点的应变时间序列是100×11的二维矩阵。

定义移动时窗用来在d取不同的值时，从离散点应变时间序列Yi里以不断移动的方式取出序列片段W。由于实际工况有50个载荷步，故设定移动时窗的长度变化范围为[t1，t2]＝[40，60]，宽度为11，以二维矩阵W表示。图7为具体的相似度计算方式示意图，移动时窗长度为45，故W是45×11的二维矩阵，每次向前移动一步后，都将移动时窗里的序列片段W与监测点应变时间序列X进行多维动态时间规整，计算两者的相似度。

在整个相似度计算过程中，以模板数据库为遍历范围，以移动时窗内序列片段W与监测点应变时间序列X之间的多维动态时间规整为计算单位，进行实际工况参数辨识。

本例辨识结果如表2所示

表2

模板工况S5与监测应变时间序列的相似度最大，为77.0，则本例的工况参数最终辨识结果为滑坡变形体中心位置为2.6m，管道监测应变时间序列对应的滑坡横向推力变化范围1.6550×10

作为本发明的一种优选方案，步骤S100具体包括；

S101、估算实际工况各参数的取值范围；

S102、将各参数在各自的取值范围内进行离散，并组合离散参数以形成多个模板工况；

S103、利用力学模型和遗传算法计算各模板工况下管道横截面的轴向最大应变，形成模板数据库；；

作为本发明的一种优选方案，步骤S101中的参数具体包括滑坡中心位置d'、滑坡推力q的变化范围、滑坡宽度L、以及管下地基土体抗力系数k。；

其中，d'的取值范围为[d

q的变化范围无需估计；

L的取值范围[L

k的取值范围[k

其中，d'为以监测点分布中心为坐标原点、以管道轴向为坐标轴x'时，滑坡中心的坐标。

作为本发明的一种优选方案，步骤S102具体包括：

在参数取值范围内，将滑坡宽度L离散为L

将管下地基土体抗力系数k离散为k

记M＝ab，将离散的L和k组成M个模板工况，即(L

其中，a和b分别为L和k的离散个数。

作为本发明的一种优选方案，步骤S103具体包括：

S1031、对于各模板工况，按一定的载荷增量增大滑坡推力，直至管道最大的轴向拉应力大于管道的屈服强度，从而管道发生破坏；

S1032、对于各模板工况，按一定的载荷增量增大滑坡推力，直至管道最大的轴向拉应力大于管道的抗拉强度，从而管道发生破坏；

S1032、对于每一载荷步，选择相应的力学模型用以计算管道横截面的轴向最大应变，得到相关的方程组；

S1033、通过遗传算法求解上述力学模型的方程组，最终得到各模板工况所有载荷步的管道横截面的轴向最大应变。

作为本发明的一种优选方案，步骤S1033中的力学模型具体包括全弹性、塑性区小于滑坡区、以及塑性区大于滑坡区三种模型；

其中，管土相互作用采用理想弹塑性本构关系，全弹性表示整个管道的管土之间均处于弹性相互作用，塑性区小于滑坡区表示管土之间的塑性相互作用区域没有超出滑坡区，塑性区大于滑坡区表示管土之间的塑性相互作用区域超出滑坡区；

根据管道挠曲微分方程可以求得三种力学模型下管道各段的挠度通解表达式，结合管道的边连条件、临界条件和伸长量条件，可以得到相应的待求方程组；

然后，由于各模型关于滑坡中心对称，取模型右半部分为研究对象，将滑坡区内的管道称为第一和第二段管道(全弹性模型只有第一段)，将滑坡区外的管道称为第三段管道；

对于全弹性模型，管道各段的挠曲微分方程、挠度通解和定解条件为：

当0≤x≤L/2时，管道第一段挠曲微分方程为

EIω

挠度通解为

当L/2≤x<∞时，管道第三段挠曲微分方程为

EIω

挠度通解为

ω

管道的定解条件为

边连条件

伸长量条件

记作

对于塑性区小于滑坡区模型，管道各段的挠曲微分方程、挠度通解和定解

条件为：

当0≤x≤Δ时，管道第一段挠曲微分方程为

EIω

挠度通解为

当Δ≤x≤L/2时，管道第二段挠曲微分方程为

EIω

挠度通解为

当L/2≤x<∞时，管道第三段挠曲微分方程为

EIω

挠度通解为

ω

管道的定解条件为

边连条件

伸长量条件

记作

对于塑性区大于滑坡区模型，管道各段的挠曲微分方程、挠度通解和定解

条件为：

当0≤x≤L/2时，管道第一段挠曲微分方程为

EIω

挠度通解为

当L/2≤x≤Δ时，管道第二段挠曲微分方程为

EIω

挠度通解为

当L/2≤x<∞时，管道第三段挠曲微分方程为

EIω

挠度通解为

ω

管道的定解条件为

边连条件

临界条件和伸长量条件

记作

其中，w

另外，c

作为本发明的一种优选方案，步骤S1033的遗传算法为：

根据塑性区半宽度Δ和管道轴力T进行迭代，迭代过程中，通过边连条件计算方程组中剩余的未知数，直至收敛。

作为本发明的一种优选方案，在步骤S1033中的遗传算法分别求解三种力学模型的具体操作步骤为：

S10331、对于全弹性模型，随机初始化T对应的种群，种群个体为T

S10332、通过边连条件计算剩余未知数，然后计算父代种群各个体的适应度，对各个体进行选择、交叉和变异三步操作，得到子代种群；

对于全弹性模型，适应度计算公式为

对于塑性区小于滑坡区模型和塑性区大于滑坡区模型，适应度计算公式为

S10333、通过边连条件计算剩余未知数，然后计算子代种群各个体的适应度，记录适应度最大的个体数值；

S10334、(f)若(c)中最大适应度小于设定的阈值，则以子代为父代，重复(a)～(c)，直至最大适应度大于设定的阈值，记录此适应度对应的个体，该个体转换为十进制后即为Δ和T的最优解，算法终止。

S10335、通过遗传算法求得方程组的未知数后，将各未知数的值代入管道的挠度通解表达式，根据挠度与横截面的最大轴向应变的关系

可以求得管道横截面的轴向最大应变通解表达式；

S10336、在模板工况中，以滑坡中心为坐标原点，以管道轴向为坐标轴x，将整个管道间隔相同距离，取N个管道横截面的轴向最大应变，将其称作离散点应变时间序列，最终求得M个模板工况的离散点应变时间序列Y，具体表示为

Y＝[Y

其中，Y

作为本发明的一种优选方案，步骤S200中时间序列形成的具体操作为：

在实际工况中，通过检测管道的各监测截面的轴向最大应变，形成监测应变时间序列X，具体表示为

其中，p为最大时刻，n为监测截面的个数，x

作为本发明的一种优选方案，步骤S300的具体步骤包括：

S301、设定移动时窗用来从每个模板工况离散点应变时间序列Y

W为t

其中，t为移动时窗所含载荷步数，即时窗长度，其变化范围为[t

S302、设模板工况滑坡中心在坐标轴x'的坐标为d，d的变化范围为[d

x＝x′-d

计算各监测截面在模板工况坐标轴x中的坐标，将这些坐标与离散点的x坐标进行比较，获得与这些坐标最近的离散点；

S303、开始时，计算d

S304、移动时窗每到达一个位置，即得到一个二维矩阵W，下面对移动时窗当前位置的W与X进行动态时间规整，定义距离d

遍历W和X里的各元素，即可得到距离矩阵

对于矩阵d，以d

S305、定义累计距离矩阵

其中Dij为d11为到dij的最小累计距离，计算公式为

到此可知，当前W和X的最小累计距离为D

S306、随着移动窗口在Yi中移动，可得相似度矩阵C

当t在[t

S307、对第i个模板工况的相似度矩阵，比较矩阵中各相似度的大小，求得该模板工况对应的最大相似度，以及相应的滑坡中心位置d，滑坡推力q变化范围；

S308、对M个模板工况S

S309、比较每个模板工况的最大相似度，可以得到所有模板工况下的模板数据库里与监测应变时间序列X最相似的序列片段，则该序列片段对应的工况参数即为辨识出的与实际工况最接近的工况参数。

以上实施例仅为本申请的示例性实施例，不用于限制本申请，本申请的保护范围由权利要求书限定。本领域技术人员可以在本申请的实质和保护范围内，对本申请做出各种修改或等同替换，这种修改或等同替换也应视为落在本申请的保护范围内。