阵列式光纤地震仪在西藏易贡湖地区的浅层探测实验研究

张衡， 张文涛， 黄稳柱， 裴顺平， 恩和得力海，李云月， 赵俊猛

1 中国科学院青藏高原研究所，青藏高原地球系统与资源环境全国重点实验室， 北京 100101 2 中国科学院大学材料科学与光电子技术学院， 北京 100049 3 中国科学院半导体研究所传感技术国家重点实验室， 北京 100083 4 新加坡国立大学土木与环境工程系， 新加坡 117576 5 普度大学地球大气行星科学系， 西拉法叶， 美国 47907

0 引言

地震层析成像作为一种有效的地下结构探测手段，近几十年来取得了飞速发展.按照震源的不同，地震层析成像主要分为主动源和被动源.主动源以爆破、气枪、电火花等人为产生的震源为主，主要用于工程勘探领域，可以获得高精度的地下浅层结构；被动源则以天然地震和背景噪声(地球内部存在的微弱振动)为主.背景噪声源通常被认为是多种震源经地下介质散射而被地震仪器记录到的一种信号.利用频率低于1 Hz的长周期信号和互相关技术提取经验格林函数的办法，科学家们可以进行大尺度的壳幔结构地震成像研究(Shapiro et al., 2005; Yao et al., 2008).得益于成像方法的进步和计算条件的不断完善、地震数据的积累和采集方式的更新，我们得到了越来越精确的地下结构图像(Liu，2022).但要想进一步提高地下结构的探测精度，尤其是极端环境下(本文中的极端环境特指青藏高原地区内无稳定电力供应以及常规地震仪器难以布设的区域)的地下结构，获取更加丰富的数据仍然是至关重要的条件.

对于野外环境极端的青藏高原地区，虽然科学家们已经利用传统电学地震仪器采集到了数据，并采用地震学方法获得了一系列重大突破(Liang et al., 2016; Owens and Zandt, 1997; Tian et al., 2015; Tilmann et al., 2003; Zhao et al., 1993; Zhang et al., 2021)，但是由于传统地震仪器在极端环境中的应用面临诸多技术瓶颈，导致该区数据采集密度仍然很低.传统地震仪内部有大量电子器件在工作时是需要供电的，而在青藏高原野外不仅供电困难且多冰雪天气，导致仪器稳定性难以保障，且无法满足长期连续工作的需求；另一方面，传统地震仪安装一般要求倾斜角度±2.5°以内，这给在该区大量架设和维护仪器带来了极大的困难，且成本高昂.由于传统水下地震仪器的价格昂贵、不易维护、难以实现数据实时监测和传输，导致整个青藏高原地区虽有数量庞大的湖泊和河流(姜加虎和黄群, 2004; 朱大岗等, 2007)，但水下地震学研究却几近空白.较之非极端环境地区，该区地下结构的探测范围和精度明显不足.因此，我们需要不断尝试使用新的数据采集技术，布控和加密数据采集点，使之填补数据空白区.

近年来，随着光电子技术和光纤传感技术的飞速发展，出现了一批新型的光纤地震仪.根据光学原理和设计方案的不同，主要分为分布式光纤地震仪(亦称为分布式声波传感技术，Distributed Acoustic Sensing, DAS)和阵列式光纤地震仪.分布式光纤地震仪基于光纤的弹光效应，通过探测由弹光效应产生的光信号的光强、该点对应的光信号的偏振状态、相位、波长等参量，可以间接获得该点处的振动频率、温度、应力改变量等参量(Rajan et al., 2016).国际上利用分布式光纤地震仪观测数据开展地震学研究已经取得一定进展，包括主动源和被动源研究(Parker et al., 2018; 宋政宏等, 2020; Yu et al., 2019; Zeng et al., 2017; 张衡等, 2021).国内的分布式光纤地震学研究相对起步较晚(刘德明和孙琪真, 2009; 徐团伟等, 2016)，现有研究主要集中于油气勘探(Yu et al., 2016; 李彦鹏等，2020; 周小慧等, 2021)和城市区域浅层探测(林融冰等, 2020; 宋政宏等, 2020).现阶段的分布式光纤地震仪尚存缺点，比如数据均为单分量，传递函数不清晰，灵敏度较低以及对光纤线缆的布设要求较高等.特别是光纤线缆的布设需要与地面有较高的耦合度，该条件在极端野外环境中(尤其是水下环境)难以得到充分保障.

与分布式光纤地震仪相比，阵列式光纤地震仪具有更低的噪声水平与更宽的测量频带，能够实现三分量测量，并且对安装条件要求不高(Ajo-Franklin et al., 2019)，这有助于在极端野外环境中获得更多的数据，从而提高地下结构成像分辨率和填补研究空白区.目前，国际上各大油气公司引领了阵列式低噪声光纤地震仪研究，法国CGG Sercel公司(Nakstad and Kringlebotn, 2008)、美国Geospace Technologies公司(Ampilov et al., 2019)、英国Stringray Geophysical(Crickmore et al., 2004)、美国Paulsson公司(Paulsson et al., 2013)等均研发了千道级以上的阵列式光纤地震仪器.近年来，美国(Zumberge et al., 2018)、意大利(Gagliardi et al., 2010)、法国(Bernard et al., 2019)、中国(Huang et al., 2018)等还开展了低噪声、宽频带光纤地震仪器研究，进一步拓展了光纤地震仪的测量频带.

中国科学院青藏高原研究所与中国科学院半导体研究所联合研制了针对极端野外环境的高精度阵列式光纤地震仪器.利用该自主研发的仪器，我们在西藏易贡湖地区获得了两条测线的数据，其中一条测线直接布设于易贡湖湖底.本文利用首次获得的西藏地区光纤地震数据，对该区域的近地表构造进行了背景噪声成像，该研究成果为光纤地震仪器在极端野外环境的应用提供了理论和试验依据.

1 阵列式光纤地震仪

图1 基于干涉式外差相位探测技术的阵列式光纤地震传感原理FC：光纤耦合器，CIR：光纤环形器，PD：光电探测器，DAQ：数据采集模块，FPGA：数字门阵列信号处理模块.Fig.1 The principle of array Fiber-Optic seismic sensing is based on interferometric heterodyne phase detection technologyFC: Fiber-Optic Coupler, CIR: Fiber-optic circulator, PD: Photoelectric detector, DAQ: Data acquisition card, FPGA: Field Programmable Gate Array.

图2 阵列式光纤地震仪Fig.2 The picture of array Fiber-Optic seismometers

图3 阵列式光纤探头测量(a) 仪器硬件测量； (b) 仪器频率响应曲线. 单位pm/g表示每一个重力加速度(g)作用下光学地震计产生的光学波长变化为1皮米(pm).Fig.3 The calibration of the fiber seismometer(a) The actual picture of comparative measurement; (b) The measurement result of the frequency response. The unit pm/g indicates that the 1 pico meter (pm) optical wavelength variation generated by the fiber seismometer for each gravity acceleration (g).

2 数据和方法

2020年9月，中国科学院青藏高原研究所与半导体研究所联合在西藏易贡湖地区开展了阵列式光纤地震仪的野外试验，试验场地为易贡湖内和非常湿润的湖岸(图4a).我们在湖底和湖岸分别布设了8个光纤地震探头(以下简称探头，探头编号见图4b)，探头间距5 m(图4b，红色测线1在湖底，绿色测线2在湖岸)，湖中测线距离湖岸大约5 m.每个探头的采样率设置为1000 Hz.光纤解调仪距离1号探头约12 m，供电设备为12 V蓄电池 (图4c).为了能够进一步测试数据质量与仪器布设环境(如探头倾斜角度以及探头与地面的耦合程度)的关系，我们仅在探头上增加了三枚固定长钉，而没有对探头进行掩埋(图4d).置于湖岸上的光纤地震仪探头直接插入泥土中；而对于湖底的光纤地震仪探头的安装，我们只是将其以添加配重的方式沉入湖底.同时，为了便于比较信号质量，我们在公路边架设了一台传统的电学地震仪器(图4b中蓝色方框处).

图4 西藏易贡湖地区野外试验(a) 研究区域； (b) 试验区测线分布；FS：光纤地震仪，五角星代表单次敲击位置，蓝色方框标示出安装在公路边的传统电学地震仪位置； (c) 野外布设光纤信号解调仪； (d) 湖岸探头布设.Fig.4 The field experiment in the Yigong lake, Tibet(a) The study area; (b) The Yigong lake and location of the profiles. FS: Fiber-Optic Seismometer, stars: The locations of two single shots, blue rectangle: The location of the traditional electronic seismometer near the road; (c) The interrogator in the field; (d) The No.12 FS deployed in the riparian.

试验期间，我们使用夯地机在湖岸边的公路进行了敲击(图4b中16号探头附近)，主要用于比较主动震源与被动震源的信号强弱，从而判断光纤地震仪器的信噪比.本次试验总记录时长167 min.对于背景噪声成像而言，增加噪声互相关函数的叠加时间长度会相应地增加信号的信噪比，但是该效果随着叠加时间长度的增加会逐渐趋于稳定(林融冰等, 2020).经过不同叠加时间长度的对比，并且由于整个数据记录过程中同时记录了主动源数据，因此本文选取了90 min包含较少强震动信号的连续数据用于背景噪声成像.图5a—d展示了光纤仪器和传统电学地震仪在水中、湖岸和公路边记录到的主动震源信号(0.5～500 Hz带通滤波)及其相应的频谱特征.图5e和5f则分别展示了光纤仪器和传统地震仪在水中、湖岸和公路边记录到的背景噪声信号(0.5～500 Hz带通滤波)及其相应的频谱特征.低于50 Hz处的背景噪声迅速增加，说明该地区的背景噪声频带主要集中于50 Hz以下.比较背景噪声和主动震源的频谱特征可以发现，主动源信号水平至少比背景噪声高出一个数量级，说明了该仪器具有较高信噪比，因此可以在水下记录到清楚的主动源信号.在距离30 m以外的水下依然能够识别明显的震动信号，甚至可以清晰地识别出夯地机的敲击次数.不仅如此，从频谱分析上可以看出，阵列式光纤仪器的信号和噪声水平与传统电学地震仪器非常接近(噪声水平在50 Hz以下较高是因为光纤仪器布置在湖边和水下，而传统仪器布设在相对易贡湖较远的公路边).

图5 湖岸与湖底的主动源数据、背景噪声数据及其对应频谱(a) 夯地机主动源数据； (b) 夯地机主动源数据频谱； (c) 单次敲击主动源数据，敲击位置见图3； (d) 单次敲击主动源数据频谱； (e) 背景噪声数据； (f) 背景噪声数据频谱.图中绿色线条为湖岸探头记录数据，红色线条为湖底探头记录数据，黑色线条为布设于公路的电学仪器记录数据.Fig.5 Active source data, ambient noise data and their corresponding frequency spectrum in the riparian and bottom of the lake(a) Active source data from the rammer. Green and red colors show the records in the riparian and at the bottom of the lake, respectively; (b) The spectrum of the rammer data; (c) A single shot, the location is shown in Fig.3; (d) The spectrum of the single shot; (e) The ambient noise data; (f) The spectrum of the ambient noise data. The green line is the data recorded by the FS in the riparian, the red line is the data recorded by the FS in the lake, and the black line is the data recorded by a traditional electronic seismometer near the road.

本次研究使用了背景噪声成像方法分别对湖岸和湖底的两条光纤数据进行处理，主要流程：(1)根据频谱图分析(图5f)，主要的噪声频段集中在100 Hz以下，因此将记录到的90 min原始数据降采样至100 Hz；(2)采用“one-bit”法将降采样后的数据进行归一化，并对数据进行0.5～50 Hz的带通滤波处理；(3)利用时频域相位加权叠加算法将长时间的互相关函数进行叠加(Bensen et al., 2007)，图6a和图6c分别展示了以探头9和12为虚拟源，其余探头为接收源的噪声互相关函数；(4)采用广义反射透射系数法(Chen, 1993)进行正演计算，获取相速度频散曲线(图6b和图6d)；(5)采用基于Particle Swarm Optimization的反演算法(Song et al., 2012)获取沿测线的S波速度结构，图7和图8展示了最终的反演结果.分层模型总厚度为60 m，每层厚度为10 m.在反演过程中，由于研究区域内缺乏湖底浅层速度模型，我们使用的初始S波速度结构参照了海底沉积层的速度结构(Essen et al., 1998).同时，由于本研究中的相速度明显低于该文献中得到的海底沉积层相速度，我们进一步降低S波初始速度至120 m·s-1.P波速度通过P波和S波的经验公式计算获得，密度则根据Brocher经验公式进行换算(Brocher, 2005).以第9号光纤地震仪为例，当反演次数达到18次时已经基本达到目标函数最优解(图7a)，且最终频散曲线反演结果与拾取频散曲线的吻合度较高(图7b).

图6 背景噪声互相关函数及其面波相速度频散曲线(a) 湖岸9号探头的噪声互相关函数(黄色透明实线表示平均视速度)； (b) 湖岸9号探头的面波相速度频散曲线；(c) 湖岸12号探头的噪声互相关函数； (d) 湖岸12号探头的面波相速度频散曲线.Fig.6 Noise cross-correlation functions and dispersion curves for phase velocity(a) The noise cross-correlation function of No.9 FS, and the yellow transparent line represents the average apparent velocity; (b) Dispersion curve of the No.9 FS; (c) The noise cross-correlation function of No.12 FS; (d) Dispersion curve of the No.12 FS.

图7 单点反演迭代次数与频散曲线反演结果(a) 目标函数与反演次数曲线图； (b) 拾取频散曲线(绿色虚线)、根据初始模型正演计算得到的频散曲线(黑色虚线)以及反演后得到的最佳频散曲线(红色实线).Fig.7 The iteration time for No.9 FS and the inversion result of the dispersion curve(a) Curve graph of objective function and iteration times of inversion; (b) Picking dispersion curve (green dashed line), the calculated dispersion curve (black dashed line) according to the initial model, and the optimal dispersion curve (red solid line) obtained from inversion.

图8 反演结果与9号探头处地表裂隙(a) 湖岸9号探头处的S波速度反演结果和反演搜索范围(红色实线表示最佳反演模型，灰色实线表示试验模型以及绿色虚线表示反演时采用的速度边界)； (b) 研究区域内沿河岸的二维速度剖面； (c) 湖岸9号探头处的地面布设情况(地表裂隙规模明显小于图3d中12号探头处的地表裂隙).Fig.8 The inversion results and the cracks near No.9 FS in the riparian(a) The S-wave velocity inversion result and the inversion search range at the No.9 FS in the riparian. The red solid line represents the optimal inversion model, the gray solid line represents the trial model, and the green dashed line represents the velocity boundary adopted in the inversion; (b) Two-dimensional velocity profile along the riparian in the study area; (c) The ground layout of the No.9 FS in the riparian. The scale of the surface fissures is significantly smaller than that of the No.12 FS in Fig.1d.

3 结果与讨论

图5a中互相关函数的面波平均视速度约为120 m·s-1，该数值与主动源得到的面波平均视速度(图5a，116 m·s-1)非常接近.由于互相关函数表现出了非常清楚的频散现象，我们得到了效果很好的相速度频散曲线.图6b中的相速度呈现出明显的正频散，可用频段为2～20 Hz.在部分工程类项目研究中发现，当低速层存在于高速层下方时，频散曲线可能会出现反转情况，即相速度随频率的增加而增加，这与我们的频散曲线(频率低于6 Hz时，图6d)和最终反演结果类似(图8b中探头12的速度结构，低速层上方存在高速层).通过正演模拟计算，我们可以在模拟频散曲线中发现类似特征(图9).值得注意的是，湖底测线的主动源面波平均视速度比湖岸测线的面波平均视速度高出约50%(图5a).我们认为产生这种现象的原因是湖底物质(如流沙或者碎石)受到湖水的压力，不会出现规模较大的裂隙，物质平均密度高于湖岸.相反，湖岸上可以直接观测到许多裂隙(图4d)，这些裂隙的存在降低了湖岸浅层物质的平均密度.大量裂隙的存在也可以解释为什么我们得出的近地表速度明显低于常规的浮土层速度(地震波在浮土层的传播速度约为200 m·s-1).

通过互相关函数的平均视速度以及频散曲线的差异(图6)，初步反映了研究区的浅层结构存在较强的横向不均匀性.图8b的最终反演结果更加直观地体现出横向结构的明显差异.在深度小于10 m的S波速普遍低于110 m·s-1，该速度明显低于浮土层速度，说明在这个位置裂隙沿测线普遍存在.该结果与地表观测到的裂隙分布高度吻合(图4d和图8c).随着深度的增加，与9号探头处的速度结构相比，第12号探头处15～35 m深度表现为低S波速异常.结合地表环境来看，我们发现9号探头处的裂隙规模明显弱于12号探头的布设区域(图4d和图8c).这种现象说明这些地表裂隙的成因可能与较深处的结构有关，而不仅仅是由于地表温度变化造成了这些裂隙规模的不同.但是，由于缺乏地表裂隙尺度的量化数据，并且频散曲线在低频处的能量扰动较大(受限于测线长度以及记录时间)，裂隙成因与地下低速异常体的具体联系有待进一步深入研究.

利用单次敲击产生的面波信号(图5c)，我们同样能够得到效果较好的频散曲线(图10a和10b)，并且岸边主动源和噪声源的频散曲线(9号探头)具有高度相似性，关于主动源信号的详细处理和解释将在后期开展具体工作.遗憾的是，利用与湖岸测线相同的噪声数据处理流程和方法，易贡湖底测线互相关函数的频散现象相对较弱(图10c)，只能得到部分探头的频散曲线，且可用频率范围较小(2～4 Hz，图10d)，因此在现阶段我们暂时未对湖底测线的数据做进一步的处理.鉴于布设于湖底的探头能够记录到清晰的主动源信号以及与湖岸测线相似的噪声频谱记录，我们认为无法得到效果良好的互相关函数与仪器本身无关，而是由于：易贡湖是具有较快流动速度的湖泊(测线布设处的平均流量为1500 m3·s-1)，并且考虑到探头在湖底的安装并未进行掩埋，因此流动的湖水对仪器本身造成了一定程度的冲击，这种湖水冲击信号的频率高于背景噪声(从主动源信号的频谱来看，这种湖水冲击信号的频率应该低于主动源，如图5b和5f)，所以对湖底原始数据的处理不能简单地使用湖岸数据的处理流程.在未来研究中，我们将在西藏地区选择有钻孔岩心采样和速度测井的地点进行观测实验，这样一方面可以将反演结果与钻孔观测结果对比，另一方面可以对结果解释加以检验.

图9 正演模拟结果(a) S波速度初始模型；每一层的模型参数显示在图右侧，第一列为P波速度，第二列为S波速度，第三列为每层厚度；(b) 模拟波形； (c) 模拟频散曲线.Fig.9 The results of forward simulation(a) The initial model of S-wave velocity; The model parameters of each layer are shown on the right side of the figure, and the first column is P-wave velocity, the second column is S-wave velocity, and the third column is the thickness of each layer; (b) The simulated waveform; (c) The simulated dispersion curve.

图10 易贡湖地区光纤数据结果(a) 水下主动源面波频散曲线； (b) 岸边主动源面波频散曲线； (c) 水下噪声互相关函数； (d) 水下噪声频散曲线.Fig.10 The results from the FS deployed in the Yigong lake(a) Dispersion curve beneath the water extracted from the active source; (b) Dispersion curve in the riparian extracted from the active source; (c) The cross-correlation function derived from the FS deployed in the lake; (d) Dispersion curve extracted from the ambient noise data.

4 结论

利用本次野外试验获得了该研究区的连续记录数据，对其中90 min的连续背景噪声数据进行处理后，提取出具有明显频散的互相关函数.采用广义反射透射系数法，拾取了频段范围在2～20 Hz间的高质量面波频散曲线.最后，参照海相沉积层的速度结构和面波相速度构建了初始模型，利用基于Particle Swarm Optimization的反演算法得到了沿湖岸测线的近地表S波速度.通过面波相速度和最终反演得到的S波速度，均反映出研究区域存在强烈的横向速度不均匀性.通过对湖底测线的数据进行分析，我们发现虽然湖底的背景噪声数据还需要进一步处理分析，但是湖底探头的主动源信号具有良好的信噪比.本次极端环境野外试验结果初步表明，自主研发的阵列式光纤地震仪器可用于特殊自然环境的浅层结构成像.在后续的研究中，我们将着重开展以下工作：(1)利用阵列式光纤仪器的三分量数据，进行更多地震成像方法的探索；(2)对湖岸和湖底的主动源数据进行细致处理，使用主动源成像方法对地下结构进行探测研究；(3)对湖底的背景噪声数据进行深入分析，找出其与湖岸数据的异同；(4)后续的试验应同时布设海底地震仪与光纤地震仪，以期取得两组数据进行对比研究.

致谢在此对中国科学院青藏高原研究所的刘红兵副研究员、王信国副研究员、中国科学院半导体研究所的张建祥博士以及其他野外工作者提供的大力帮助表示诚挚的感谢.感谢编辑以及评审人对文章提出的建设性意见.