用椭圆衰减关系模型计算任意场点烈度及地震动参数的数值方法

1、文章编号:1001-8956(201101-0021-08中图分类号:P3153文献标识码:A用椭圆衰减关系模型计算任意场点烈度及地震动参数的数值方法石建梁1,闫庆民1,葛秋莹2(陕西省地震局,陕西西安710068;西安中勘工程有限公司,陕西西安710016摘要:用椭圆衰减关系计算地震动或地震烈度,比较经典的做法是建立极坐标系,采用极坐标方程的方法。给出了一种在直角坐标系中即可求解的用椭圆衰减关系计算地震动或地震烈度的数值计算方法。另外,地震动椭圆衰减关系,应满足一定的模型要求,才能更客观地反映实际的地震动衰减物理过程。本算法的附加功能能够鉴别非标准的椭圆衰减关系。关键词:标准椭圆衰减关系;超

2、越方程;数值解;二分法;椭圆震级极限轴长;直角坐标系依据我国的地震烈度等震线统计资料或强震记录,用椭圆衰减关系模型描述地震所造成的烈度等震线分布或地震震动参数(加速度、速度或反应谱等等震线分布,相对于用点圆衰减关系模型表述更为客观合理。但对于某一给定的地震,计算其对任意场点造成的影响烈度或地震动值时,椭圆模型方法远比点圆模型方法复杂,需要确定场点所在位置的影响烈度或地震动值的极坐标方程。为此本文中给出了一种数值算法,比较好地解决了这个问题,即在直角坐标系中求解,并用历史地震和现代仪器记录分别对本算法进行了验证。本算法中嵌入的附加模块,还能够鉴别非标准的椭圆衰减关系。1求解公式椭圆衰减关系的通用

3、形式为长轴log Y a =C 1+C 2M +C 3M 2+C 4log R a +C 5exp (C 6M ;(1短轴log Y b =D 1+D 2M +D 3M 2+D 4log R b +D 5exp (D 6M (2式中,Y a 、Y b 是烈度或地震动值;M 是震级;C 1 C 6、D 1 D 6是回归统计得到的系数,R a 、R b 是震中距。对于烈度衰减关系,log Y a 项为I a ,log Y b 项为I b ,C 6、D 6为0,C 5、D 5分别为确定的统计值。log 一般为lg 或ln 。对于某一给定的地震,其震中经纬度、震级M 和方向角(椭圆长轴与正东方向的夹角

4、已知。对任一确定的场点,其经、纬度已知。取震中为原点,建立以正东、正北方向为X 、Y 轴的笛卡尔右手直角坐标系。通过坐标系第25卷第1期2011年3月内陆地震INLAND EARTHQUAKEVol25No1Mar2011*收稿日期:2010-02-09;修回日期:2010-04-01作者简介:石建梁(1965 ,男,陕西西安人,高级工程师,工学硕士,毕业于中国地震局工程力学研究所地震工程与防护工程专业,主要从事工程地震与岩土工程勘察及人工地基工程质量检测工作E-mail :shjl65yahoocomcn的旋转变换,将X、Y轴逆时针旋转角度,可得到场点在以椭圆长短轴为X、Y轴的新的直角坐标系

5、中坐标(x c、y c。于是可得以下两个方程x2 c R2a +y2cR2b=1;(3C 1+C2M+C3M2+C4logRa+C5exp(C6M=D1+D2M+D3M2+D4logRb+D5exp(D6M(4在以上两个方程中,有两个未知数R a和R b,求出任一个,利用(1或(2式,即可求得该场点的烈度或地震动值。设一个代数Q,当(1或(2式中的log为lg时,Q取10;当(1或(2式中的log为ln时,Q 取e(271815。则由(4式得Rb =(Ra+C5e C6MC44Q C1D1+(C2D2M+(C3D3M24D5e D6M(5令A=Q C1D1+(C2D2M+(C3D3M2D4;B

6、=C5e C6M;E=D5e D6M;F=C4D4。则(5式变为Rb =(Ra+BF AE;(6对(6式取平方得R2b =(Ra+B2F A22AE(Ra+BF+E2;(7代(3式入(7式中得y2 c R2aR2a x2c=(Ra+B2F A22AE(Ra+BF+E2;(8这样就得到以R a为未知数的一元实系数超越方程,即f(Ra =y2cR2aR2ax2c(Ra+B2F A2+2AE(Ra+BFE2=0(92超越方程的数值解超越方程有没有根和有多少根的数学问题到目前还没有完全解决,而且也没有解析解。常用的数值解法有二分法、截弦法、优选法、插值法、梯度法和拟牛顿法等。每种方法各有其特点和要求,

7、到现在还做不到使某一种算法对任何方程都是有效的。所以应根据所要求根的方程的具体特点,例如根据方程的特性(是否含有复系数等、未知数的多少、是否有复数根、初值点的好坏等灵活使用这些算法。这里使用的是二分法,它是针对实系数超越方程的一种区间套法。在求根时,要求给出求根区间的端点值,用不断紧缩区间的办法求出方程在某一个区间内的全部单重实根。式(9显然是实系数的超越方程,且未知数R a或R b是待求的椭圆的长轴或短轴,必然是方程(9的实根,故选用二分法较为合适。22内陆地震25卷二分法的具体步骤如下1。(1设定根存在的区间,即给出根之下界AA 与上界BB ,给出计算步长H ,给出函数值f (R a 的精

8、度E 1及根的精度E 2。(2自AA 开始a =AA ,以H 为基本步长,向右跨出宽为H 的小区间a ,b 。当跨得的小区间两端函数值同号时继续向右跨;异号时,则用二分法找出里面的一个实根,即反复计算f (a +b2,直至满足求根原则。(3求根原则是函数值满足给定的函数值的精度E 1或相邻两次根的计算值的差值的绝对值满足根的精度E 2。(4重复以上步骤,直至小区间的上界大于BB 。下面针对我们所面对的具体超越方程(9,关于二分法的输入参数AA 、BB 、H 、E 1及E 2做如下分析。(1下界AA 。对于给定的震源与场点,其间的距离RR (即震中距可以求得。对于椭圆长轴R a ,必有R a R

9、R ,故取AA =RR 。(2上界BB 。我们知道,标准的地震烈度或地震动参数椭圆衰减关系的一个重要特征是,在衰减关系的起始点是重合的,在远场区近似圆,在近场区为椭圆。即对于确定的震级M ,椭圆模型的长轴有个上限长度DD (图1中交点所对应的震中距。当震中距为0或DD 时,地震动或烈度值相等;当震中距为0至DD 的中间时,地震动或烈度再无相等的值。震中距 大于图1标准椭圆衰减关系模型适用的震中距范围图Fig1The suitable epicenter distance extent for standard ellipsoid attenuation model321期石建梁等:用椭圆衰减关

10、系模型计算任意场点烈度及地震动参数的数值方法DD ,椭圆模型失效,如图1所示。这里我们把DD 定义为椭圆震级极限轴长。DD 可以通过求解一个新的超越方程(10得到。令(4式中R a =R b ,得到g (R a =A (R a +B FR a E =0(10该方程的求解方法与方程(9的相同,不同之处仅仅是上限的选取。这里结合地震工程的实际特点,上限取为1000km 。(3步长H 。衰减关系中的R a 与R b 以及计算流程中的长度的单位均是km ,故取H =001,即10m 。作者收集了10数条地震烈度或加速度峰值衰减关系,设定不同的震源点与场点位置,设定不同的震级M 与方位角,从而构成了数百

11、种关于超越函数f (R a =y 2c R 2aR 2a x 2c(R a +B 2F A 2+2EA (R a +B FE 2(11的组合。在参数上下限AA 与BB 范围内,其函数图像无一例外如图2所示。图3为在方程根附近的函数图像。从图2、3中可以看出,首先实际地震影响的烈度或地震动参数值的椭圆长轴R a ,短则数公里,长则数百公里,所以方程求根的区间宽度为数10至数百公里,以001km 的步长计算,精度是足够的。其次,函数图像没有呈现类似地震记录那样的随机、高频震荡的形式,而是平缓的递减函数,故取H =001,分辨率是充分足够的。关于f (R a 是递减函数,证明如下 。图2超越方程的函

12、数图像(根的下限以上范围Fig2The functional image of surpass equation (above lower limit of roots 函数f (R a 的一阶导数为f '(R a =y 2c x 2c(R 2c x 2c22FA (R a +B F 1A (R a +B FE (1242内陆地震25卷由椭圆模型的性质,其参数必满足F 1,B E 0,且A 一般不小于02,则当R a RR 时,f '(R a 0,即函数f (R a 为单调减函数。也就是说取H =001,分辨率是充分足够的,且方程往往只有1个根 。图3超越方程的函数图像(根的附

13、近范围Fig3The functional image of surpass equation (near roots (4精度E 1及E 2。取超越方程算法中常用的E 1=105及E 2=105或E 1=104及E 2=104,完全满足我们这种地震工程问题的精度要求。另外,用二分法求解超越方程根时,输入参数中有个“欲求根的个数”问题。尽管通过以上论证及作者所进行的近百个算例说明,在前述R a RR 范围内的实际根的数量都只有1个,但为了方法与概念上的完备与充分,也为了能够应对超越方程(9出现一个以上根的情况,在算法即所编制的FORTRAN 程序中,“欲求根的个数”一项都取为3。对于可能出现的

14、1个以上根的情况,采用了比较椭圆长短轴长度的筛选办法。若某个根所对应的长轴的长度小于将其代入(3式中求得的短轴的长度,则此根为增根,舍去。若仍有一个以上的根,则将不同的等震线长、短轴长度及由此计算出的地震动参数结果全部打印出来,结合实际情况,进行进一步分析和舍取。需要补充说明的是,超越方程(9的实根不止1个。由于函数f (R a 在R a =x c 处存在奇异点,往往会出现一个与场点坐标x c 绝对值大小一样的增根,它明显不满足方程(4。这也说明将方程(9的求解下界取为RR 的合理与必要。计算方法的FORTRAN 程序流程如图4所示。3算例(12009年11月5日陕西临潼高陵交界处发生M 44

15、地震,震中烈度 +。其度521期石建梁等:用椭圆衰减关系模型计算任意场点烈度及地震动参数的数值方法26 内 陆 地 震 25 卷 图4 Fig 4 计算方法的 FORTRAN 程序流程图 The FORTRAN process program for this numerical solution method 区长轴 NWW 向， 长约 24 km; 短轴近 NNE 向， 长约 16 5 km。 2 长安台记录的基岩地震动最大加速度见表 1 。对长 安 台所在 场 点， 用 文献 中 给 出的 中国西部衰减关系式 lnY a = 5 912 025 + 1 836 588 M 2 846 5

16、8ln R a + 3 400exp( 0 451 M ; lnY b = 2 509 012 + 1 360 759 M 1 791 51ln R b + 1 046exp( 0 451 M 计算得到的基岩水平地震动最大加速度见表 1 所示。 表1 数值类型 记录值 计算值 ( 13 ( 14 基岩场地实际地震记录与计算结果的比较 东西向 2 0 南北向 1 3 垂直向 1 6 水平向 Gal 2 4 ( 矢量合 2 7 ( 最大值 从表 1 数 据 看 出 ， 计算结果与实际地震记录值基本吻 合。 考 虑 到 M4 4 地震 不 是 强 震， 而 衰减关系是由强震记录统计、 映射得到的，

17、表中的偏差是可以理解的。 ( 2 1556 年陕西华县 8 烈度衰减关系( 汪素云 4 1 3 级地震， 等震线如图 4 所示 ， 长轴方向角约 35 。 用 中国 西 部 4 计算， 得到相关场地的地 面 烈 度 与 宏 观 烈 度 ( 表 2 。考 虑 到 烈 度值 ; ( 15 ( 16 跳跃性的不连续特点， 计算结果与宏观烈度资料十分吻合， 相应的基岩加速度峰值也合理。 I a = 5 253 + 1 398 M 4 164lg( R a + 26 I b = 2 019 + 1 398 M 2 943lg( R b + 8 2009 陕西省地震局 2009 年 11 月 5 日临潼高

18、陵交界 4 4 级地震现场工作报告， 2009 陕西省地震局 陕西省地震局强震台网资料， 1期 石建梁等: 用椭圆衰减关系模型计算任意场点烈度及地震动参数的数值方法 27 图5 Fig 5 1556 年华县 8 1 级地震等震线图 4 The isoseismic lines of Huaxian earthquake in 1556 4 扩展功能 考虑到工程界许多正在使用中的所谓椭圆 衰减关系， 由于 各 种原 因 并 不 都 是 如 衰减关系 ( 12 那样严格符合椭圆 模 型 的 概 念， 式( 11 、 为了 鉴 别 并 且 使 这 些 衰减关系 能 够继 续 用 以 估 算所要求的目

19、标值， 本方法做了以下功能扩展。 表2 场点地名 华县 渭南 临潼 潼关 朝邑 蓝田 西安 蒲城 富平 1556 年华县 8 1 级地震宏观烈度与计算结果的比较 4 场点地名 临晋 夏县 咸阳 户县 汾城 绛县 灵宝 铜川 鳞游 宏观烈度 /度 计算烈度 /度 8 0 7 4 8 0 7 9 7 0 7 0 7 8 7 7 7 0 宏观烈度 /度 计算烈度 /度 10 5 9 7 9 0 8 7 8 9 8 8 8 4 8 4 8 3 ( 1 输出椭圆震级极限轴长 DD 值， 若 DD 不存在， 或 DD 50 km ( 衰减关系是针对 强 震 而言的。实际资料表明， 即使 4 级地震其 椭圆

20、 震级 极 限 轴 长 也 不 可能 小 于 50 km ， 或在不同 的震级下 DD 为同一个值， 都说 明 该 椭圆 衰减关系 不 完 备， 不 是 标 准 的 椭圆 衰减关系。 对此， 本方法在所编制的 FORTRAN 程序中分别加以识别和处理， 仍能够按 照 原 衰减关系 继 续 计 算， 28 内 陆 地 震 25 卷 并输出 DD 值， 表明不是标准的椭圆衰减关系。 ( 2 输出震中距 RR = 0 时， 分 别 按长短 轴 衰减关系计 算 得 到 的 烈 度或 地震动 参 数 值， 若 二者不相等， 则该椭圆衰减关系不完备， 不符合标准椭圆衰减关系模型。 ( 3 当震中距 RR

21、DD， 即椭圆衰减关系衰变 为 圆 衰减关系 时， 分 别将 RR 输 入 长 轴 与 短 轴衰减关系进行计算， 取其算术平均值作为欲求的烈度或地震动值。 5 结语 本算法符合椭圆模型的数学、 物理概念， 满足烈 度或 地震动 椭圆 衰减关系 的地震 学 规 律， 过程简单直接。算例表明， 与实际地震影响记录吻合。并能够鉴别非标准的椭圆衰减关系， 有 扩展功能。 参考文献: 1 郭富印 FORTRAN 算法汇编 M 北京: 国防工业出版社， 1981 2 卢寿德 GB177412005 “工程场地地震安全性评价” M 北京: 中国标准出版社， 2006 宣贯教材 3 中国地震局震害防御司 中国

22、历史强震目录( 公元前 23 世纪至公元 1911 年 M 北京: 地震出版社， 1995 4 俞言祥 中国东部和西部地区水平向基岩加速度反应谱衰减关系 J 震灾防御技术， 2006 ， 1 ( 3 : 4652 AN ALGORITHM FOR ARBITRARY ENGINEERING SITE EARTHQUAKE INTENSITY OR MOTION PARAMETER USING ELLIPSOID ATTENUATION MODEL SHI Jianliang1 ，YAN Qingmin1 ，GE qiuying2 ( 1 Shaanxi Earthquake Administration，Xi an 710068 ，China; 2 Xi an Central Reconnaissance Engineering Co ，Ltd，Xi an 710016 ，China Abstract: The