Ability of decreasing noise and the characteristics of nearsurface wave field around Dazhai borehole in Pu'er
-
摘要: 基于2011年建立的云南普洱大寨深井台站，开展了噪声压制及附近波场特征研究.通过计算该台站的噪声功率谱概率密度函数，显示该井下台站对1 Hz以上的高频噪声具有明显的压制效果, 最高能降低40 dB，其降噪能力优于其它井下台阵，推断与该台站附近的场地条件有关.基于地表与井下地震记录的差异，应用正则化反卷积干涉方法进一步研究该台站附近的波场特征.以地表记录为参考，对井下记录进行反卷积，获取两台站之间的格林函数，直接识别出了原始记录上无法区分的上行入射波与下行地表反射波，然后利用两震相的到时差建立了一个浅层地震波速度模型，与理论模拟的结果一致.研究结果表明，相对于地表观测，井下台站在压制噪声和近地表地震波传播特征研究等方面具有很大的优势，同时该研究对其它地区开展深井观测具有参考意义.Abstract: In 2011, the Dazhai borebole station was established in Pu'er, Yunnan Province, to carry out some researches on noise suppression and the wavefield characteristics. We compute the power spectral density of continuous noise recordings within ten months to quantify the reduction in seismic noise with the station depth. The results demonstrate that the downhole station has obvious effect on decreasing noise up to 40 dB within the frequency band more than 1 Hz, which might be associated with the local site conditions. The Green's functions are then obtained between the two receivers by deconvoluting the down-hole recordings with the surface ones. From the deconvolution wavefield, the arrivals of up-going incident and down-going surface-reflection waves are directly identified that cannot be distinguished from the original seismic recordings. The two phases are utilized to set up a shallow seismic wave velocity model, which is consistent with the theoretical results. A comparison with the surface observation suggests that the borehole stations have great advantages in reducing noise and studying the propagation characteristics of near-surface seismic waves. The research here is of great significance for future down-hole observations in other sites.
-
Keywords:
- borehole /
- noise reduction /
- regularized deconvolution /
- incident wave /
- surface-reflected wave
-
引言
随着人类社会城镇化和工业化的高速发展，各种交通、工业等城市噪音严重干扰了传统的地震观测结果，微震和来自地下的微弱地球物理信息混杂在高噪音背景下，无法被有效地识别和利用(Fukao, Ishibashi, 1996；徐纪人，赵志新，2009).大量研究表明，城市噪声随着观测台站深度的增加而迅速衰减，因此相比于地表观测，井下观测能够有效地压制噪音，得到信噪比较高的地震记录，是在高噪音区域开展地震观测的有效手段(Hauksson et al, 1987; Aster, Shearer，1991; Boese et al, 2015).
井下观测在美国、日本和我国台湾等地区被用于地震监测能力提升和科学研究(Asai et al, 2005; Okubo et al, 2005; Zoback et al, 2011; Ma et al, 2012)，而我国大陆作为一个地震多发区，井下台阵仍处于建设发展的初级阶段(李凤杰，1989；冯德益等，1990；徐纪人，赵志新，2009).井下观测多用于盆地地区，例如在首都圈布设了大量井下台站，但由于覆盖层较厚，导致台基多为沉积层，且井深较浅，多集中在两百米左右，提高信噪比的效果并不明显(王林瑛等，2008；张尉等，2009).
井下观测除了具有压制噪声的作用外，在研究近地表地震波传播特性方面也独具优势.在地表及地下不同深度布设地震计，即建立井下垂直台阵，以此探测不同深度介质对地震波传播的影响，即可为研究近地表介质速度特性、土层放大和衰减效应等方面提供直观的依据，这对于强地面震动和地震灾害研究具有重要的意义(Fukushima et al, 1992; Parolai et al, 2009; Nakata, Snieder, 2012; Nakata, 2013; Fukushima et al, 2016).但是目前国内对利用井下数据进行浅层地震波传播特性的研究也很少(Chong, Ni, 2009；刘渊源等，2011).同时，由于缺少同一地点地表台站，只能联合附近其它地点的地表台站进行结构反演等研究，而无法精确地探讨井下台阵对于降噪以及近地表介质特性对地震波传播的影响等方面的效果.
为深入研究井下观测对降噪的作用及基岩地区井中观测的效果，我们于2011年在云南普洱大寨布设一口深约410 m的井，并在地表和井深375 m处分别配备了一套短周期地震计以对比地表与井下观测资料的差异.本文拟以大寨井下台站观测数据作为分析对象，展开噪声水平和井下方位角校正等研究，然后通过反卷积干涉方法重建不同接收点之间相对简单的格林函数，从而直接识别出一些无法在原始波场上区分的震相信息，以期获取该地区近地表地震波传播特征及地下介质特性，为未来其它垂直井下台阵研究提供参考依据.
1. 观测系统
兰坪—思茅盆地地处云南省西部, 呈北西—南东向展布, 南宽北窄, 是西南“三江”地区的大型中新生代盆地.由于受到印度板块与欧亚大陆之间的碰撞和持续挤压作用，该地体内部形成了一系列逆冲和走滑断裂带，北起兰坪，南至景洪，普遍发育轴向近北西或北北西向的褶皱带；地体外部则发育多条大规模左旋走滑断裂带(仝亚博等，2014).大寨井下台站布设于该盆地南部的普洱地区(图 1)，地表和井下375 m处分别配置了同一型号短周期FSS-3DBH地震计，响应频带范围为0.5—50 Hz，两个深度的观测记录由同一6通道24位数据的EDAS-24IP采集器采集.
![]()
图 1 大寨深井台站位置及其记录到的地震分布黑色圆圈为台站附近的地震分布，黑色实心圆圈为文中反卷积计算用到的地震事件Figure 1. Location of the Dazhai borehole station (triangle) and distribution of local earthquakes Circles represent the location of local earthquakes and solid circles are the earthquakes used in this paper.F1 : Wuliangshan fault; F2 : Lancangjiang fault; F3 : Honghe fault大寨井孔附近的地层结构如图 2所示.可见：近地表有22.43 m的黏土覆盖层，下覆岩层为泥岩和砂岩.该台站自2011年9月起投入运行，记录了多个地震事件.本文选取2011年9月至2012年10月期间的数据，这期间的数据连续性较好，且受干扰较小, 井孔附近600 km范围内共发生96个地震事件，震级集中在1.5≤ML < 4.0和4.0≤MS < 6.0(图 1).
![]()
图 2 大寨井孔结构两个黑色三角形分别代表位于地表和井下的地震计Figure 2. Dazhai borehole columnar diagramSolid triangles indicate sensors at different depths2. 井下地震计方向校正
井下地震计安置的关键问题之一是保证安装方位的准确性，因此数据分析之前，需对其进行方位角校正.对于垂直井下台阵而言，需要1个能够精确定位的地表地震计，基于地表与井下记录的相似性原理，以地表地震计方位作为标准，展开两者之间的相关性分析，从而实现对井下地震计的方位校正(Seale, Archuleta, 1989; Parolai et al, 2009; Grigoli et al, 2012).本研究所用的两个地震计为同一厂家、同一型号设备，便于开展相关性分析.
利用地表台站校正井下台站方位角的一个基本假设是两个台站记录的波形是相似的.由于不同深度地震计记录到的噪声差异较大，使用高信噪比的地震事件获取的结果比利用噪声数据进行计算更加可靠(吕永清等，2007)；而且，地震从地下一定深度传至地表的过程中，某些频段尤其是高频段的地震波会发生散射衰减，从而造成井下与地表记录在高频段的相关性较低，因此需进行低通滤波以去除显著差异.对于滤波频带的选取，目前并没有统一的标准.谢剑波(2014)针对所选取的6个远震事件的数据，将频谱分为两种情况进行处理：对于震中距相对较小、频率成分稍高的云南地震，采用0.1—0.5 Hz通带；对于其它5个事件，采用0.06—0.2 Hz通带.吕永清等(2007)在校正方位角之前，对地震事件进行了0.05—10 Hz的带通滤波处理.因此，为了确保校正结果的稳定性，本文从96个地震事件中选取信噪比较高的56个事件，去除仪器响应后计算这56个地震事件的井下与地表水平分量的标准傅里叶谱比，结果如图 3所示.短周期FSS-3DBH地震计的频率特性曲线如图 4所示，去除仪器响应后，能够校正仪器原来的低频振幅压制效果和相位，从而达到一定程度的扩频效果.李少睿等(2016)利用与本文相同的短周期地震计进行研究，认为该仪器对0.1 Hz的信号仍然有很好的响应, 因此本文选取频率大于0.1 Hz的信号进行研究.从图 4可以看出，曲线峰值对应的频率值为2.4 Hz, 此峰值可能是由于上行波和下行波传播过程中发生干涉所致.由于本文的地震震级相对较小，且频率成分较高，因此低通滤波需要以2.4 Hz或者更低频率作为截止频率.因此，本文分别选取2.0 Hz, 1.0 Hz作为截止频率进行滤波，以2011年11月2日发生在距离台站71.21 km的ML3.7地震事件为例，低通滤波结果显示1.0 Hz以下的地表与井下记录相似度较高(图 5)，这可能与图 4中所显示的振幅谱比在1.0 Hz以下几乎无变化有关.因此，进行校正前对各事件进行1.0 Hz的低通滤波.
![]()
图 3 所选取56个地震事件的地表记录与井下记录南北分量(a)和东西分量(b)的傅里叶谱比Figure 3. Fourier spectral ratio of local earthquakes using NS (a) and EW (b) components between surface and borehole recordings![]()
图 4 短周期FSS-3DBH地震计的频率特性曲线Figure 4. Frequency characteristic curves of the short-period seismometer FSS-3DBH![]()
图 5 ML3.7地震事件南北向分量的原始波形(a)及其2.0 Hz(b)和1.0 Hz(c)低通滤波结果Figure 5. Original waveforms (a) and low-pass filtered waveforms with the cut-off frequencies of 2.0 Hz (b) and 1.0 Hz (c) of the NS component of the ML3.7 eventBlack lines indicate the borehole recordings, and grey lines indicate the surface recordings假设井下地震计的角度偏差为α，则可以通过
(1) 进行校正而获得真实的水平向记录, 式中BNS和BEW分别表示校正前井下的南北分量和东西分量, α表示旋转的角度(自北向东旋转为正)，TNS和TEW分别表示旋转后井下的南北分量和东西分量.
旋转后的记录应该与地表记录相应方向的记录之间的差异最小.本文中分别采用均方根误差和相关系数描述波形之间的差异，其中均方根误差定义为
(2) 式中SNS和SEW分别表示地表南北分量和东西分量.
互相关函数为
(3) 式中T和S分别表示旋转后的井下和地表同一分量，S*是S的共轭，τ为两者之间的延迟时间.
为得到井下地震计的真实方位角，我们按照式(1)将井下地震记录在-180°—180°范围内以一定增量(此处为1°)进行旋转，然后分别按照式(2)和式(3)计算旋转后的记录与地表记录之间的均方根误差和互相关函数，并将均方根误差最小和互相关函数最大的角度作为井下地震计的方位角.
利用不同事件获取的最佳旋转角度分布情况如图 6所示.通常中值和均值最接近真实值，且两者数值非常接近，因此此处选用中值和均值的平均值即-6°, 作为井下地震计水平分量应该旋转的最佳角度(表 1).
![]()
图 6 基于均方根误差方法(a)和互相关函数方法(b)得到的不同地震事件南北分裂的最佳旋转角度统计结果黑色圆圈表示不同事件对应的角度，虚线表示不同角度的中值及均值Figure 6. The optimal rotation angles corresponding to NS components of different events by RMS error (a) and cross-correlation (b)Open circles indicate different angles corresponding to different events, and dashed lines indicate the median and mean value of all rotation angles表 1 最佳旋转角度的统计结果Table 1. The statistics of optimal rotation angles方法 均值 中值 两者平均值 均方根误差法 -6.66° -7° -6.83° 互相关函数法 -5.36° -5° -5.18° 图 7给出了ML3.7地震事件地表及井下南北分量旋转前后的波形对比, 可见旋转后井下与地表记录的相似性稍有提高，但变化不大.这说明本文所使用的井下地震计在安装时，其方位角误差较小. Wang等(2016)关于中国东北地震台阵方位角的研究结果表明, 10°的方位角误差最多能够引起5%的振幅变化，而本文研究的井下台站方位角误差很小，因此后面的分析中我们不再对井下地震记录进行特殊处理.
![]()
图 7 ML3.7地震事件井下南北分量按照最佳旋转角度旋转前(a)、后(b)的波形对比Figure 7. Waveforms before (a) and after (b) rotating the NS component of the ML3.7 event according to the optimal alignment angleBlack lines indicate the borehole recordings, and grey lines indicate the surface recordings3. 噪声分析
通过计算功率谱密度(power spectral density, 简写为PSD)可以得到台站的背景噪声水平. Peterson(1993)定量地分析了全球各地75个固定地震台站的地震背景噪声功率谱密度，建立了全球地震噪声模型，包括新高噪声模型(new high noise model，简写为NHNM)和新低噪声模型(new low noise model，简写为NLNM).为全面地了解大寨井下地震台站的降噪效果，我们利用2012年1—10月的连续波形数据计算地表和井下台站不同时间的噪声功率谱密度，并在此基础上获得了台站噪声的概率密度函数.
本文参照葛洪魁等(2013)对山西临汾地区流动台站的背景噪声加速度功率谱密度的计算方法，将各台站水平分量和垂直分量的连续地震波形记录截成长度为1小时(3600 s)的数据段，相邻数据段重合50%.这样每个台站的每一分量每天有47个数据段，10个月的数据被分成约1万4000个数据段.对去除仪器响应后的1小时连续数据进行如下处理：
1) 为尽可能地减小PSD方差，将1小时的数据段分成14小段，每小段数据的长度为1000 s，设定步长为200 s；
2) 将Welch窗口叠加到数据段上，以压制数据序列进行快速傅里叶变换所产生的旁瓣效应，平滑快速傅里叶变换的计算结果；
3) 计算时间序列的速度功率谱密度；
4) 重复步骤(2)—(3), 计算所有14小段的速度功率谱密度，取平均值，得到1小时长度数据段的速度功率谱密度随频率的分布；
5) 将速度功率谱密度转化为加速度功率谱密度后，分别对各个台站的结果进行1/8倍频程滤波，得到平滑的且在对数坐标上均匀分布的加速度功率谱密度，进一步将其变换为以dB为单位，以便于与全球地震噪声模型(Peterson, 1993)进行对比分析.
上述计算的是单条1小时数据的加速度功率谱密度随周期的分布.但是，由于外界环境的变化和突发事件(如地震等)的影响，功率谱密度将随时间不断变化.功率谱谱密度在不同周期下某一时间段内取某值的概率，即概率密度函数(probability density function，简写为PDF), 能够较直观地表现出地震观测台站的噪声水平，已经成功地应用于诸多台阵观测质量的评估(McNamara, Boaz, 2005; Díaz et al, 2010; 葛洪魁等, 2013; Boese et al, 2015).该方法在计算过程中不需要专门挑选连续平静的记录，而是将诸如地震等各种扰动统一进行处理，能够全面地反映台站噪声水平的动态变化范围.
我们对1/8倍频程滤波后的加速度功率谱密度进行统计.设T为统计的中心周期，d表示任意PSD值(单位为dB)，以1 dB为间隔进行统计，则(T，d)对应的概率密度函数为
(4) 式中，NTd是以T为中心周期的所有功率谱密度个数，Nd是周期T处在[d，d+1]范围内的功率谱密度个数.
统计所有的中心周期，得到功率谱密度在-190—-85 dB范围内的分布.图 8给出了该台站垂直和水平分量的加速度功率谱概率密度函数.噪声水平的统计特征通常用众数、均值或中值表示，其中：众数是功率谱密度最可能的取值，最能反映台站的平均噪声水平，但是在5 Hz以上的频段由于受到人为噪声规律变化的影响容易出现跳变；平均值容易受异常值的影响，而严重偏离众数统计；中值统计曲线平滑，基本与众数统计曲线重合，且不易受系统突跳的影响.因此，一般取中值统计代表台站相应分量的平均噪声水平(McNamara, Boaz，2005；葛洪魁等，2013；Boese et al, 2015).
![]()
图 8 井下(a)和地表(b)台站三分量的概率密度函数Figure 8. Probability density function of three components from borehole (a) and surface (b) station以PDF的中值作为台站的平均噪声水平，分析两个台站不同分量的噪声水平差值随频率的变化情况(图 9).由于台站的地震计记录周期范围为0.5—50 Hz, 而本文对噪声数据去除了仪器响应，所以此处考虑仪器频带范围附近有效频带(0.1—20 Hz)内的噪声情况.可以看出:相较于地表台站，当频率大于1 Hz时，井下台站具有显著的降噪效果；随着频率增加，降噪效果日趋显著；当频率大于10 Hz时，降噪效果随着频率增加而逐渐降低；而频率小于1 Hz时，两台站的噪声水平几乎相当.因此，大寨井下台对于降低频率大于1 Hz的噪声具有良好的效果, 对于5 Hz附近的噪声压制效果最好，可降低噪声40 dB之多.
![]()
图 9 地表和井下台站不同分量的概率密度函数中值差值Figure 9. Difference of the median value between probability density function of three-components recorded by surface and borehole stations4. 井孔附近介质结构
4.1 反卷积方法
傅里叶谱比和反卷积干涉是两种最常用的用于分析井下垂直台阵数据的方法.前者多用于研究场地效应，而反卷积方法能够针对复杂的输入波场，通过重建不同接收点之间相对简单的格林函数，从而直接识别出一些在原始波场上难以区分的震相信息，以此获取地下介质特性参数等.因此，本文利用反卷积干涉法展开井下与地表之间的地下结构研究.
不同深度z1和z2的地震记录之间的频率域反卷积表示为
(5) 式中，u(z1, w)和u(z2, w)分别表示深度z1和z2处地震记录的傅里叶谱，w表示角频率.
将式(5)所获取的反卷积波场变换至时间域后，能够识别出在垂直分布的不同台站之间传播的上行波和下行波(Fukushima et al, 1992; Snieder, Safak, 2006; Nakata et al, 2014; Fukushima et al, 2016).地震P波或S波从地下一定深度入射至井孔后，上行传至地表，即为上行的入射波，在地表发生反射，重新传至井孔中，即为下行的地表反射波.在反卷积波场上，入射波应出现在负时间轴上，而地表反射波出现在与入射波到时对称的正时间轴上，即为从地表反射至井孔内的下行波.
通过上行波和下行波的旅行时间和射线传播路径，即可获取井下与地表之间的平均波速剖面.由于射线路径的复杂性，选取近似垂直入射的地震事件，将传播模型简化，根据台站之间的垂直距离及上行、下行波的旅行时间，便可估算出不同深度台站之间的平均速度结构.
由于井下与地表记录傅里叶谱差异，式(5)无法获取稳定的反卷积波场，因此，一种规则化参数的反卷积方法被发展起来(Tikhonov，Arsenin, 1977; Bertero, Boccacci, 1998; Mehta et al, 2007a, b; Parolai et al, 2009; Nakata, Snieder, 2012)：
(6) 
(7) 式中，ε是一种规则化参数，不同的值代表z2深度处的地震计记录到的能量谱均值所占的不同比例，该参数会对结果的分辨率产生重要的影响.
另外还有一种基于Landweber迭代理论的反卷积算法也可以实现对反卷积格林函数的约束，从而获得稳定的波场(Bertero, Boccacci, 1998；Parolai et al, 2009)：
(8) 
(9) 
(10) 式中，Wn(w)为Landweber因子，n为迭代次数, τ为约束因子，对噪声有一定的约束作用.按照Bindi等(2000)和Parolai等(2009)的设定，取τ=1/|umax(z2, w)|2.在该约束因子一定的情况下，随着迭代次数的增加，计算结果的精度也随之提高，但所受到噪声的影响也会随之逐渐增大.为了使计算结果稳定，迭代次数的选择尤为关键.
4.2 反卷积结果分析
考虑到方位角对结果可能会产生一定的影响，因此在进行计算前，首先将每个地震事件的水平向记录旋转成沿事件—台站分布的径向和切向分量.对地震记录进行0.5 Hz的高通滤波后，分别选取包含P波和S波到时的时间窗口，以地表记录作为参考，对井下记录进行反卷积计算.
为了寻找最优化的正规化迭代次数，以表 2中的2号地震事件为例，比较不同正则化参数和迭代次数下的反卷积波场.正则化参数ε分别取地表记录平均能量谱的0.05%, 0.1%，0.5%, 1%, 5%, 10%, 20%, 50%, 迭代次数n自10逐渐变化至10000.图 10给出了该地震水平分量的反卷积波场，可见，随着ε的降低或者n的增加，所获取波场的分辨率提高, 但同时也受到更多噪声的干扰.比较结果显示，取n为1000，ε设为地表记录平均能量谱的10%，即可获得相对精度较高且受到噪声干扰较小的稳定反卷积波场.
表 2 用于反卷积计算的事件列表Table 2. List of the events used for deconvolution事件序号 发震时间 ML 震中距/km 年-月-日 时：分：秒 1 2011-11-02 19:33:05.1 3.7 71.21 2 2011-11-09 09:04:23.4 3.2 124.63 3 2011-11-16 23:47:37.5 2.6 121.54 4 2011-12-28 15:16:32.3 3.5 90.64 5 2012-02-22 00:41:39.7 2.3 145.54 6 2012-03-12 01:24:19.1 2.4 68.49 7 2012-04-13 23:06:45.7 2.7 78.48 8 2012-04-29 01:33:01.3 2.7 154.29 9 2012-07-28 08:16:15.5 3.2 57.63 10 2012-07-30 00:05:32.5 4.2(MS) 41.47 注：地震目录引自中国地震台网中心(2015). ![]()
图 10 2011年11月9日ML3.2地震水平分量反卷积波场(a)不同正则化参数ε下的反卷积波场, 线条颜色由浅变深，表示ε逐渐减小；
(b)不同迭代次数n下的反卷积波场, 线条颜色由深变浅，表示n逐渐增大Figure 10. Deconvolution wavefield of the event ML3.2 with different regularization parameters(a) Deconvolution variation with regularization parameter ε decreasing from light to dark lines;
(b) Deconvolution variation with Landweber iterations n increasing from dark to light lines对比两种反卷积方法可知，除了反卷积波场的频率成分稍有不同以外，所获取的直达波及地表反射波结果基本一致，因此，本文随后只基于正则化参数ε的反卷积方法进行计算即可，并取ε为地表记录平均能量谱的10%.图 11给出了表 2中事件2的反卷积波场结果，如虚线所标示，可以识别对称的入射波和地表反射波，同时，也有其它峰值出现.图 12给出了文中所用地震事件的径向、切向分量及垂直分量的反卷积结果，可见：在垂直分量的反卷积波场上，入射波与地表反射波是一对时间对称、振幅相当的峰值；在水平分量的反卷积波场上，入射波较易识别，而地表反射波相对较弱，甚至无法明显识别.这种现象与所选地震的震级大小、震中距及方位角并无直接的关联，应与地下介质特性有关.因此，接下来对该井口周围的地下介质特性展开分析.
![]()
图 11 表 2中事件2(ML3.2)的原始记录(a)及其相应的反卷积波场(b)红线表示用于进行反卷积计算的时间窗口，黑色虚线是井下与地表之间相互对称的上行波和下行波Figure 11. The deconvolution wavefield (b) from the original records (a) of No.2 event listed in Table 2The red lines indicate the time window used for deconvolution, and the dashed lines show the up-going wave and down-going wave from bottom to surface![]()
图 12 表 2中地震事件的原始井下记录(a)及其对应的反卷积波场(b)(a)中红色实线表示用于反卷积计算的时间窗口, (b)中红色虚线分别表示入射波和地表反射波的位置Figure 12. Aligned deconvolution wavefield (b) from the original borehole records (a) of all the events listed in Table 2Red lines indicate the time window used for deconvolution, and vertical dashed red lines indicate the arrivals of incident and surface-reflected waves4.3 等效波速模型
参考测井所得的地质资料(图 2)，依据各地层岩性及埋深的不同，首先估算一个初始的剪切波速度模型范围，如图 13a所示(《工程地质手册》编委会, 2007).随着埋深的增加，剪切波速度逐渐增加, 其中:顶部的黏土层速度较低，处于150—300 m/s的范围内，与下面的沙质泥岩岩层存在一个明显的分界；第二层的沙质泥岩波速范围为500—1000 m/s，与第三层的石英砂岩也存在明显分界.总体来说，按照剪切波速的差异，该井孔的结构大致可分为3层，包括黏土层、沙质泥岩层、泥砂岩层，因此会有两个明显的波阻抗界面，此外地下介质中还存在多个小的波阻抗界面.
![]()
图 13 大寨井孔附近的S波(a)和P波(b)速度模型实心三角形表示两个地震计的深度Figure 13. Shear velocity model (a) and P wave velocity (b) around the boreholeTriangle indicates sensors at different depths, and dashed lines indicate the estimated S velocity本文用的井下台站属于浅层井孔，仅375 m深, 而选用的地震事件发生深度均大于10 km，且距离台阵不超过150 km, 相对于地震发生的层位，井下台站所处的地层速度非常低.因此，由斯奈尔(Snell)定律可以将入射至该台站的射线路径视为垂直入射(Fukushima et al，1992；Chong，Ni, 2009).垂直分量反卷积波场的上行波和下行波到时分别为-0.1 s和0.1 s, 而水平分量反卷积波场的上行波和下行波到时分别为-0.34 s和0.34 s，则可获取地表至井下地震计所在深度之间的平均剪切波速度(图 13).对于P波，整个地层的平均速度为3750 m/s；对于S波，平均速度为1100 m/s.
为了验证所获取速度模型的准确性，利用Zhu和Rivera(2002)的频率-波数法计算包括上下行地震波的理论地震图，与所获取的反卷积格林函数中上下行波的到时差进行对比.该方法是基于水平分层模型合成理论地震图，计算时对频率和波数分别进行积分，采用传播矩阵计算地震的位移场分布.以井下至地表作为第一层，结合全球地壳速度模型Crust1.0，本文给出了计算理论地震图所用的分层模型，如表 3所示.
表 3 合成理论地震图所用到的地壳模型Table 3. Velocity model for calculating the synthetic seismograms层数 vP/(km·s-1) vS/(km·s-1) 层厚/km QP QS 1 3.75 1.10 0.375 1000 500 2 6.10 3.55 13.85 1000 500 3 6.30 3.65 12.34 1000 500 注：QP，QS表示P波和S波的地下介质品质因子. 基于表 3给出的地壳模型，针对表 2中事件1(ML3.7, 震中距为71.21 km)，以切向分量作为分析对象，计算了台站记录到的0, 100, 200, 300, 400, 500和600 m深度处的理论地震图，如图 14所示.可以看出：100 m以上深度的台站记录并不能很好地区分开直达S波与地表反射波；对于100 m至600 m深度的台站记录，随着井下台站深度的增加，直达波与地表反射波更加容易区分.理论地震图(图 14a)与实际地震图(图 14b)的S波到时基本一致，后者由于实际场地的复杂性，无法从中直接分离出直达波和反射波震相.而从理论地震图得到的两震相到时差与本文通过反卷积计算得到的结果(图 14c)一致，说明本文的速度模型基本可靠.但是，由于缺少中间深度的地震计，很难从反卷积波场上明确其它峰值所对应界面的反射波，无法对井下与地表之间的地层进行更精细的分层.
![]()
5. 讨论与结论
云南普洱大寨井下台站是一个地表及井下联合观测的典型范例，对其观测资料的研究充分证明了井下台站在降低噪声、研究近地表土层介质特性、地震波的传播特征等方面均具有显著的效果，是城市地震监测系统发展的方向，同时，也为未来布设其它井下台阵提供了参考.本文研究分析主要得到以下结论：
1) 基于井下与地表地震计记录的相关性原理，以地表地震计作为参考，对井下地震计的方位角进行了校正，结果表明在垂直井下台阵中，地表记录与井下记录的相似性原理是校正井下地震计方位的有效依据.
2) 利用10个月的连续波形数据，通过计算地表与地下记录的噪声功率谱概率密度函数，定量地分析了该台站的噪声水平，结果表明，对于频率小于1 Hz的微震噪声，两个深度地震计记录到几乎相同的噪声水平，可见增加地震计的埋深，对于压制该频段的噪声水平并无明显的效果，这是由于该频段的噪声波长与地震计深度数量级相当所致(Carter et al, 1991; Boese et al, 2015).同时，该井下地震计在降低频率大于1 Hz的噪声上有明显的效果，尤其对于5 Hz附近的噪声，可降低40 dB之多；当频率大于10 Hz时，随着频率增加，降低噪声的水平逐渐降低.与Boese等(2015)关于新西兰一垂直井下台站噪声水平的研究结果相比，该井下台站的降噪能力表现更优，这与该井下地震计布设的场地条件有一定关系.
3) 利用正则反卷积方法对距离台站150 km内的10个近震事件展开了分析，构建了375 m深度至地表之间的格林函数，并从中识别出入射波及地表反射波，从而获得了浅层的平均速度结构.以其中一个地震事件为例，利用Zhu和Rivera(2002)的频率-波数法计算了直达波和地表反射波的理论到时，与反卷积结果获取的实际到时差基本一致，说明平均速度结构是可靠的.
由钻孔地质资料估算的初始速度模型可以看出，该地区的地下浅层结构中，存在多个较大的波阻抗界面，例如覆盖层及其下面的层位之间有明显的速度差异，当地震波在从井孔基部传播至地表的过程中，发生了多次反射，而传播至地表后，从地表反射回井下地震计的能量明显减小.这是导致反卷积波场上入射波与地表反射波之间还存在多个峰值，而且地表反射波的振幅相对弱于入射波，甚至无法识别的原因.但是，由于只有地表和深度375 m处两个地震计，缺少中间深度的地震计，很难验证这些峰值到底对应于地下哪些波阻抗界面.
这种反卷积方法能够构建不同深度垂直井孔之间的上行波及下行波场，可以用于研究不同层位的地震波传播特征，并获取钻孔区的垂直地震剖面；如果该方法应用于连续长时间的井下观测中，还可以研究井区附近介质特性的变化.因此，本文的研究方法在未来的垂直井下台阵分析中具有广泛的应用前景.
云南省普洱市地震局为本文提供了数据，德国地学研究中心Stefona Parolai教授和Bojana Petrovic博士为本文计算方法和程序提供了很多帮助，作者在此一并表示感谢. -
![]()
图 1 大寨深井台站位置及其记录到的地震分布
黑色圆圈为台站附近的地震分布，黑色实心圆圈为文中反卷积计算用到的地震事件
Figure 1. Location of the Dazhai borehole station (triangle) and distribution of local earthquakes Circles represent the location of local earthquakes and solid circles are the earthquakes used in this paper.
F1 : Wuliangshan fault; F2 : Lancangjiang fault; F3 : Honghe fault
![]()
图 2 大寨井孔结构
两个黑色三角形分别代表位于地表和井下的地震计
Figure 2. Dazhai borehole columnar diagram
Solid triangles indicate sensors at different depths
![]()
图 3 所选取56个地震事件的地表记录与井下记录南北分量(a)和东西分量(b)的傅里叶谱比
Figure 3. Fourier spectral ratio of local earthquakes using NS (a) and EW (b) components between surface and borehole recordings
![]()
图 4 短周期FSS-3DBH地震计的频率特性曲线
Figure 4. Frequency characteristic curves of the short-period seismometer FSS-3DBH
![]()
图 5 ML3.7地震事件南北向分量的原始波形(a)及其2.0 Hz(b)和1.0 Hz(c)低通滤波结果
Figure 5. Original waveforms (a) and low-pass filtered waveforms with the cut-off frequencies of 2.0 Hz (b) and 1.0 Hz (c) of the NS component of the ML3.7 event
Black lines indicate the borehole recordings, and grey lines indicate the surface recordings
![]()
图 6 基于均方根误差方法(a)和互相关函数方法(b)得到的不同地震事件南北分裂的最佳旋转角度统计结果
黑色圆圈表示不同事件对应的角度，虚线表示不同角度的中值及均值
Figure 6. The optimal rotation angles corresponding to NS components of different events by RMS error (a) and cross-correlation (b)
Open circles indicate different angles corresponding to different events, and dashed lines indicate the median and mean value of all rotation angles
![]()
图 7 ML3.7地震事件井下南北分量按照最佳旋转角度旋转前(a)、后(b)的波形对比
Figure 7. Waveforms before (a) and after (b) rotating the NS component of the ML3.7 event according to the optimal alignment angle
Black lines indicate the borehole recordings, and grey lines indicate the surface recordings
![]()
图 8 井下(a)和地表(b)台站三分量的概率密度函数
Figure 8. Probability density function of three components from borehole (a) and surface (b) station
![]()
图 9 地表和井下台站不同分量的概率密度函数中值差值
Figure 9. Difference of the median value between probability density function of three-components recorded by surface and borehole stations
![]()
图 10 2011年11月9日ML3.2地震水平分量反卷积波场
(a)不同正则化参数ε下的反卷积波场, 线条颜色由浅变深，表示ε逐渐减小；
(b)不同迭代次数n下的反卷积波场, 线条颜色由深变浅，表示n逐渐增大Figure 10. Deconvolution wavefield of the event ML3.2 with different regularization parameters
(a) Deconvolution variation with regularization parameter ε decreasing from light to dark lines;
(b) Deconvolution variation with Landweber iterations n increasing from dark to light lines![]()
图 11 表 2中事件2(ML3.2)的原始记录(a)及其相应的反卷积波场(b)
红线表示用于进行反卷积计算的时间窗口，黑色虚线是井下与地表之间相互对称的上行波和下行波
Figure 11. The deconvolution wavefield (b) from the original records (a) of No.2 event listed in Table 2
The red lines indicate the time window used for deconvolution, and the dashed lines show the up-going wave and down-going wave from bottom to surface
![]()
图 12 表 2中地震事件的原始井下记录(a)及其对应的反卷积波场(b)
(a)中红色实线表示用于反卷积计算的时间窗口, (b)中红色虚线分别表示入射波和地表反射波的位置
Figure 12. Aligned deconvolution wavefield (b) from the original borehole records (a) of all the events listed in Table 2
Red lines indicate the time window used for deconvolution, and vertical dashed red lines indicate the arrivals of incident and surface-reflected waves
![]()
图 13 大寨井孔附近的S波(a)和P波(b)速度模型
实心三角形表示两个地震计的深度
Figure 13. Shear velocity model (a) and P wave velocity (b) around the borehole
Triangle indicates sensors at different depths, and dashed lines indicate the estimated S velocity
![]()
表 1 最佳旋转角度的统计结果
Table 1 The statistics of optimal rotation angles
方法 均值 中值 两者平均值 均方根误差法 -6.66° -7° -6.83° 互相关函数法 -5.36° -5° -5.18° 
下载: 导出CSV
表 2 用于反卷积计算的事件列表
Table 2 List of the events used for deconvolution
事件序号 发震时间 ML 震中距/km 年-月-日 时：分：秒 1 2011-11-02 19:33:05.1 3.7 71.21 2 2011-11-09 09:04:23.4 3.2 124.63 3 2011-11-16 23:47:37.5 2.6 121.54 4 2011-12-28 15:16:32.3 3.5 90.64 5 2012-02-22 00:41:39.7 2.3 145.54 6 2012-03-12 01:24:19.1 2.4 68.49 7 2012-04-13 23:06:45.7 2.7 78.48 8 2012-04-29 01:33:01.3 2.7 154.29 9 2012-07-28 08:16:15.5 3.2 57.63 10 2012-07-30 00:05:32.5 4.2(MS) 41.47 注：地震目录引自中国地震台网中心(2015).
下载: 导出CSV
表 3 合成理论地震图所用到的地壳模型
Table 3 Velocity model for calculating the synthetic seismograms
层数 vP/(km·s-1) vS/(km·s-1) 层厚/km QP QS 1 3.75 1.10 0.375 1000 500 2 6.10 3.55 13.85 1000 500 3 6.30 3.65 12.34 1000 500 注：QP，QS表示P波和S波的地下介质品质因子.
下载: 导出CSV
-
冯德益, 张少泉, 卫鹏飞, 王俊国. 1990.深井观测地震波典型记录与分析应用[M].北京:地震出版社: 1-2. Feng D Y, Zhang S Q, Wei P F, Wang J G. 1990. Typical Records and Application of Deep Borehole Seismic Waves[M]. Beijing: Seismological Press: 1-2 (in Chinese).
葛洪魁, 陈海潮, 欧阳飚, 杨微, 张梅, 袁松湧, 王宝善. 2013.流动地震观测背景噪声的台基响应[J].地球物理学报, 56(3): 857-868. doi: 10.6038/cjg20130315 Ge H K, Chen H C, Ouyang B, Yang W, Zhang M, Yuan S Y, Wang B S. 2013. Transportable seismometer response to seismic noise in vault[J]. Chinese Journal of Geophysics, 56(3): 857-868 (in Chinese). doi: 10.6038/cjg20130315
《工程地质手册》编委会. 2007.工程地质手册[M].第4版.北京:中国建筑工业出版社: 1-1099. Editorial Committee of Engineering Geology Manual. 2007. Engineering Geology Manual[M]. 4th ed. Beijing: China Architecture and Building Press: 1-1099 (in Chinese).
李凤杰. 1989.深井地震波观测研究[M].北京:学术期刊出版社: 1-4. Li F J. 1989. Study on Deep Borehole Seismic Wave Observation[M]. Beijing: Academic Journal Press: 1-4 (in Chinese).
李少睿, 毛国良, 王党席, 罗治国. 2016.井下地震计方位角检测技术应用研究[J].地球物理学报, 59(1): 299-310. doi: 10.6038/cjg20160125 Li S R, Mao G L, Wang D X, Luo Z G. 2016. Research on the application of borehole seismometer azimuth detection technology[J]. Chinese Journal of Geophysics, 59(1): 299-310 (in Chinese). doi: 10.6038/cjg20160125
刘渊源, 崇加军, 倪四道. 2011.基于井下摆天然地震数据测量首都圈近地表波速结构[J].地震学报, 33(3): 342-350. http://www.dzxb.org/Magazine/Show?id=27826 Liu Y Y, Chong J J, Ni S D. 2011. Near surface wave velocity structure in Chinese capital region based on borehole seismic records[J]. Acta Seismologica Sinica, 33(3): 342-350 (in Chinese). http://www.dzxb.org/Magazine/Show?id=27826
吕永清, 蔡亚先, 程骏玲. 2007.确定地震计安装方位的相干性分析法[J].大地测量与地球动力学, 27(4): 124-127. http://www.cqvip.com/QK/87801X/200801/26453603.html Lü Y Q, Cai Y X, Cheng J L. 2007. Orientation for seismometer with coherence analyzing method[J]. Journal of Geodesy and Geodynamics, 27(4): 124-127 (in Chinese). http://www.cqvip.com/QK/87801X/200801/26453603.html
仝亚博, 杨振宇, 王恒, 张旭东, 安纯志, 徐颖超, 赵越. 2014.中国西南思茅地体中部白垩纪古地磁结果及陆内地壳变形特征[J].地球物理学报, 57(1): 179-198. doi: 10.6038/cjg20140116 Tong Y B, Yang Z Y, Wang H, Zhang X D, AnC Z, Xu Y C, Zhao Y. 2014. The Creatceous paleomagnetic results from the central part of the Simao terrane in the southwest part of China and its tectonic implication[J]. Chinese Journal of Geophysics, 57(1): 179-198 (in Chinese). doi: 10.6038/cjg20140116
王林瑛, 郭永霞, 刘芳, 蒋长胜. 2008.文安地震前后首都圈分区波速比时变特征[J].地震学报, 30(3): 240-253. http://www.dzxb.org/Magazine/Show?id=26581 Wang L Y, Guo Y X, Liu F, Jiang C S. 2008. Temporal vP/vS variation characteristics in different zones of China's capital area before and after 2006 Wen'an earthquake[J]. Acta Seismologica Sinica, 30(3): 240-253 (in Chinese). http://www.dzxb.org/Magazine/Show?id=26581
谢剑波. 2014.地震记录的时间域反褶积、仿真及在地震计方位角相对测量中的应用[J].地球物理学报, 57(1): 167-178. doi: 10.6038/cjg20140115 Xie J B. 2014. Deconvolution, simulation of seismic records in the time domain and application in the relative measurements of seismometer orientation[J]. Chinese Journal of Geophysics, 57(1): 167-178 (in Chinese). doi: 10.6038/cjg20140115
徐纪人, 赵志新. 2009.深井地球物理观测的最新进展与中国大陆科学钻探长期观测[J].地球物理学进展, 24(4): 1176-1182. doi: 10.3969/j.issn.1004-2903.2009.04.003 Xu J R, Zhao Z X. 2009. Recent advance of borehole geophysical observation and Chinese continental scientific drilling long-term observatory at depth[J]. Progress in Geophysics, 24(4): 1176-1182 (in Chinese). doi: 10.3969/j.issn.1004-2903.2009.04.003
张尉, 陈棋福, 丘学林, 陈颙. 2009.首都圈数字地震台网对微弱爆破信号的检测能力[J].地球物理学报, 52(3): 681-690. http://www.irgrid.ac.cn/handle/1471x/197822 Zhang W, Chen Q F, Qiu X L, Chen Y. 2009. Weak explosion signal detection by the Beijing metropolitan digital seismic network[J]. Chinese Journal of Geophysics, 52(3): 681-690 (in Chinese). http://www.irgrid.ac.cn/handle/1471x/197822
中国地震台网中心. 2015. 历史查询[EB/OL]. [2015-08-07]. http://www.ceic.ac.cn/history. China Earthquake Network Center. 2017.History search[EB/OL]. [2015-08-07]. http://www.ceic.ac.cn/history (in Chinese).
Asai Y, Okubo M, Ishii H, Aoki H, Yamauchi T, Kitagawa Y, Koizumi N. 2005. Co-seismic strain-steps associated with the 2004 off the Kii Peninsula earthquakes-Observed with Ishii-type borehole strainmeters and quartz-tube extensometers[J]. Earth Planets Space, 57: BF03352568. doi: 10.1186/BF03352568.pdf
Aster R C, Shearer P M. 1991. High frequency borehole seismograms recorded in the San Jacinto fault zone, southern California, Part 1: Polarization[J]. Bull Seismol Soc Am, 81(4): 1081-1100. https://pubs.geoscienceworld.org/bssa/article-lookup/81/4/1057
Bertero M, Boccacci P. 1998. Introduction to Inverse Problems in Imaging[M]. Bristol: IOP Publishing: 1-347.
Bindi D, Parolai S, Spallarossa D, Catteneo M. 2000. Site effects by H/V ratio: Comparison of two different procedures[J]. J Earthq Eng, 4(1): 97-113. https://www.cambridge.org/core/books/the-italic-people-of-ancient...
Boese C, Wotherspoon L, Alvarez M, Malin P. 2015. Analysis of anthropogenic and natural noise from multilevel borehole seismometers in an urban environment, Auckland, New Zealand[J]. Bull Seismol Soc Am, 105(1): 285-299. doi: 10.1785/0120130288
Chong J J, Ni S D. 2009. Near surface velocity and QS structure of the Quaternary sediment in Bohai basin, China[J]. Earthquake Science, 22(5): 451-458. doi: 10.1007/s11589-009-0451-1
Carter J A, Barstow N, Pomeroy P W, Chael E P, Leahy P J. 1991. High-frequency seismic noise as a function of depth[J]. Bull Seismol Soc Am, 81(4): 1101-1114. https://pubs.geoscienceworld.org/bssa/article-lookup/81/4/1101
Díaz J, Villaseor A, Morales J, Pazos A, Córdoba D, Pulgar J, García-Lobón J L, Harnafi M, Carbonell R, Gallart J, TopoIberia Seismic Working Group. 2010. Background noise characteristics at the IberArray broadband seismic network[J]. Bull Seismol Soc Am, 100(2): 618-628. doi: 10.1785/0120090085
Fukao Y, Ishibashi K. 1996. Damages Caused by Osaka-Kobe-Awaji Large Earthquake and Earthquake Prediction[M]. Tokyo: Iwanami Press (in Japanese).
Fukushima Y, Kinoshita S, Sato H. 1992. Measurement of Q-1 for S waves in mudstone Chikura, Japan: Comparison of incident and reflected phases in borehole seismograms[J]. Bull Seismol Soc Am, 82(1): 148-163.
Fukushima R, Nakahara H, Nishimura T. 2016. Estimating S-wave attenuation in sediments by deconvolution analysis of KiK-net borehole seismograms[J]. Bull Seismol Soc Am, 106(2): 552-559. doi: 10.1785/0120150059
Grigoli F, Cesca S, Dahm T, Krieger L. 2012. A complex linear least-squares method to derive relative and absolute orientations of seismic sensors[J]. Geophys J Int, 188(3): 1243-1254. doi: 10.1111/gji.2012.188.issue-3
Hauksson E, Teng T L, Henyey T. 1987. Results from a 1500 m deep, three-level downhole seismometer array: Site response, low Q values, and fmax[J]. Bull Seismol Soc Am, 77(6): 1883-1904. https://www.researchgate.net/publication/317931488_Hydraulic...
Ma K F, Lin Y Y, Lee S J, Mori J, Brodsky E E. 2012. Isotropic events observed with a borehole array in the Chelungpu fault zone, Taiwan[J]. Science, 337(6093): 459-463. doi: 10.1126/science.1222119
McNamara D E, Boaz R I. 2005. Seismic Noise Analysis System, Power Spectral Density Probability Density Function: Stand-Alone Software Package[R]. USGS Open-File Report: 2005-1438.
Mehta K, Snieder R, Grazier V. 2007a. Extraction of near-surface properties for a lossy layered medium using the propagator matrix[J]. Geophys J Int, 169(1): 271-280. doi: 10.1111/gji.2007.169.issue-1
Mehta K, Snieder R, Grazier V. 2007b. Downhole receiver function: A case study[J]. Bull Seismol Soc Am, 97(5): 1396-1403. doi: 10.1785/0120060256
Nakata N, Snieder R. 2012. Estimating near-surface shear wave velocities in Japan by applying seismic interferometry to KiK-net data[J]. J Geophys Res, 177: B01308. http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.705.4195
Nakata N. 2013. Combination of Hi-net and KiK-net data for deconvolution interferometry[J]. Bull Seismol Soc Am, 103(6): 3073-3082. doi: 10.1785/0120130101
Nakata N, Snieder R, Behm M. 2014. Body-wave interferometry using regional earthquakes with multidimensional deconvolution after wavefield decomposition at free surface[J]. Geophys J Int, 199(2): 1125-1137. doi: 10.1093/gji/ggu316
Wang X, Chen Q F, Li J, Wei S J. 2016. Seismic sensor misorientation measurement using P-wave particle motion: An application to the NEC saids array[J]. Seismol Res Lett, 87(4): 901-911. doi: 10.1785/0220160005.
Okubo M, Asai Y, Aoki H, Ishii H. 2005. The seismological and geodetical roles of strain seismogram suggested from the 2004 off the Kii peninsula earthquakes[J]. Earth Planets Space, 57(4): 303-308. doi: 10.1186/BF03352567
Parolai S, Ansal A, Kurtulus A, Strollo A, Wang R J, Zschau J. 2009. The Ataky vertical array (Turkey): Insights into seismic wave propagation in the shallow-most crustal layers by waveform deconvolution[J]. Geophys J Int, 178(3): 1649-1662. doi: 10.1111/gji.2009.178.issue-3
Peterson J. 1993. Observations and Modeling of Seismic Background Noise[R]. USGS Open-File Report: 93-322.
Seale S H, Archuleta R J. 1989. Site amplification and attenuation of strong ground motion[J]. Bull Seismol Soc Am, 79(6): 1673-1696.
Snieder R, Safak E. 2006. Extracting the building response using seismic interferometry: Theory and application to the Millikan library in Pasadena, California[J]. Bull Seismol Soc Am, 96(2): 586-598. doi: 10.1785/0120050109
Tikhonov A N, Arsenin V Y. 1977. Solutions of Ill-Posed Problems[M]. WashingtonD C:V.H.Winston & Sons: 1-258.
Zhu L P, Rivera L A. 2002. A note on the dynamic and static displacements from a point source in multilayered media[J]. Geophys J Int, 148(3): 619-627. doi: 10.1046/j.1365-246X.2002.01610.x
Zoback M D, Hickman S H, Ellsworth W L. 2011. Scientific drilling into the San Andreas fault zone: An overview of SAFOD's first five years[J]. Scientific Drilling, 11: 14-28. doi: 10.5194/sd-11-14-2011
-
期刊类型引用(4)
1. 万文涛，陈畅，王赟，穆朝民，贺永胜，汪超. 地面与地下深部地震背景噪声对比分析. 地球物理学报. 2024(02): 793-808 . 
百度学术
2. 徐纪人，李海兵，曾祥芝，许健生，赵志新. 中国井下地震观测研究回顾与展望——从井下到东海深井垂直地震台阵. 地震学报. 2024(06): 919-935 . 本站查看
3. 陈畅，王赟，郭高源，操玉文，李帅，张东明，菅一凡，汪超. 几种旋转地震仪在深部地下巷道的观测对比. 地球物理学报. 2022(12): 4569-4582 . 百度学术
4. 徐纪人，李海兵，曾祥芝，赵志新. 江苏东海深井观测地震波形及其信噪比研究. 地震学报. 2022(06): 1007-1018 . 本站查看
其他类型引用(1)
计量
- 文章访问数: 1443
- HTML全文浏览量: 678
- PDF下载量: 29
- 被引次数: 5
