基于应变率的分布式光纤声波传感全波形反演研究

刘辉，李静,2*，迟本鑫

1 吉林大学地球探测科学与技术学院，长春 130026 2 地球信息探测仪器教育部重点实验室 (吉林大学)，长春 130026 3 中国科学院精密测量科学与技术创新研究院, 武汉 430077

0 引言

地震波速度是反映地下构造和岩石物性的重要参数之一，获得精确的浅地表速度结构是地震数据反演成像的核心研究内容.近年来，随着光纤技术的迅速发展，一种地震数据采集新技术——分布式光纤声波传感系统(Distributed Acoustic Sensing, DAS)受到广泛关注.DAS系统可以永久布设于地下，易实现大区域、高密度观测，具有较强抗干扰能力，随需求选择DAS空间采样率，实现超高密度的动态监测，在天然地震观测、城市活动监测以及油气勘探领域具有很大潜力.

DAS系统最早在1990年由Dakin(Dakin and Lamb, 1990)提出，Posey等(2000)研发了基于瑞利散射的光纤应变传感器，使其具有初步应用可行性.之后如Optasense、Silixa等公司相继实现了DAS系统研发.随着仪器系统的发展，Mestayer等(2011)进行了DAS系统的地球物理监测试验，证明其在井中观测的可行性.Mateeva等(2014)与Daley等(2016)分别利用DAS系统实现了VSP储层监测与CO2封存监测.Lindsey等(2017)成功采用DAS技术对天然地震和诱发地震的敏感性进行了测试.之后DAS系统逐步推广到近地表地球物理研究领域(Parker et al., 2014).Dou等(2017)利用DAS记录长时间的交通噪声实现了近地表剪切波成像.近年来，国内也陆续开展了相关研究，张丽娜等(2020)和周小慧等(2021)分别讨论了DAS系统在天然地震学研究和油气勘探领域面临的挑战和发展趋势.Zeng等(2017)对比了DAS观测数据与传统检波器数据提取的频散曲线，证明了二者可以提供相似的地下信息.宋政宏等(2020)和林融冰等(2020)分别利用DAS采集主动源与背景噪声信号拾取频散曲线进行近地表横波速度成像.

速度结构反演是DAS数据处理的重要内容，常规基于层状假设的面波频散反演方法横向分辨率不足.全波形反演(Full Waveform Inversion, FWI)自20世纪80年代提出以来(Lailly, 1983; Tarantola, 1984, 1986)，在时间域(Bleistein et al., 1987; Zhang and Luo, 2013)和频率域(Pratt and Worthington, 1988; Song et al., 1995; Sirgue and Pratt, 2004; Zhang et al., 2021)得到了长足的发展.基于波动方程的全波形反演方法利用完整的波场记录实现对地下介质的精确成像.常规的全波形反演方法使用检波器采集的质点速度，而DAS记录沿光纤的轴向应变率(或应变)，因此无法将全波形反演直接应用于DAS数据.现阶段研究根据应变率与速度关系将DAS信号转换为检波器信号(Daley et al., 2016; Egorov et al., 2018; Zhu et al., 2020)，但在天然地震学领域当远震地震波垂直到达水平DAS阵列时，转换会使低波数数据不稳定(Lindsey et al., 2020)，而在油气勘探领域，会强烈地放大低波数噪声信号(Egorov et al., 2018).

基于上述问题，本文首先利用弹性波有限差分算法对DAS信号进行正演模拟，对比DAS信号与传统检波器信号响应特征及频散特征差异.然后通过修改伴随状态源方程，推导实现了直接基于应变率参数的DAS数据全波形反演方法.模型测试对比了直接反演与转换为检波器信号的全波形反演方法的反演效果，以及不同标距长度对DAS数据反演分辨率的影响.最后将本文提出方法应用于美国实测DAS数据中并得到了稳定的反演结果.

1 方法原理

1.1 DAS基本原理

光与直径远小于光波波长的粒子发生碰撞时会产生弹性散射，散射光的波长和频率不变称为瑞利散射，DAS系统基于光纤中的瑞利散射效应获取光纤周围的地震波场信号.DAS系统包含地面仪器与埋在地下的光纤，地面仪器向光纤中发射特定频率与宽度的光脉冲，光脉冲与光纤中的不均匀散射体发生弹性碰撞产生瑞利散射效应，散射光向四面八方传播.当地震振动使光纤发生局部应变时，会导致光纤中散射体空间位置的改变，进而改变沿光纤反向传播的背向瑞利散射光相位，通过地面仪器解调获得光纤周围的地震波场信息(马国旗等，2020).使用分布式光纤传感器获取的数据与使用地震检波器这类点传感器获取的数据不同.传统检波器测量单个点的质点速度信息，而DAS响应为标距长度内的轴向平均应变率或应变(Daley et al., 2016)，应变信号在时间域求导即可转换为应变率，本文主要以应变率为例进行讨论:

(1)

1.2 DAS信号与传统地震检波器信号转换

1.2.1 检波器信号转换为DAS信号

应变率与速度关系式：

(2)

Daley等(2016)根据应变与速度关系提出了另外一种转换方式：沿x方向传播的平面波，u(x,t)=A(x)ei(kx-wt).

(3)

1.2.2 DAS信号转换为检波器信号

将DAS信号在时间域积分转换为应变ε，根据速度与应变关系(Zhu et al., 2020)：

(4)

(5)

1.3 地震检波器数据弹性波全波形反演

二维各向同性介质一阶速度—应力弹性波方程：

(6)

其中v是速度，σ是应力，ρ是密度，λ和μ是拉梅常数，s为源函数，将方程(6)写成矩阵形式：

(7)

其中A为正演算子，f为源，建立目标函数：

(8)

其中(,)表示内积，m为模型参数，w(m)=(vx,vz,σxx,σzz,σxz)T为模型正演数据，d为检波器观测数据，利用伴随状态法(Mora, 1987; Plessix, 2006; Li et al., 2019)求解目标函数的梯度：

(9)

由于源不随模型参数发生变化，将等式(7)对模型参数m求导：

(10)

移项可得：

(11)

将式(11)代入式(9)，有

(12)

A(m)Tw′=Δd(m)，(13)

其中A(m)T为伴随状态法的反传算子，w′=(v′x,v′z,σ′xx,σ′zz,σ′xz)T，Δd(m)为反传源.完整的伴随状态方程为

(14)

最后根据式(12)、(13)得到目标函数对λ、μ和ρ的梯度：

(15)

根据VP、VS与λ、μ、ρ之间关系，可以获得VP、VS的梯度：

(16)

1.4 基于应变率DAS数据全波形反演推导

根据1.2.1节应变率与速度关系，建立DAS数据的目标函数：

(17)

其中m为模型参数，w(m)=(vx,vz,σxx,σzz,σxz)T，ddas为DAS观测数据，J为空间导数矩阵：

(18)

利用伴随状态法(Mora, 1987; Plessix, 2006)求解目标函数的梯度：

(19)

由于源不随模型参数变化，等式(11)依旧成立，将式(11)代入式(19)：

(20)

其中Δddas为DAS观测数据与模型数据的残差，令{A(m)-1}TJTΔddas(m)=w′，移项，得

A(m)Tw′=JTΔddas(m)，(21)

其中A(m)T为伴随状态法的反传算子，w′=(v′x,v′z,σ′xx,σ′zz,σ′xz)T，JTΔddas(m)为反传源.DAS全波形反演的伴随状态方程为

分别由式(15)、(16)计算得到λ、μ、ρ和VP、VS的梯度.对比传统全波形反演，两者的主要区别在反传的源.由于检波器可以采集多分量，所以传统全波形反演既可以用单分量反演，也可以同时用多个分量反演，但DAS仅接收沿光纤的轴向应变率信号，所以仅反传轴向分量的残差.DAS全波形反演基本流程图如图1所示，基本步骤可以概括为：

图1 DAS全波形反演基本流程图Fig.1 Workflow of DAS full waveform inversion

(1)建立初始模型;

(3)计算理论数据与观测数据的残差，并利用伴随状态法求解梯度 (Li et al., 2017);

(4)利用抛物线拟合法(Köhn, 2011)求最佳步长;

(5)根据梯度和最佳步长更新模型;

(6)重复(2)—(5)步，直到满足迭代停止条件(误差小于阈值或达到最大迭代次数)，获得最佳模型.

1.5 频率域子波估计

(23)

(24)

其中fk是振幅调节因子，其表达式为：

(25)

(26)

师：很好！我们描述一个几何体的时候，如果这个几何体的各个面都是平面图形，我们便可以通过这些平面图形的形状、大小和位置关系描述；如果有的表面不仅仅是平面图形，我们也可以通过这些面的展开图，或者从不同的方向去看，或者这些立体图形的截面来描述．

(27)

2 DAS频散曲线反演模型测试

建立如图2所示的四层夹低速横波速度模型，纵波速度根据固定泊松比给出(VP=1.732×VS)，密度固定为1.8 g·cm3，模型大小为25 m×170 m，低速层深度范围为5～10 m，离散网格步长为1 m.源和接收器均位于地表，震源使用主频30 Hz雷克子波，检波器共170个，DAS道间距为1 m，标距长度g=1 m.图3为层状介质模型测试结果，图3a和3b分别为利用弹性波有限差分模拟传统检波器水平速度分量和DAS信号单炮数据记录，图3c为检波器与DAS第40道波形对比，DAS(实线)代表应变率信号而检波器(虚线)代表质点速度，两者的振幅和相位均存在差异.图3d和3e分别为传统检波器和DAS记录利用拉东变换获得的频散曲线，两者也有很好的相似性，说明本文所述方法可以很好地进行常规检波器信号和DAS信号的转换.利用Dix快速反演方法(Haney and Tsai, 2015)对检波器和DAS数据进行一维频散曲线反演测试，反演结果如图3f所示，横波速度反演结果具有很好的一致性，说明利用DAS数据进行传统的频散曲线反演能获得稳定的横波速度结果.

图2 层状介质横波速度模型Fig.2 S-velocity layer model

图3 层状模型测试结果(a) 传统检波器地震记录； (b) DAS地震记录； (c) 检波器与DAS第40道地震记录对比； (d) 检波器信号频散曲线； (e) DAS信号频散曲线； (f) 一维频散反演结果.Fig.3 Layer model test results(a) Geophone data; (b) DAS data; (c) Comparisons of geophone and DAS data of 40th trace; (d) Dispersion curve of geophone data; (e) Dispersion curve of DAS data; (f) 1D dispersion S-velocity inversion results.

3 全波形反演模型测试

3.1 DAS数据全波形速度反演

图4 (a) 真实横波速度模型； (b) 初始速度模型Fig.4 (a) True S-velocity model; (b) Initial S-velocity model

图5 (a) 检波器地震记录与DAS记录转换后波形对比图； (b) 检波器地震记录t-k能量谱； (c) DAS记录转换后t-k能量谱Fig.5 (a) Comparison of the geophone data (black line) and converted DAS data (red line); (b) Geophone data in the t-k domain; (c) Converted DAS data in the t-k domain

图6a为利用检波器水平速度分量的全波形反演结果，其与真实模型能很好吻合.而将DAS数据转换为检波器数据后的全波形反演结果(图6b)与真实模型误差较大.图6c是采用本文提出的直接应变率全波形反演方法获得的结果，与常规检波器数据全波形反演结果具有很好的一致性.对上述反演结果进行局部抽道对比(图6d)，常规DAS数据转换为检波器数据后的全波形反演结果与真实模型拟合存在一定误差，而DAS全波形反演方法能获得与常规检波器全波形反演相当的效果.图7为三种反演模式下的数据波形对比，图7a和图7b为常规检波器与本文提出的DAS数据直接全波形反演波形对比，观测数据(实线)与反演数据(虚线)拟合较好.采用传统DAS数据转换后的反演结果显示(图7c)，反演数据与观测数据误差较大.图8a和图8b为归一化目标函数曲线与归一化模型误差曲线，转换后的全波形反演(虚线)最终的目标函数值与模型误差均较大，而传统全波形(黑线)与DAS全波形(点线)的目标函数更收敛，模型误差也更小.

图6 (a) 检波器地震数据全波形反演结果； (b) DAS记录转换为检波器记录的全波形反演结果； (c) DAS全波形反演结果； (d) X=60 m处反演结果对比图Fig.6 (a) Geophone data FWI result; (b) Conventional converted DAS data FWI result; (c) Proposed DAS FWI result; (d) Comparison of inversion results at X=60 m

图7 观测数据(实线)与反演数据(虚线)波形拟合对比(a) 常规检波器FWI; (b) 本文提出的直接DAS数据FWI; (c) 转换DAS数据FWI.Fig.7 Comparison of the observed data (solid line) and inverted data (dashed line)(a) Geophone data FWI; (b) Proposed DAS data FWI； (c) Converted DAS data FWI.

3.2 DAS数据全波形反演标距测试

进一步测试标距长度对DAS全波形反演的影响，建立图9a所示含多个低速异常体的横波速度模型，震源使用主频15 Hz雷克子波，DAS道间距为1 m.图9b为初始横波速度模型；图9c、9d、9e分别为标距1 m、8 m、16 m的DAS全波形反演结果；图9f为深度6 m处各标距长度的反演结果对比.可以看出当标距长度较小时，DAS全波形反演能稳定成像浅部各个低速异常体，随着标距长度的增加，受平均效应影响水平分辨率降低.通过上述测试验证了标距长度对DAS数据反演分辨率的影响.

图8 (a) 归一化目标函数曲线； (b) 归一化模型迭代误差曲线Fig.8 (a) Normalized data misfit residual; (b) Normalized model misfit residual

图9 (a) 真实横波速度模型； (b) 初始速度模型； (c) 标距长度为1 m的DAS全波形反演结果； (d) 标距长度为8 m的DAS全波形反演结果； (e) 标距长度为16 m的DAS全波形反演结果； (f) 深度6 m处不同标距长度反演结果对比图Fig.9 (a) True S-velocity model; (b) Initial S-velocity model; (c) Proposed DAS FWI with a gauge length of 1 m; (d) Gauge length of 8 m; (e) Gauge length of 16 m; (f) Comparison of 1D horizontal inversion results of different gauge lengths at a depth of 6 m

4 实测DAS数据全波形反演测试

实测数据来自于美国加利福利亚州的加纳山谷井下阵列站，工区示意图如图10所示，蓝线为DAS光纤分布区域，光纤铺设于加州74号高速公路旁，埋深约15～46 cm，使用Silixa有限公司的Intelligent Distributed Acoustic Sensor(iDAS)系统采集应变率信号，道间距为1 m，震源为扫频振动源，位于五角星处，主频范围为0～10 Hz，记录时间63 s.

图10 DAS工区示意图 (Zeng et al., 2017)Fig.10 Schematic diagram of DAS work area (Zeng et al., 2017)

为进一步验证DAS全波形反演方法的有效性与可靠性，对长线Ι内的184道数据进行反演测试.原始DAS记录如图12a所示，处理流程如图11所示，首先将采样频率降为200 Hz以减少计算量与内存占用，然后去趋势、均值，对数据进行分段，每段长度为1.5 s，接着使用1～20 Hz带通滤波以消除高频噪声干扰，时间域“one-bit”归一化，以及谱白化处理，最后通过互相关计算以及时频域相位加权叠加(time-frequency phase-weighted stack, tf-PWS)(Schimmel and Gallart, 2007)获得虚拟炮集记录(图12b).从虚拟炮集记录中可以看到清晰的面波，通过拉东变换获取频散能量图(图12c)，再根据能量最大值点提取频散曲线(黑色虚线).根据测线空间位置将多个虚拟炮集记录的频散曲线组合成二维相速度剖面(图12d)，可见相速度在高频横向变化较为平滑而低频起伏变化明显.使用Dix快速反演方法(Haney and Tsai, 2015)反演频散曲线，图14a为将多个一维频散曲线反演结果组合成的拟二维剖面，可以看出地下为层状结构，基岩层横向存在起伏变化.

图11 实测DAS数据处理流程图Fig.11 Workflow of DAS data processing

图12 (a) 长线Ι原始DAS记录； (b) 虚拟炮集记录； (c) 频散能量图; (d) 二维相速度剖面Fig.12 (a) Long line Ι DAS record; (b) Virtual shot gather of cross-correlation; (c) Dispersive curve spectrum; (d) 2D phase velocity profile

利用上述虚拟炮集记录进行DAS全波形反演测试，炮间距为8 m，共23个炮集记录，初始模型为速度随深度递增的层状模型，图13为利用1.5节震源校正滤波估计的第10炮震源子波，试探子波为8Hz雷克子波.采用频率多尺度反演策略(张文生等，2015)，选取0～10 Hz,0～15 Hz和0～20 Hz三个尺度依次进行反演，上一尺度的反演结果作为下一尺度的初始模型.最终反演结果如图14b所示，对比频散曲线反演结果(图14a)可见基岩层的起伏具有很好的一致性.放大的波形拟合图(图15)显示，观测(实线)与反演(虚线)的波形拟合较好.进一步提取一维频散反演(图14a)与DAS-FWI反演(图14b)结果频散曲线，并对比观测记录频散曲线(图16).DAS-FWI结果(图16b)相对于频散反演结果(图16a)与观测记录频散曲线(黑色虚线)误差更小，验证了DAS全波形反演方法获得的横波速度结果相比于一维频散反演速度有更高的反演精度.

图13 子波估计结果Fig.13 Estimated source wavelet

图14 (a) 一维频散曲线反演结果组合成的拟二维横波速度剖面； (b) DAS全波形横波速度反演结果Fig.14 (a) 2D S-velocity profile composed of the 1D dispersion curve inversion; (b) DAS FWI S-velocity tomogram

图15 第10炮DAS全波形反演波形拟合图Fig.15 The DAS FWI waveform fitting of shot 10th

图16 观测记录频散能量图(a) 频散反演结果频散曲线； (b) DAS-FWI结果频散曲线.Fig.16 Dispersive spectrum of observed data(a) Dispersion curve of 1D inversion data; (b) Dispersion curve of DAS-FWI data.

设置如图17a所示的棋盘模型开展DAS反演方法分辨率测试.每个棋盘小方格的大小为6 m×12 m，根据实测数据反演横波速度平均值约为300 m·s-1，速度扰动为±5%，震源为8 Hz雷克子波，其他参数与实测数据相同.DAS全波形反演结果显示，大多数棋盘目标能被准确地恢复(图17b).对于左右两侧边缘异常恢复效果较差，因此对于实测数据反演结果的这部分区域有待进一步确认.

图17 (a) 横波速度棋盘模型； (b) DAS全波形反演结果Fig.17 (a) S-velocity Checkerboard test model; (b) DAS-FWI S-velocity tomogram

该区域近地表层由厚18～25 m的湖床冲积物组成，存在的土壤类型有粉土质砂、砂、粘土质砂和粉土质砾石等(Youd et al., 2004).因此我们对反演结果的综合解释是第1层由粉土质黏土、和黏土组成，深度为0～6 m，横波速度为140～220 m·s-1；第2层由砂质粉土和粘土质砂组成，深度为6～12 m，横波速度为220～300 m·s-1; 第3层为粉土质砂和砂组成，深度为12～18 m，横波速度为300～380 m·s-1；最后一层存在横向起伏变化，由砾质砂和粉土质砾石组成，深度为18～25 m，横波速度为380～440 m·s-1.通过上述测试表明DAS全波形反演方法可以适用于实测DAS数据，并获得稳定的反演结果.

图18 (a) 两层速度模型； (b) 检波器记录； (c) DAS记录Fig.18 (a) Two-layer S-velocity model; (b) Geophone data; (c) DAS data

5 讨论

DAS理论上仅对沿光纤的轴向应变(应变率)敏感，当光纤埋于地表时，其对沿光纤传播的直达波与面波较为敏感，但反射波被接收时与光纤存在一定夹角，会对观测到的记录产生影响.为此建立如图18a所示的两层速度模型进行测试，VP=1.732×VS，模型区域为100 m×200 m，离散网格步长为1 m，震源使用主频30 Hz雷克子波，位于地表X=100 m处，道间距为1 m.图18b和18c分别为利用弹性波有限差分正演获得的常规地震检波器记录速度分量与DAS记录应变率分量，从检波器和DAS的观测记录中均能看到直达P波，瑞利面波，反射PP、PS、SS波，直达S波混夹在面波中.但DAS记录中反射波能量较弱，特别是PP波，在近偏移距时基本消失.说明在反射观测系统下DAS与检波器相比对反射波振幅响应较弱.

实测信号中往往存在大量噪声，振幅较弱的反射波信号容易被淹没，采用螺旋式缠绕光纤等方式提高DAS系统的灵敏度以及信噪比将会是未来一个重要的研究方向.

6 结论

本文通过弹性波有限差分方法对DAS地震信号进行数值模拟，提出了直接基于应变率的DAS数据全波形反演方法，并对理论数据和实测数据开展了一维频散曲线反演以及二维全波形反演测试.结果表明：DAS信号与检波器信号频散曲线总体趋势一致，仅低频存在些许差异；将DAS信号转换为检波器信号后，会出现一些低波数波形干扰，而利用DAS全波形反演方法无需转换，能获得更稳定的反演效果；随着标距长度的增加，受平均效应影响，DAS全波形反演分辨率降低.

在实际采集的野外数据中，往往存在大量的噪声，长标距有助于压制噪声.高质量的数据是反演的前提，建议在保证数据质量的前提下，选择尽可能小的标距长度，以提高速度成像分辨率.

致谢感谢匿名审稿人提出的宝贵修改意见，感谢美国威斯康星大学麦迪逊分校提供的实测DAS数据.