基金项目:中国地震局地震科技星火计划项目“青藏块体东北缘前兆系统结构分析与异常识别”(XH12072)与陕西省科学技术发展计划项目“汾渭断裂带地壳运动的非均匀负位错模拟与应变积累研究”(2011JM5016)联合资助.
(Second Crust Deformation Monitoring and Application Center,CEA,Xi'an 710054,Shannxi,China)
fixed-point precursor abnormity; spreading function; inversion; earthquake risk area
备注
基金项目:中国地震局地震科技星火计划项目“青藏块体东北缘前兆系统结构分析与异常识别”(XH12072)与陕西省科学技术发展计划项目“汾渭断裂带地壳运动的非均匀负位错模拟与应变积累研究”(2011JM5016)联合资助.
1998年张北6.2级地震前,山西、河北及其附近地区出现的定点形变、流体前兆异常,具有随时间推移从震区外围大范围区域向震区附近迁移、集中的趋势特性。通过构建异常信号传播函数,借助网格搜索—拟牛顿最小二乘法反演汇聚区域即推测可能的地震危险区,所得汇聚区域边缘距实际震中仅25~51 km、汇聚时刻与发震时刻仅差数天,且震前中短期阶段汇聚速度加快。
Precursory anomalies of deformation and fluid measured in fixed-point appear the trend features of moving and concentrating from distance areas to the near of seismic area before the Zhangbei MS6.2 earthquake in 1998.Constructing abnormal signal spreading function,and using grid searching-quasi Newton least square method,we inverse the converged region of precursory anomalies which is the predicting possible risk region.The result shows that the distance between the margin of calculated converged region and the real epicenter is only 25-51 km.The difference between the converging time and the original time is only several days.Moreover,the converging rate speeds up at the middle-short term before Zhangbei MS6.2 earthquake.
引言
我国大陆中强以上地震尤其强烈地震前中—短期阶段,震区及其周围一定范围内或多或少会出现由在定点前兆台站观测到一些异常变化(孟国杰,黎凯武,1999; 石绍先等,2004; 张立等,2005; 徐辉等,2007; 刘翔等,2008; 张晓清等,2008; 李惠玲等,2009; 林洁萍,李佳,2009; 牛安福等,2009)。观测前兆异常的手段、测项不同(形变、流体、地电等),时间上可能同步、准同步,也可能分中期、短期或短临等不同阶段(孟国杰,黎凯武,1999; 石绍先等,2004; 尹亮等,2005; 刘翔等,2008; 张晓清等,2008; 李惠玲等,2009; 牛安福等,2009)。有些地震前还出现明显的、同类或多类手段测项异常,如由震区外围稍大范围区域向震区或其附近局部构造区域迁移或包围、集中的趋势特征(Scholz,1977; 汪成民等,1990; 孟国杰,黎凯武,1999; 龙海英等,2007; 徐辉等,2007; 张晓清等,2008; 李惠玲等,2009),与梅世蓉和朱岳清(1993)对唐山地震等的研究中的硬包体模式比较相符——即震源是地层中镶嵌着的一个硬包体,异常首先在震源的外围出现,在孕震过程的后期硬包体产生弱化,震源区或其附近则出现异常。
如果能构建数理模型模拟上述异常的时、空迁移轨迹、规律,则可能对中强以上尤其强震地点,甚至发震时间预测有所帮助。本文利用1998年1月10日张北6.2级地震前数月至一、二年山西北部、河北及其附近地区定点形变、流体前兆异常资料,探索构建异常信号传播函数,模拟其随时间推移从震区外围数百公里区域向震区附近迁移、集中的时空演化过程,借助网格搜索—拟牛顿最小二乘法反演、寻找前兆异常汇聚区域,即推测可能的潜在地震危险区。考虑硬包体模式(梅世蓉,朱岳清,1993),本文将异常汇聚(收缩)区设定为一个具有圆心、一定半径的圆形空间区域,较之单点假设更合理,比椭圆等其它形状简化,便于模型构建和反演计算。
1 信号传播模型构建与反演方法
为模拟中—短期(震前数月至一、二年)定点台站前兆异常向某区域(假设为潜在的地震危险区)集中、汇聚(收缩)或迁移的时空演变规律,设定式(1)为含6个参数的目标函数。6个参数依次为汇聚区域圆心经度、纬度的高斯投影平面坐标x0、y0; 异常信号传播速度a(负值); 到达汇聚区边缘的时刻t0; 借助Taylor展开式获得的类似加速度的参量b; 汇聚区域的半径R。其中异常台站数m≥6,其经纬度的高斯投影坐标为xk、yk,k=1,…,m; tk为异常开始时刻,k=1,…,m。
f(x0,y0,a,t0,b,R)=1/2∑mk=1[((x0-xk)2+(y0-yk)2)1/2-R-a(tk-t0)-b/2(tk-t0)2]2.(1)
即利用异常台站至汇聚区边缘由时间、速度和加速度确定的函数(式(1)右括号中后两项),近似异常台站与汇聚区边缘距离( 式(1)右括号中前两项),最终寻找合理的参数使目标函数(1)最小。而拟牛顿最小二乘法(Matsu'ura等(1986)以及张希等(2003)的文献中不考虑参数初始中误差的情况)在参数初值相对接近最优解时,比遗传算法收敛效率更高。
假设x为待定模型参数向量; y、f(x)、e分别为观测值向量、理论计算函数向量与残差向量; 定义目标函数为T=1/2eTP-1e,残差e=y-f(x)y-AkXk,为观测值与第k次迭代解xk(线性化模型,通过对目标函数求偏导获得)理论值向量之差,P为观测值方差矩阵,本文P取单位矩阵,即认为各站点观测值等权,则迭代解按式(2)计算
xk+1=xk+αk(ATkP-1Ak)-1rk.(2)
其中,0<αk≤1; rk=ATkP-1[y-f(xk)].
拟牛顿最小二乘法终止准则为两次迭代目标函数之差(本文规定绝对值10-10以下)与偏导数(本文计算结果的绝对值均在10-8以下)趋于零。具体计算时,参数a,t0,b,R的初值根据异常迁移特征、规律如估算平均速率等信息确定; 在方圆数百公里较大范围内寻找最优的前兆异常汇聚区的圆心位置,其高斯坐标初值(x0,y0)由给定经、纬度范围内等间隔网格点投影转换而确定,相当于对数十对网格点坐标初值分别进行拟牛顿最小二乘反演,搜索所有反演所得参数中使目标函数(1)最小的值作为最终结果。
2 张北震前实测资料反演分析
1998年1月10日,晋、冀、蒙交界附近河北省境内张北地区发生6.2级地震(41.10°N,114.30°E)。孟国杰和黎凯武(1999)给出了该次地震前东侧(河北地区)约300 km范围内6个定点形变异常,分别是兴隆台水准(1996年8月起出现异常,震中距276 km; 宽城台水准于1996年1月起出现异常,震中距290 km; 易县台水准于1997年1月起出现异常,震中距229 km; 张家口台地倾斜1997年10月起出现异常,震中距60 km; 宝昌台地倾斜于1996年7月起出现异常,震中距120 km; 赤城台地倾斜于1997年11月起出现异常,震中距129 km。从其时、空分布来看,前兆异常先在震区外围出现,总体上具有由远至近迁移、集中的趋势特性。
李惠玲等(2009)则给出了相近时段、张北震前山西北部3个流体前兆异常,即大同镇川井水位在1997年6月起出现异常,震中距175 km; 而大同机车厂在1997年12月出现水氡异常,反演时是出现得最早的一个水位异常; 静乐娘子神井水位在1995年7月起出现异常,震中距350 km; 定襄七岩泉1996年4月起出现水氡异常,震中距318 km。时间上有由南向北迁移,即有与河北地区相似的异常从外围向震区附近靠近的趋势。异常台站分布见图1,具体异常观测曲线参见孟国杰和黎凯武(1999)、李惠玲等(2009)的研究,本文略。
将这两个地区共9个前兆异常(震前两年半内)进行综合反演。首先,从异常出现时间上看,张家口台站和赤城台站的异常出现得最晚(于1997年10或11月才出现),可能是因为距汇聚区更近,本文以这两个台站的中点(40.86°N,115.34°E)为中心,将汇聚区圆心网格搜索范围向东、南、西、北各延伸2°,即在其周围4°×4°(约300 km)范围内取8×8个网格点(经高斯投影后作为汇聚区圆点坐标初值)进行搜索; 汇聚时刻初值为零,规定1998年1月1日为零,以月为单位即1997年12月为-1.0,这样兴隆水准异常出现时间记为-17.0,其它台站异常出现时间依此类推。其次,考虑到山西地区静乐台至大同台距离236.9 km,异常出现时间间隔23个月,汇聚速度值约为-10.3 km/月; 河北地区宽城台至张家口台距离303.7 km,异常出现时间间隔21个月,速度值约为-14.5 km/月,两地区平均汇聚速度-12.4 km/月,取反演所用汇聚速度初值为-12.4 km/月。此外,加速度初值定为零; 鉴于6级地震的破裂尺度约为10km(Wells,Coppersmith,1994),一般认为孕震区约为破裂尺度的数倍(Keilis-Borok,Rotwain,1990; 彭克银等,2003),因此将区域半径初值定为50 km。模拟所得汇聚区1圆心(40.73°N,114.99°E)位于震中东南约72 km处,区域半径为46.6 km,汇聚速度为-9.3 km /月; 到达汇聚区边缘的时刻为1998年1月9日前后(图1、表1),与张北地震发震时刻非常接近,汇聚区边缘距实际震中约25 km。
如果只考虑对震前一年半内,保证至少6个异常台站,基本属于震前中短期阶段,即对1996年下半年以来出现异常的6个台站进行反演,速度初值调整为-16.0 km /月(山西地区只剩大同台; 河北地区兴隆台至张家口台距离223.5 km,异常出现时间间隔14个月,速度值约为-16.0 km/月),所得模拟汇聚区2的圆心(40.50°N,115.16°E),位于震中东南约98公里处,汇聚区域半径47.3 km; 速度为-14.7 km /月; 到达汇聚区边缘时刻为1998年1月18日前后(图1、表2),与张北地震发震时间也比较接近,汇聚区边缘距实际震中约51 km。而汇聚区1圆心与汇聚区2圆心相距仅29 km,大多数参数量值、特性接近,故可以近似地认为其反映的是一个潜在危险区。
图1 张北震前河北—山西及其附近异常台站分布及模拟汇聚区1、2的位置
Fig.1 Distribution of abnormity stations and location of No.1 and No.2 inverted converged regions in and near Hebei-Shanxi before the Zhangbei earthquake如果只考虑河北及其附近6个定点形变异常(流体异常较少),速度初值调整为-14.5 km/月(宽城台至张家口台距离303.7 km,异常出现时间间隔21个月,速度值约为-14.5 km/月),其它参数初值不变,所得模拟汇聚区3圆心(40.99°N,115.62°E)位于震中以东约112 km处,汇聚区域半径46.4 km; 汇聚速度为-12.2 km /月; 到达汇聚区边缘时刻为1997年12月31日前后(图2、表3),与实际发震时间也比较接近。汇聚区边缘距实际震中约66 km,距离汇聚区1、2的圆心分别约61 km、68 km,也有一定预测意义。但就本文震例来看,综合两种资料异常所得汇聚区距实际震中更近、地点预测效果相对更好。
注:6个台站中里公式(1)中括号内差值的绝对值小于或等于拟合中误差的占3/6.3 结论与讨论
(1)1998年1月10日张北6.2级地震前两年半时段内山西北部、河北地区出现的定点形变、流体前兆异常,具有随时间推移从震区外围大范围区域向震区附近迁移、集中的趋势特性。通过构建具有一定速度、加速度,向某一区域收缩的传播函数,借助网格搜索—拟牛顿最小二乘法反演所得汇聚区域边缘距实际震中仅25~51 km、汇聚时间与发震时间比较接近(最短数天); 且震前中短期阶段(一年半)相对震前中期阶段(两年多)汇聚速度加快,具有一定程度的中短期甚至短临地点与时间预测意义,以及对多手段、多测项前兆异常的联合计算意义。
(2)本文仅是针对1998年张北地震事件模拟其显著异常时空迁移规律的数值方法探索中的一种,是否通用于其它地震前兆,还需今后更深入、广泛地研究。
- 李惠玲,高云峰,张登科,等.2009.1998年1月10日河北张北6.2级地震前山西北部流体观测前兆异常分析[J].山西地震,(4):12-14.
- 林洁萍,李佳.2009.2007年6月3日云南宁洱6.4级地震异常特征分析[J].内陆地震,23(1):36-41.
- 刘翔,赵盼盼,和宏伟,等.2008.滇西地区前兆观测资料的分类、分级与强震关系研究[J].地震研究,31(2):103-108.
- 龙海英,聂晓红,孙甲宁,等.2007.新疆3次6级强震前的短期前兆异常[J].内陆地震,21(4):304-310.
- 梅世蓉,朱岳清.1993.唐山地震孕育模式研究的科学思路和结果[C]//梅世蓉.梅世蓉地震科学研究论文选集.北京:地震出版社,381-390.
- 孟国杰,黎凯武.1999.张北6.2级地震前的地壳形变特征[J].地震,19(3):261-266.
- 牛安福,张凌空,闫伟,等.2009.汶川地震前南北地震带中北段地形变变化特征的研究[J].地震,29(1):100-107.
- 彭克银,尹祥础,和锐.2003.用临界点理论讨论应变能加速释放现象和孕震区尺度[J].中国地震,19(4):425-430.
- 石绍先,曹刻,徐彦,等.2004.云南楚雄地区成组强震与单个强震的前兆差异性研究[J].地震研究,27(2):111-118.
- 汪成民,孙振璈,车用太.1990.大同-阳高地震的地下水异常初析[J].地震,10(4):37-44.
- 徐辉,白亚萍,康庆强,等.2007.岷县5.2级地震前兆异常分析[J].高原地震,19(3):25-28.
- 尹亮,李东生,周丽萍,等.2005.民乐、山丹6.1级地震异常及台站实现短临预报的思考[J].西北地震学报,27(21):112-116.
- 张立,牛安福,汤曙恩,等.2005.大姚地震前连续形变短期前兆异常[J].地震,25(2):75-82.
- 张希,江在森,王琪,等.2003.1999~2001年青藏块体东北缘地壳水平运动的非震反位错模型及变形分析[J].地震学报,25(4):374-381.
- 张晓清,张晓香,秦松涛,等.2008.短期前兆综合异常方法在青藏高原北部地区的应用——以兴海地震为例[J].西北地震学报,30(2):179-183.
- Keilis-Borok V I,Rotwain I M.1990.Diagnosis of time of increased probability of strong earthquakes in deferent regions of the world[J].Phys Earth Planet Inter,61(1-2):57-72.
- Matsu'ura M,Jackson D D,Cheng A.1986.Dislocation model for aseismic crustal deformation at Hollister,California[J].JGR,91(B12):12 661-12 674.
- Scholz C H.1977.A physical interpretation of the Haicheng earthquake prediction[J].Nature,265(5607):121-124.
- Wells D L,Coppersmith K J.1994.New empirical relationships among magnitude,rupture length,rupture width,rupture area,and surface displacement[J].BSSA,84(4):974-1 002.