论文部分内容阅读
摘要：利用减组合滤波方法对2008—2013年云南省内7次M≥55地震前3个月10个台站的水管倾斜观测资料进行处理，消除潮汐改正成分和零漂，探寻短临同步异常变化。结果表明，减组合滤波方法提取的高频成分，在震前存在同步异常，可尝试用于短临跟踪工作中；利用“Nakai拟合模型”求解得到潮汐响应率和线性漂移速率在短临阶段均未表现出易识别的异常，其反映的信息尚需进一步挖掘。
关键词：减组合滤波；Nakai拟合；水管倾斜；短临异常
中图分类号：P31563文献标识码：A文章编号：1000-0666（2017）01-0082-06
0引言
云南是我国破坏性地震频发的区域之一，仅2014年就发生3次M≥6地震。随着数字化形变观测的发展，研究人员可以从整点值、分钟值提取更多的短临前兆异常信息（刘仲全，1997，2001；付虹，刘丽芳，2003；Harrison，1976；刘序俨，张雁滨，1991；李杰等，2003；张晶等，2004），本文尝试寻找水管倾斜形变资料短临阶段潮汐及非潮汐变化特征，并浅析其与地震的关系。
地壳形变是地震发生过程中最直接的伴随现象（陈德福，1993），它与地震的关系最为直接，因此历来受到国内外地震预报探索者们的高度重视，也常常作为地震预测的重要判据之一，尤其在云南，地壳形变作为中短期预测依据已有很多应用实例，比仅用于趋势异常判定的模拟观测资料更具现实意义。识别形变观测资料异常的难点在于排除干扰，有的干扰能定量排除，有的只能定性排除，对于这些干扰因素规律性的探索及其观测值之间函数关系的研究，引起了许多学者的关注。
由于异常识别问题过于复杂，各个观测点表现出明显的地区性差异，因此这些研究未能获得显著的进展（吴翼麟，周景明，1981）。为了降低干扰排除的难度，对几个台站同步异常的提取同时进行观测分析是比较简便并且可信的方法，可排除局部环境所造成的影响，增加异常的可信度，同时假设异常信息来自震源体，在多个台站同步出现异常信号也容易从机理上得到解释。
目前，云南省共建有形变观测台站22个，其中，有水管倾斜仪的台站有11个，最南端的勐腊台因为距离研究地震的样本较远，将其剔除，不参与分析。多年的观测实践表明，水管倾斜仪记录的资料连续性、稳定性较好，能记录到清晰固体潮及年变趋势。而固体潮是唯一能准确计算出理论值的干扰（郗钦文，1982；郗钦文，侯天航，1986），对于已知的固体潮，可以采用均值法、别尔采夫滤波和组合滤波等方法剔除潮汐成分，提取高频信息。也可以利用固体潮理论值对台站观测值进行拟合，得到反映地壳介质的特征参量，通过对其分析研究，有可能了解到地壳局部岩石弹性性质的变化。
本文主要利用观测质量相对较好的水管倾斜仪资料，采用减组合滤波分析方法消除潮汐改正成分和零漂，对2008—2013年云南地区M≥55地震前3个月的资料进行处理，寻求同步异常特征，从时间尺度和空间尺度加以分析。同时利用倾斜固体潮观测模型求解潮汐响应率和线性漂移速率，从这些反应地壳介质的参数上，寻求震前的同步异常变化。
1数据选取及预处理
本文选用2008—2013年以来云南省内7次M≥55地震作为研究震例（表1，彝良为双震），选用震前3个月的水管倾斜资料整点值作为研究对象，并对数据进行挑选及预处理。
（1）对所研究的地震，逐个检查震前3个月水管倾斜整点值数据，要求记录连续完整，剔除缺数个数超过24小时（一天）的数据。
（2）查阅观测日志，对“十五”预处理数据库中的数据进一步筛选。剔除观测质量差、可信度低的数据。对无观测日志且数据变化不正常的也进行剔除。由于本研究是建立在预处理数据的基础上，数据的完整性、真实性对结果影响很大。因此在处理过程中，根据日志记录，只对明确由于停电、调仪器、归零、人为干扰产生观测误差的数据以及远震造成的干扰分不同情况进行了去台阶、擦除、补差值等处理。其余数据均不做改动，尽量保持数据真实性，并对每一個测项数据情况以及处理过程都做了详细记录，便于日后检查分析。
地震研究40卷第1期李智蓉等：临近地震前云南地区水管倾斜潮汐及非潮汐变化分析（3）经过前两步，可得减组合滤波的研究数据；在Nakai拟合中，考虑到形变观测受气压影响较大，对挑选出的数据，还要进一步挑选出气压正常且连续的水管倾斜观测资料作为最终拟合的研究样本。
2理论模型与计算方法
21减组合滤波分析方法对形变固体潮观测的整点值数据序列{y（j）}（j=1，2，…，n），采用下列减组合滤波（蒋骏等，2000）公式求取中心时刻的滤波值：
式中，Z（i）=y（i）-y（-i）为减组合。
减组合滤波能同时消除潮汐改正和零漂，最终得到高频信息，运用此方法对较短时段的观测资料处理可有效地识别出观测资料中存在的高频异常。
22Nakai拟合模型
倾斜固体潮理论值NS向分量ξ和EW向分量η可用下列公式计算（郗钦文，1982；张晶等，2003），下标s代表太阳的作用，未标注代表月亮的作用：
其中，ρ/ρ0=1+0001 679 26cos2φ，θ为地心天顶距，φ为测点地理纬度，τ为地方时角，δ为赤纬，R为天体至地心的距离，C为天体至地心的平均距离。
在倾斜潮汐观测资料中，除潮汐成分外，还包含有长期漂移项及气压项的影响，观测值的数学模型如下（杨林章等，1994）：Yt=RV（t+Δt）+KPP（t+ΔtP）+K0+K1t+K2t2（4）式中，Yt为倾斜观测值，R为观测值相对于理论值V的响应率，Δt为观测值相对理论值的时间滞后，KP为气压系数，P为观测点上的气压值，ΔtP为倾斜相对气压的滞后时间，K0、K1、K2分别为常数项、线性漂移速率、非线性加速率。 对式
以两天为单位进行最小二乘拟合计算，可求得潮汐和非潮汐参数和每小时的潮汐残差值。
3结果分析
31减组合滤波结果分析（1）异常的提取：①减组合滤波后所得残差，规定变化幅度超过2倍均方差视为异常，出现一次异常记为1；②逐日统计各测项异常数量之和，得到每日异常次数Ai，其中i为距离发震的时间；③将研究时段内所有异常次数累加得S；④考虑到异常出现的台站数（同步异常要求至少两个台站同时出现）和频次，研究规定当异常测项数大于等于3，且Ai≥2×（S/L）（异常频次日均值的2倍，其中L为窗长）时，认为在震前第i天出现同步异常。在统计同步异常过程中，未对异常方向作区分，不同台站各方向测项出现一次异常即记为1，因此所得结论中不能区分异常方向信息。表2列出了各台站在2008年8月21日盈江59级地震前出现同步异常时间段的统计结果。
（2）计算结果特征分析：7次震例同步异常分布时间、连续性各异，具体分布情况如图2所示。在日常跟踪中，连续性异常（异常持续时间超过3天）比较容易识别，因此本研究着重分析了连续性同步异常出现的时间段及频次，结果如表3所示。除中甸地震外，其余6次地震均出现连续性同步异常，持续时间为4～10 d不等，且出现的频次占各异常总频次的40%～100%；盈江M58、宁蒗M57、彝良M57、洱源M55地震的连续时段是重叠的，盈江M59、姚安M60地震的连续时段虽未重叠，但均分布在震前40～55 d左右。
（3）仅针对本次研究而言，连续性同步异常出现率达86%（6/7），且持续时间均在10 d以内，该方法可以尝试用于短临预测跟踪工作。
研究震例连续性同步异常时段/d占总异常频次的百分比震源机制类型序列类型2008-08-21盈江M5946～56625%（10/16）走滑前主余型2009-07-09姚安M6041～4550%（5/10）走滑主余型2011-03-10盈江M5871～79692%（9/13）走滑主余型2012-06-24宁蒗M5772～76100%（5/5）走滑主余型2012-09-07彝良M5724～2740%（4/10）走滑双震型2013-03-03洱源M5522～2656%（5/9）正断双震型2013-08-31中甸M59无0正断前主余型32Nakai拟合结果分析
（1）异常的提取：利用Nakai拟合模型，以2 d（48 h）为单位进行最小二乘拟合计算，可求得潮汐响应率、时间滞后、气压系数、气压滞后、线性漂移速率以及非线性加速率等参数，并求出每小时的潮汐残差值。由于模型中每项参数彼此之间是相关联的，因此文中挑选了潮汐响应率、线性漂移速率进行统计研究。
对于不同台站不同测项，规定所求得的各参数变化超出2倍均方差视为异常。针对每个震例，以2 d为单位，统计异常测项数目，将异常测项数N≥3视为有同步异常出现（保证至少有2个台站出现异常，即测项至少为3项）。
（2）计算结果特征分析：7次震例参数同步异常分布情况如图3所示，对于不同地震，潮汐响应率和线性漂移速率在震前3个月异常出现时段各异，持续时间短，在短临跟踪中难以识别。将潮汐观测值与拟合值作差（拟合精度各测项各时段均不相同，最小为096×10-3角秒，最大为46532×10-3角秒，其均值为2760×10-3角秒，具体不详述）可得每个测项潮汐残差时间序列图。残差计算结果中有一些大幅度的变化信息（图4a），这部分信息可能与仪器问题、远震或干扰信号等有关，根据每条残差曲线具体情况进行限幅后（图4b），分析残差成分，并未发现有用的规律。
4结论和讨论
利用减组合滤波方法提取高频成分，Nakai拟合求取潮汐和非潮汐参数，试图寻找短临阶段同步异常变化特征。由于数字化资料灵敏性高，易受周围各种已知和未知因素的干扰，且数据预处理情况会对结果造成一定影响，尤其是插值对高频成分的影响。因此在研究中为保证结果的可靠性，除了预处理过程需格外谨慎外，在异常提取中应选择至少3个测项（最少两个台站）同时段出现异常以保证异常的信度。经过分析研究，可得以下初步认识：
（1）减组合滤波后数据分析结果表明：所研究的7组地震中，除中甸地震外，其余6组地震均存在异常信息，连续性同步异常出现率达86%（6/7），异常出现在震前40～50 d之间。从机理上分析减组合滤波方法提取的是高频信息，表明这些55级以上地震，震前短期阶段有高频信息增多的现象出现，且这些高频信息可以被多台同时记录到。中甸地震前没有异常，可能与震中附近没有台站和高频信息传播的方向有关。从其他6次地震前得到的结果，可以认为，当有3个以上测项同时记录到持续的用减组合滤波方法提取的高频信息时，地震孕育可能进入了短临阶段，因此可以用该方法进行地震发震时间的短临跟踪。
（2）Nakai拟合结果在所选震例前未反映出有规律的异常变化，且潮汐残差值中夹杂一些大幅变化的干扰信息，因此目前該方法难以用于前兆短临异常的识别。
感谢中国地震局地震预测研究所张晶研究员的多次指导和刘琦提供的“钻孔应变处理程序V11”，感谢审稿专家的宝贵建议。
参考文献：
陈德福1993中国地震倾斜潮观测台网[J].内陆地震，1（3）：211-224
付虹，刘丽芳2003云南地区中短期前兆场异常与强震关系研究[J].地震研究，26（增刊1）：95-100
蒋骏，李胜乐，张雁滨，等2000地震前兆信息处理与软件系统（EIS2000）[M].北京：地震出版社
李杰，韩海华，马玉香，等2003数字化形变观测资料异常识别方法的应用及评价[J].东北地震研究，19（1）：25-33
刘序俨，张雁滨1991排除形变观测数据中降水干扰的数学物理方法的研究[J].地壳形变与地震，11（1）：36-40
刘仲全1997永胜台定点形变与丽江70级地震[J].地震研究，20（3）：273-277
刘仲全2001姚安65级地震前云南倾斜场变化研究[J].地震研究，24（4）：301-306
吴翼麟，周景明1981利用周期分析法研究观测资料中的干扰因素[J].地壳形变与地震，（3）：54-59
郗钦文，侯天航1986固体潮汐与引潮常数[J].中国地震，（2）：30-41
郗钦文1982固体潮理论值计算[J].地球物理学报，25（增刊1）：632-643
杨林章，何世海，郗钦文，等1994用潮汐体应变加卸载响应比研究岩石弹性性质的变化[J].中国地震，10（增刊1）：90-94
张晶，牛安福，高福旺，等2003数字化形变观测提取的地震短临异常特征[J].地震，23（1）：70-76
张晶，牛安福，高福旺，等2004新疆石河子、巴楚震前数字化形变短临异常[J]，地震，4（1）：82-87
HARRISON J C1976Cavity and topographic effects in tilt and strain measurement[J].JGR，81（2）：319-328
关键词：减组合滤波；Nakai拟合；水管倾斜；短临异常
中图分类号：P31563文献标识码：A文章编号：1000-0666（2017）01-0082-06
0引言
云南是我国破坏性地震频发的区域之一，仅2014年就发生3次M≥6地震。随着数字化形变观测的发展，研究人员可以从整点值、分钟值提取更多的短临前兆异常信息（刘仲全，1997，2001；付虹，刘丽芳，2003；Harrison，1976；刘序俨，张雁滨，1991；李杰等，2003；张晶等，2004），本文尝试寻找水管倾斜形变资料短临阶段潮汐及非潮汐变化特征，并浅析其与地震的关系。
地壳形变是地震发生过程中最直接的伴随现象（陈德福，1993），它与地震的关系最为直接，因此历来受到国内外地震预报探索者们的高度重视，也常常作为地震预测的重要判据之一，尤其在云南，地壳形变作为中短期预测依据已有很多应用实例，比仅用于趋势异常判定的模拟观测资料更具现实意义。识别形变观测资料异常的难点在于排除干扰，有的干扰能定量排除，有的只能定性排除，对于这些干扰因素规律性的探索及其观测值之间函数关系的研究，引起了许多学者的关注。
由于异常识别问题过于复杂，各个观测点表现出明显的地区性差异，因此这些研究未能获得显著的进展（吴翼麟，周景明，1981）。为了降低干扰排除的难度，对几个台站同步异常的提取同时进行观测分析是比较简便并且可信的方法，可排除局部环境所造成的影响，增加异常的可信度，同时假设异常信息来自震源体，在多个台站同步出现异常信号也容易从机理上得到解释。
目前，云南省共建有形变观测台站22个，其中，有水管倾斜仪的台站有11个，最南端的勐腊台因为距离研究地震的样本较远，将其剔除，不参与分析。多年的观测实践表明，水管倾斜仪记录的资料连续性、稳定性较好，能记录到清晰固体潮及年变趋势。而固体潮是唯一能准确计算出理论值的干扰（郗钦文，1982；郗钦文，侯天航，1986），对于已知的固体潮，可以采用均值法、别尔采夫滤波和组合滤波等方法剔除潮汐成分，提取高频信息。也可以利用固体潮理论值对台站观测值进行拟合，得到反映地壳介质的特征参量，通过对其分析研究，有可能了解到地壳局部岩石弹性性质的变化。
本文主要利用观测质量相对较好的水管倾斜仪资料，采用减组合滤波分析方法消除潮汐改正成分和零漂，对2008—2013年云南地区M≥55地震前3个月的资料进行处理，寻求同步异常特征，从时间尺度和空间尺度加以分析。同时利用倾斜固体潮观测模型求解潮汐响应率和线性漂移速率，从这些反应地壳介质的参数上，寻求震前的同步异常变化。
1数据选取及预处理
本文选用2008—2013年以来云南省内7次M≥55地震作为研究震例（表1，彝良为双震），选用震前3个月的水管倾斜资料整点值作为研究对象，并对数据进行挑选及预处理。
（1）对所研究的地震，逐个检查震前3个月水管倾斜整点值数据，要求记录连续完整，剔除缺数个数超过24小时（一天）的数据。
（2）查阅观测日志，对“十五”预处理数据库中的数据进一步筛选。剔除观测质量差、可信度低的数据。对无观测日志且数据变化不正常的也进行剔除。由于本研究是建立在预处理数据的基础上，数据的完整性、真实性对结果影响很大。因此在处理过程中，根据日志记录，只对明确由于停电、调仪器、归零、人为干扰产生观测误差的数据以及远震造成的干扰分不同情况进行了去台阶、擦除、补差值等处理。其余数据均不做改动，尽量保持数据真实性，并对每一個测项数据情况以及处理过程都做了详细记录，便于日后检查分析。
地震研究40卷第1期李智蓉等：临近地震前云南地区水管倾斜潮汐及非潮汐变化分析（3）经过前两步，可得减组合滤波的研究数据；在Nakai拟合中，考虑到形变观测受气压影响较大，对挑选出的数据，还要进一步挑选出气压正常且连续的水管倾斜观测资料作为最终拟合的研究样本。
2理论模型与计算方法
21减组合滤波分析方法对形变固体潮观测的整点值数据序列{y（j）}（j=1，2，…，n），采用下列减组合滤波（蒋骏等，2000）公式求取中心时刻的滤波值：
式中，Z（i）=y（i）-y（-i）为减组合。
减组合滤波能同时消除潮汐改正和零漂，最终得到高频信息，运用此方法对较短时段的观测资料处理可有效地识别出观测资料中存在的高频异常。
22Nakai拟合模型
倾斜固体潮理论值NS向分量ξ和EW向分量η可用下列公式计算（郗钦文，1982；张晶等，2003），下标s代表太阳的作用，未标注代表月亮的作用：
其中，ρ/ρ0=1+0001 679 26cos2φ，θ为地心天顶距，φ为测点地理纬度，τ为地方时角，δ为赤纬，R为天体至地心的距离，C为天体至地心的平均距离。
在倾斜潮汐观测资料中，除潮汐成分外，还包含有长期漂移项及气压项的影响，观测值的数学模型如下（杨林章等，1994）：Yt=RV（t+Δt）+KPP（t+ΔtP）+K0+K1t+K2t2（4）式中，Yt为倾斜观测值，R为观测值相对于理论值V的响应率，Δt为观测值相对理论值的时间滞后，KP为气压系数，P为观测点上的气压值，ΔtP为倾斜相对气压的滞后时间，K0、K1、K2分别为常数项、线性漂移速率、非线性加速率。 对式
以两天为单位进行最小二乘拟合计算，可求得潮汐和非潮汐参数和每小时的潮汐残差值。
3结果分析
31减组合滤波结果分析（1）异常的提取：①减组合滤波后所得残差，规定变化幅度超过2倍均方差视为异常，出现一次异常记为1；②逐日统计各测项异常数量之和，得到每日异常次数Ai，其中i为距离发震的时间；③将研究时段内所有异常次数累加得S；④考虑到异常出现的台站数（同步异常要求至少两个台站同时出现）和频次，研究规定当异常测项数大于等于3，且Ai≥2×（S/L）（异常频次日均值的2倍，其中L为窗长）时，认为在震前第i天出现同步异常。在统计同步异常过程中，未对异常方向作区分，不同台站各方向测项出现一次异常即记为1，因此所得结论中不能区分异常方向信息。表2列出了各台站在2008年8月21日盈江59级地震前出现同步异常时间段的统计结果。
（2）计算结果特征分析：7次震例同步异常分布时间、连续性各异，具体分布情况如图2所示。在日常跟踪中，连续性异常（异常持续时间超过3天）比较容易识别，因此本研究着重分析了连续性同步异常出现的时间段及频次，结果如表3所示。除中甸地震外，其余6次地震均出现连续性同步异常，持续时间为4～10 d不等，且出现的频次占各异常总频次的40%～100%；盈江M58、宁蒗M57、彝良M57、洱源M55地震的连续时段是重叠的，盈江M59、姚安M60地震的连续时段虽未重叠，但均分布在震前40～55 d左右。
（3）仅针对本次研究而言，连续性同步异常出现率达86%（6/7），且持续时间均在10 d以内，该方法可以尝试用于短临预测跟踪工作。
研究震例连续性同步异常时段/d占总异常频次的百分比震源机制类型序列类型2008-08-21盈江M5946～56625%（10/16）走滑前主余型2009-07-09姚安M6041～4550%（5/10）走滑主余型2011-03-10盈江M5871～79692%（9/13）走滑主余型2012-06-24宁蒗M5772～76100%（5/5）走滑主余型2012-09-07彝良M5724～2740%（4/10）走滑双震型2013-03-03洱源M5522～2656%（5/9）正断双震型2013-08-31中甸M59无0正断前主余型32Nakai拟合结果分析
（1）异常的提取：利用Nakai拟合模型，以2 d（48 h）为单位进行最小二乘拟合计算，可求得潮汐响应率、时间滞后、气压系数、气压滞后、线性漂移速率以及非线性加速率等参数，并求出每小时的潮汐残差值。由于模型中每项参数彼此之间是相关联的，因此文中挑选了潮汐响应率、线性漂移速率进行统计研究。
对于不同台站不同测项，规定所求得的各参数变化超出2倍均方差视为异常。针对每个震例，以2 d为单位，统计异常测项数目，将异常测项数N≥3视为有同步异常出现（保证至少有2个台站出现异常，即测项至少为3项）。
（2）计算结果特征分析：7次震例参数同步异常分布情况如图3所示，对于不同地震，潮汐响应率和线性漂移速率在震前3个月异常出现时段各异，持续时间短，在短临跟踪中难以识别。将潮汐观测值与拟合值作差（拟合精度各测项各时段均不相同，最小为096×10-3角秒，最大为46532×10-3角秒，其均值为2760×10-3角秒，具体不详述）可得每个测项潮汐残差时间序列图。残差计算结果中有一些大幅度的变化信息（图4a），这部分信息可能与仪器问题、远震或干扰信号等有关，根据每条残差曲线具体情况进行限幅后（图4b），分析残差成分，并未发现有用的规律。
4结论和讨论
利用减组合滤波方法提取高频成分，Nakai拟合求取潮汐和非潮汐参数，试图寻找短临阶段同步异常变化特征。由于数字化资料灵敏性高，易受周围各种已知和未知因素的干扰，且数据预处理情况会对结果造成一定影响，尤其是插值对高频成分的影响。因此在研究中为保证结果的可靠性，除了预处理过程需格外谨慎外，在异常提取中应选择至少3个测项（最少两个台站）同时段出现异常以保证异常的信度。经过分析研究，可得以下初步认识：
（1）减组合滤波后数据分析结果表明：所研究的7组地震中，除中甸地震外，其余6组地震均存在异常信息，连续性同步异常出现率达86%（6/7），异常出现在震前40～50 d之间。从机理上分析减组合滤波方法提取的是高频信息，表明这些55级以上地震，震前短期阶段有高频信息增多的现象出现，且这些高频信息可以被多台同时记录到。中甸地震前没有异常，可能与震中附近没有台站和高频信息传播的方向有关。从其他6次地震前得到的结果，可以认为，当有3个以上测项同时记录到持续的用减组合滤波方法提取的高频信息时，地震孕育可能进入了短临阶段，因此可以用该方法进行地震发震时间的短临跟踪。
（2）Nakai拟合结果在所选震例前未反映出有规律的异常变化，且潮汐残差值中夹杂一些大幅变化的干扰信息，因此目前該方法难以用于前兆短临异常的识别。
感谢中国地震局地震预测研究所张晶研究员的多次指导和刘琦提供的“钻孔应变处理程序V11”，感谢审稿专家的宝贵建议。
参考文献：
陈德福1993中国地震倾斜潮观测台网[J].内陆地震，1（3）：211-224
付虹，刘丽芳2003云南地区中短期前兆场异常与强震关系研究[J].地震研究，26（增刊1）：95-100
蒋骏，李胜乐，张雁滨，等2000地震前兆信息处理与软件系统（EIS2000）[M].北京：地震出版社
李杰，韩海华，马玉香，等2003数字化形变观测资料异常识别方法的应用及评价[J].东北地震研究，19（1）：25-33
刘序俨，张雁滨1991排除形变观测数据中降水干扰的数学物理方法的研究[J].地壳形变与地震，11（1）：36-40
刘仲全1997永胜台定点形变与丽江70级地震[J].地震研究，20（3）：273-277
刘仲全2001姚安65级地震前云南倾斜场变化研究[J].地震研究，24（4）：301-306
吴翼麟，周景明1981利用周期分析法研究观测资料中的干扰因素[J].地壳形变与地震，（3）：54-59
郗钦文，侯天航1986固体潮汐与引潮常数[J].中国地震，（2）：30-41
郗钦文1982固体潮理论值计算[J].地球物理学报，25（增刊1）：632-643
杨林章，何世海，郗钦文，等1994用潮汐体应变加卸载响应比研究岩石弹性性质的变化[J].中国地震，10（增刊1）：90-94
张晶，牛安福，高福旺，等2003数字化形变观测提取的地震短临异常特征[J].地震，23（1）：70-76
张晶，牛安福，高福旺，等2004新疆石河子、巴楚震前数字化形变短临异常[J]，地震，4（1）：82-87
HARRISON J C1976Cavity and topographic effects in tilt and strain measurement[J].JGR，81（2）：319-328