混凝土坝强震观测资料的ARX模型分析

郑东健,罗德河,仇建春

(1.河海大学水文水资源与水利工程科学国家重点实验室,江苏 南京 210098;2.河海大学水利水电学院,江苏 南京 210098; 3. 珠江水利科学研究院,广东 广州 510611)

我国是世界上遭受地震灾害比较严重的国家之一。为了减少地震带来的损害,需要科学地认识强震作用及其对水工结构工程的影响,对水工结构进行强震观测就是解决上述问题的重要途径。迄今为止,世界范围内的混凝土坝已有不少因遭遇强烈地震而破坏的实例。1962年3月19日我国广东新丰江水库发生了MS6.1级强烈地震,使其混凝土大头坝头部发生水平裂缝[1];1967年12月11日印度的Koyna重力坝遭遇了MS6.5级水库诱发地震,使其下游折坡点附近的上、下游坝面出现大量的水平裂缝[2];1971年2月9日和1994年1月17日,美国Pacoima拱坝分别遭遇了MS6.6级和MS6.8级强烈地震,造成左岸坝体与重力墩之间横缝张开[3-4];1990年6月21日伊朗的Sefid Rud大头坝遭遇了MS7.6级大震,坝体受到了严重损伤,上、下游坝面出现大量的水平裂缝[5];1999年9月21日台湾石岗闸坝遭遇MS7.3级强烈地震,北段三跨泄洪道断塌,断裂处错动达8 m[6];2008年5月12日我国沙牌碾压混凝土坝遭遇MS8.0级汶川大地震,拱坝仅坝顶电梯井房发生破坏,其他部位表现良好,没有明显的损伤开裂等现象发生。以上工程震害相对明显,可以通过目视检查发现。但结构内部、水下等隐蔽部位损伤情况常难以目视检查,而一般损伤情况下渗流和变形观测资料也变化不明显,如何及时把握地震过程结构状态变化是目前面临的巨大挑战。随着我国西南地震高烈度区一批高坝大库的完建运行,及时分析强震观测资料,在线快速评估地震过程大坝结构状态变化,对确保大坝安全有重要意义。目前DL/T 5416—2009《水工建筑物强震动安全监测技术规范》第9条加速度记录的处理分析,规定根据强震观测系统“各个测点的记录和预存的抗震设计计算求得的各个测点能抗御的最大加速度值,对大坝进行安全评估”(以下简称“规范法”),并规定“应对加速度记录进行常规处理分析”,其内容包括校正加速度记录、速度和位移时程、5个阻尼值(0,0.02,0.05,0.1,0.2)反应谱和傅里叶谱等。SL 486—2011《水工建筑物强震动安全监测技术规范》第7条加速度记录的处理分析,规定与DL/T 5416—2009类似,这些条款对水工建筑物的震后评价起重要作用。但由于工程的复杂性以及运行环境的影响,设计模型和参数与工程实际情况会存在差异,根据设计计算的测点能抗御最大加速度评估大坝震后安全有其不足。为此,本文直接应用混凝土坝强震观测资料,在强震加速度记录常规处理分析基础上,研究了强震资料时间序列ARX(autoregressive exogenous)模型,通过跟踪ARX模型的时变参数变化,提出混凝土坝结构状态震后评价方法,并采用遗忘因子提高模型效率,实现混凝土坝强震观测资料的在线快速评估。

1 ARX模型

ARX模型可以看作系统的动力学方程,其参数反映了系统的动力特性[7]。混凝土坝强震观测资料是混凝土坝在受地震作用下的振动响应。当地震作用的强度较小时,一般不会对结构的状况产生影响,此时结构系统获得的动力响应较为平稳,强震观测时间序列为平稳时间序列;当地震作用强度较大时,混凝土坝可能会进入非弹性阶段并出现损伤,观测资料为非平稳时间序列。用强震观测资料时间序列建立ARX模型,可以实现用较少参数反映序列所蕴含的结构动力特性,识别混凝土坝结构状态的变化。

本文所采用时间序列ARX模型如下:

式中:y(t)为混凝土坝强震测点采集的信号;u(t)为混凝土坝基础或周边场地强震仪采集的信号;ξ(t)为与误差和模型不确定性有关的白噪声信号;t为时间；ai、bj为系数;na、nb为系统模型的阶次,可根据AIC准则[8-9]确定。

为判断混凝土坝结构在不同时间段的结构参数,在式(1)基础上,将结构参数表示为时间的函数,即:

y(t)=φT(t)θ(t)+ξ(t) (t=1,2,…,N)

(2)

其中φT(t)=(-y(t-1),-y(t-2),…,

-y(t-na),u(t),u(t-1),…,u(t-nb))

θ(t)=(a1(t),a2(t),…,ana(t),b0(t),b1(t),

b2(t),…,bnb(t))T

(3)

其中Y=(y(1),y(2),…,y(N))T

φ=(φT(1)φT(2) …φT(N))T

(4)

P(t)=[I-K(t)φT(t)]P(t-1)

式中：I为单位矩阵。

在采集到结构系统输入输出数据的同时,为实现结构时变参数的实时在线识别,跟踪混凝土坝结构系统不确定性规律的参数变化且提高计算效率,可引入遗忘因子f(t),结合递推最小二乘算法实现时变参数的识别,其主要公式如下[10]:

(5)

式中：f(t)为遗忘因子,计算过程中需满足0

P(0)=αI

(6)

(7)

式中:α为充分大的正实数(104～1010);ε为零向量或充分小的正实向量。

2 实例验证

为验证ARX模型时变参数对混凝土坝结构状态评价的有效性,采用印度Koyna地震波对西部某混凝土拱坝进行结构动力有限元分析,并将拱坝典型位置动力响应作为仿真强震资料。该坝为混凝土双曲拱坝,水库正常蓄水位1 880 m,电站装机容量3 600 MW,坝顶高程1 885.00 m,坝基最低建基面高程1 580.00 m,最大坝高为305.0 m,坝顶宽度16.0 m,坝底厚度63.0 m。图1为拱坝的有限元模型,模型上游、下游和地基取1倍的坝高。模型采用ABAQUS软件8节点六面体实体单元C3D8R进行网格划分,整个模型共有61 909个单元,68 448个节点,其中坝体单元4 708个。

图1 某拱坝有限元模型

为简化计算,将坝基简化为相同材料参数。结构动力分析时采用动弹模,即在静弹模的基础上提高30%作为动弹模进行分析[14],采用无质量地基进行模拟。为了模拟地震作用过程中坝体损伤的发生和发展,采用混凝土塑性损伤模型(CDP模型)进行分析,并通过损伤因子直观的得到结构的损伤位置与程度。坝体和坝基动弹性模量分别取32 GPa和20 GPa,密度分别取2.4 t/m3和2.7t/m3,泊松比分别取0.167和0.25;坝体混凝土材料塑性参数膨胀角取30°,偏心率取0.1,双轴极限抗压强度与单轴受压极限强度之比fb0/fc0取1.16,拉伸子午面上与压缩子午面上的第二应力不变量之比K取0.666 7,黏性参数取0.000 5。

混凝土塑性损伤模型的损伤因子采用张劲公式法[15],该方法将规范提供的混凝土本构模型与CDP模型统一起来,并提出损伤因子dk可表示为

(8)

式中:t、c分别为拉伸和压缩情况;β为塑性应变和非弹性应变之比,在受压情况时取值为0.35～0.7,在受拉情况时取值为0.5～0.95;εin为混凝土在拉压情况下的非弹性阶段应变；E0为坝体初始(无损)弹性模量。dk=0表示无损伤,dk=1表示完全损伤。

进行混凝土拱坝损伤有限元模拟时,输入的地震波采用有完整记录的印度Koyna地震波,取其顺河向地震波作为输入波,地震加速度时程曲线如图2所示。数值仿真模拟的采集时长取10 s,采样频率取50 Hz。为了分析不同地震强度下拱坝结构状态变化,将加速度峰值分别调整为0.1g和0.2g的两条Koyna地震波时程曲线,作为结构的地震激励输入,用ABAQUS有限元软件模拟相应地震强度下拱坝损伤的发展,并分析其与ARX模型参数变化的对应关系。

图2 水平向地震加速度时程曲线

为了提取各地震工况拱坝的动力响应时间序列,在拱冠坝顶、1/2坝高、左右岸1/4拱坝处分别设置动力响应采集测点,编号为1号～4号,具体位置如图3所示。建立ARX时间序列模型时,采用下游河床坝基位置的一个测点响应信息作为ARX时间序列模型的响应输入。

图3 拱坝强震测点仿真布置

下面通过模拟两种不同工况下的损伤情况,分析各测点加速度响应时间序列建立的ARX模型参数和结构状态变化的对应关系。

a. 峰值加速度为0.1g。当结构输入加速度峰值为0.1g的地震激励时,有限元模拟仿真的结果表明,坝体总体上保持完好,没有发生损伤。坝基测点、1号、2号、3号测点加速度时程曲线见图4。将坝基位置测得的加速度响应信息作为输入信息,各个测点的加速度响应信息作为输出建立ARX模型(na,nb分别取6和2)。图5给出了1号测点ARX模型的时变参数辨识结果(2号～4号测点类似)。从图5可以看出,测点的响应时间序列计算的系统参数值在整个地震期间保持稳定不变,说明结构系统状态没有发生变化,即拱坝在小震情况下没有发生损伤,与仿真分析结果一致。由于递推迭代原因,迭代初始阶段参数变化可能较大,可通过预先加1～2 s白噪声激励(幅值与地震波初始振幅接近)来消除。

图4 测点响应时间序列

图5 ARX模型参数辨识结果

图6 峰值加速度为0.20g时坝体损伤状态

b. 峰值加速度为0.20g。为了使拱坝出现损伤,将峰值加速度为0.2g的Koyna地震波作为地震动作用输入,结构的破坏情况如图6(a)所示。从图中可以看出拱坝的损伤区域集中在拱冠、两岸坝肩及坝踵处,拱冠位置上下游面均发生破坏。其中拱坝的中上部位区域破坏的较为严重。在地震作用的后期,破坏区域逐渐向两侧和坝体下方扩展。本文模拟的大坝损伤结果与文献[16]结果较为一致。三个损伤区域的拉伸损伤因子时程图如图6(b)所示,可以清晰地看出各区域的损伤情况,此时结构发生破坏的时间不集中,并且损伤因子值呈阶梯状增长。如坝踵部位,4.26 s首先发生损伤,到4.9 s时损伤因子值为0.268,但损伤程度较小,到6.2 s时损伤因子值突变为0.775;左岸坝肩1/4拱位置在7.44 s时损伤因子突变为0.849;拱冠位置4.62 s时损伤初始,6.52 s时损伤因子值为0.191。将坝基位置测得的加速度响应信息作为输入信息,各个测点的加速度响应信息作为输出建立ARX模型。ARX模型时变参数的辨识结果如图7所示，可以看出,1号、3号和4号测点ARX模型参数在4～5 s之间发生变化,尤其1号测点突变明显;在6～7 s间各测点ARX模型参数再一次出现较明显变化,2～4号测点突变最明显;随后各测点模型参数仍有不同程度变化,但幅度有所减小。对比图6(b)和图7可见模型参数突变时段损伤因子值也发生了突变,且不同位置测点识别出的结构损伤发生时间与仿真模型损伤因子值突变时间基本一致,表明混凝土坝强震资料的ARX模型参数具有较好的结构状态识别能力,且坝顶等动力响应较大的位置测点对损伤初始判别更为敏感。

图7 ARX模型参数辨识结果

3 结 语

我国许多高混凝土坝建在高地震烈度区,在经受各种程度地震动荷载作用下,这些高坝是否出现损伤直接关系到大坝安全和震后应急措施。本文直接利用实测强震观测资料,通过建立ARX模型，利用模型参数反映结构动力特性的原理，结合遗忘因子递推最小二乘算法,可以实现混凝土坝震后结构状态的快速评估。混凝土拱坝地震作用仿真分析表明,ARX模型参数的时变过程与损伤因子的变化过程较为一致,验证了时变参数识别方法的有效性。强震观测资料的ARX模型分析法具有快速和便于在线评价的特点,可以作为规范法的补充。