CN104459784B - 基于单台、双台和双事件数据二维Lg波Q值层析成像方法 - Google Patents

Abstract

本发明公开基于单台、双台和双事件数据二维Lg波Q值层析成像方法，包括如下步骤：(a)区域地震资料收集及预处理，获得Lg波真振幅和信噪比，设定信噪比的门限值，选择用于Lg波Q层析成像的Lg谱数据体；(b)从单台数据中提取双台和双事件数据；(c)基于振幅数据的Lg波Q值层析成像。中国东北地区与朝鲜毗邻，对该区域的层析成像结果可用于数确测定朝鲜核爆的体波震级，准确估计核爆当量。另外，中国东北地区富含丰富的油气、矿产、热液资源，这些都与地壳Q值有一定的对应关系；地壳与上地幔的物质交换对地壳衰减的影响强烈，反过来，利用Q值成像模型可进一步研究中国东部岩石圈减薄。

Description

基于单台、双台和双事件数据二维Lg波Q值层析成像方法

技术领域

本发明涉及一种地震Lg波Q值测量技术，具体是一种基于单台、双台和双事件数据二维Lg波Q值层析成像方法。可应用于地下核爆地震学监测。

背景技术

近些年来，在我国周边地区的一些国家，例如印度、巴基斯坦、朝鲜相继进行了地下核试验。这些试验区至我国地震台网的距离在区域地震范围内。虽然这些国家核试验的当量较小，仍可激发可明显观测到的Lg波。利用地震记录快速地识别地下核爆事件，并准确地测定其爆炸当量，对于我国国土安全具有重要意义。准确估计核爆当量需要获得高精度的事件震级。再利用Lg波振幅测量震级时，需要知道台站到事件路径上Lg波Q值的大小。现有欧亚大陆地区Lg波Q值分布图像的分辨率很低，尤其是台站分布较为稀疏的区域，无法用于高精度的Lg波体波震级测定。

发明内容

针对现有技术中存在的不足，本发明的目的在于提供一种基于单台、双台和双事件数据二维Lg波Q值层析成像方法，使得能够利用中国数字地震台网宽频带地震仪的记录，测定中国东北及邻近地区的Lg波Q值，可以为该地区区域震级的测定提供基础数据。

本发明的技术方案是这样实现的：基于单台、双台和双事件数据二维Lg波Q值层析成像方法，包括如下步骤：(a)地震资料预处理，获得Lg波真振幅和信噪比，设定信噪比的门限值，选择用于LgQ层析成像的Lg谱数据体；(b)从单台数据中提取双台和双事件数据；(c)基于振幅数据的Lg波Q值层析成像。

本发明的有益效果是：本发明的基于单台、双台、双事件的二维Lg波Q值层析成像方法，使得能够利用中国数字地震台网宽频带地震仪的记录，测定中国东北及邻近地区的Lg波Q值，可以为该地区区域震级的测定提供基础数据。另外，中国东北地区富含丰富的油气、矿产、热液资源，这些都与地壳Q值有一定的对应关系；地壳与上地幔的物质交换对地壳衰减的影响强烈，反过来，利用Q值成像模型可进一步研究中国东部岩石圈减薄。

附图说明

图1-1 LgQ成像区域为深色的地形部分，扩展的地图显示了CNDSN(方块)和GSN(三角)的地震台站的位置。图中的十字丝表示所用事件的震中位置。图中还标出了大兴安岭(DMs)，长白山(CMs)和渤海湾盆地(BBB)的所在区域。

图1-2a射线数与震中距离的关系；事件总数为125，射线总数为1720。

图1-2b射线数与台站的关系。

图1-3显示Lg波传播的示例地震图。这些地震记录来自事件20040324的地表速度记录。归一化的垂直分量速度记录以台站到震中距离为序排列，台站代码，震中距和最大振幅标注在记录左端，图中短棒标出了典型的视群速度值。内附地图显示了Lg波传播路径。

图1-4模拟的20040324事件的WWSSN短周期垂直分量地表位移记录图。归一化的垂直分量位移记录以台站到震中距离为序排列，台站代码，震中距和最大振幅标注在记录左端，图中短棒标出了典型的视群速度值。

图1-5a-g Lg波振幅谱确定示例：(a)实际观测的地震记录图；(b)除去仪器响应的地表速度记录；阴影部分标出了截取的Lg波列及等时间长度的P波前噪声；(c)Lg波地表位移；(d)P波前噪声；(e)Lg波列与噪声时间序列的振幅谱；(f)参考频率点上的信噪比；(g)参考频率点上的真振幅值。

图1-6a-e LgQ层析成像的分辨率分析示例。(a)由中国东北及邻近区域的1HzLg谱所确定的事件-台站几何，图中绘出了射线，标出台站(白色三角)及事件(十字丝)；(b)基于1°×1°网格剖分的正、负相间Q值扰动叠加的初始Q值模型；(c)输入的任意的源函数；(d)重建的源函数；(e)重建的Q值扰动量。

图1-7a-d不同网格尺度的Q值重建分辨率。从a至d，网格单元分别为0.7°×0.7°,1°×1°,1.5°×1.5°,2°×2°。

图1-8a-b(a)在0.05至10Hz之间的参考频率上的射线数；(b)迭代误差曲线。

图1-9中国东北及邻近地区大地构造分布图(改自Xie,2000and Ge and Ma,2007)。研究区域包含3个大型板块：西伯利亚克拉通(I)，中国东北拼贴板块(II)，和华北克拉通(III)。拼贴板块由褶皱带(II1～7)和抽板块构成(II8～14)。

图1-10a-f典型的Lg波Q值分布及其分辨率。图中标出了参考频率。Q值分布图像采用统一的色标，列于最下端。

图1-11a-b不同地质单元的宽频带范围Q值统计结果。灰色十字丝表示成像反演的Q值，黑色实心圈表示在参考频率点上的统计Q值，误差棒显示了统计误差。统计Q值能够描绘Q值随频率变化。Q0值标示于每个框图中。

图1-12不同地质单元的统计Q值随频率的变化。

图1-13-(1-11)每个框图中为同一事件不同频率的Lg波源函数及其拟合曲线。稠密的浅灰色曲线指的是由Bootstrap方法得到的拟合源模型，中间的黑色曲线为最优源模型。在每个框图内标出了标量地震矩M₀和拐角频率f_c。框图对应的事件标于上方，与表3-2中的事件一一对应。

图1-14a-b成像结果与合成的源模型的比较。

图1-15 M₀和f_c的关系。实线为假定应力降为常数时的线性拟合，虚线是采用不同的应力降及-3的指数率所得的M₀和f_c的线性关系。

图1-16a-c标量地震矩与全球震级的关系。

图1-17显示了单台、双台、双事件的位置关系。

具体实施方式

一、区域地震资料

从1995年10月到2007年8月间，在中国东北及邻近地区发生的125个区域事件，在20个台站记录到的1720个宽频带垂直分量数字地震图。这些资料用于Lg波Q值成像，其相应的台站代码和位置列于表3-1。9个台站，MDJ,INCN,BJT,HIA,YAK,TLY,ULN,YSS和MAJO，归属于全球地震台网(GSN)，由美国国家地质调查局(USGS)和国际地震学联合会(IRIS)管辖，其它11个台站为中国国家数字地震台网(CNDSN)的台站，由中国地震台网中心(CENC)管辖。GSN和CNDSN的台站都配置宽频带仪器，在0.03至8.0Hz的带宽内具有平坦的速度仪器响应曲线，采样间隔为20,40或50s-1。表3-2列出了事件的参数。研究中用到的台站和事件示于图1-1，其中三角和方块表示台站，十字丝表示事件，图中深色的地形图为探测区域。

表3-1台站代码与位置

许多关于地壳结构的调查表明，研究区域内Moho面埋深存在较大的变化，例如，在渤海湾盆地最浅，约29km，在大兴安岭北部地区深度接近50km。因此，我们仅选择深度小于50km的地壳事件。选择事件的震级范围在3.0到6.0之间，以避免较大震级的事件破裂复杂性的影响。最小的震中距离为192km，因而在观测记录上Lg波波列十分清晰。通常，需要对每一条垂直分量地震图进行可视化的处理，删除饱和波形，噪声干扰较大的波形，以及具有错误到时的波形。这些问题通常是由过大的震中距离，震级太小，多事件叠加等引起的。尽管数据选择的处理过程会导致大量的数据移除，仍然可以获得可用于Q值成像的数据系列。图1-2显示了事件对于震中距离的分布和每个事件的记录数。结合图1-2、图1-3分别可以看出此区域事件对于震中距离和每个事件的记录数和Lg波波列具有很明显的视速度值，3.5km/s。

表3-2区域事件列表

续表3-2区域事件列表

*Beijing regional network 北京区域地震台网

Sakhalin regi onal network Sakhalin区域地震台网

China earthquake network center 中国地震台网中心

§Baykal regional network Baykal区域地震台网

作为示例，图1-3显示了由发生于2004年3月24日的区域地震事件激发的垂直分量地震图。波形为归一化的垂直分量地动速度曲线，在每个地震道的左端列出了最大振幅值，其单位为μm/s。在每道波形记录上标出了视速度值，单位是km/s。在图中可以看出，Lg波波列具有很明显的视速度值，3.5km/s。在内附的地形图中画出了数据与事件的大圆路径。

二、地震资料预处理

地震Lg波的Fourier振幅谱是进行地壳Lg波Q值成像的基础数据。收集每一个事件的地震记录并计算其Fourier振幅谱是进行Lg波Q值层析成像的准备工作。数据预处理分为4个步骤：(1)分离Lg波列与P波前噪声；(2)计算Lg波和噪声的Fourier振幅谱；(3)拾取谱振幅值；(4)去噪处理。

为了识别Lg波列并选择有效的Lg波时窗，使用一种递归的零相位四阶Butterworth滤波器对地震记录进行滤波处理，其拐角频率设定为1Hz，从而得到增强的Lg震相以便于确定其群速度窗口。或者在宽频带地震图上卷积一个全球范围标准地震台网(WWSSN)的短周期仪器响应，并通过地震图可视化来选择Lg波群速度窗。因为WWSSN仪器响应的高通属性，可识别的离散的Lg波在WWSSN地震图上清晰可见。图1-4例举了WWSSN Lg波形，这些波形是在图1-3给出的地震图上卷积WWSSN仪器响应得到的。通过对全部记录的可视化检测，可以确定Lg波列的群速度窗口为3.6-3.0或3.5-2.9km/s。噪声序列一般在P波初至前与Lg波列等时窗内提取。

对于Lg波列和噪声，无论是群速度窗还是时间窗，一旦确定，可以在宽频带地震图上轻松获得有效波列的持续时间。一般在波列两端各增加时窗长度的10％用于余弦镶边消减。对镶边的信号进行快速Fourier变换(FFT)，得到Lg波和噪声的振幅谱。图1-5显示了一个Lg波Q值成像数据预处理的实例。例中显示了Lg波波形，20％的余弦镶边的用法，Lg波及噪声的振幅谱，以及每个参考频率点上的信噪比。由图1-5e可以看出，Lg波振幅谱在1.0Hz附近的频带内具有足够的信噪比，信噪比与震源尺寸，仪器响应的截止频率和震中距离有关。

Lg波有效频带的低频端(0.01-0.1Hz)Lg谱对地震矩(M₀)起主要的约束作用。在0.05至10.0Hz的频带范围内选取参考频率点，并使之在对数域中为等间隔，即，Δlog₁₀f_ref＝0.04。对于每个参考频率(f_ref)，与其对应的Lg波振幅谱(A(f_ref))采用均方根的方法来求取，即，其中f_i∈[f_ref-Δf,f_ref+Δf]，这里的Δf可由Δlog₁₀(Δf)＝0.02来确定。这样，在0.05至10.0Hz的频段内选择58个参考频率点。图1-5g显示了对于参考频率点的Lg谱取值。这些数值较好地反映了Lg波振幅谱的变化。

一些台站的记录具有较大的背景噪声，特别是震中距较大的台站记录，如图1-3所示。进行噪声校正可以有效提高谱振幅测量的精度。假定Lg时窗内的波列为Lg波与噪声的能量叠加，并认为Lg信号与背景噪声的谱振幅是不相关的，因此可以利用下式估计每个频率点上的Lg信号，

其中A表示谱振幅和窗口长度，下标sig、obs和noi分别指真振幅、观测值，以及噪声。

由图1-5g可见，0.2Hz附近有3个参考频率点上没有Lg波真振幅取值，是因为相应的频率点上信噪比的取值小于2.0，如图1-5f。

经过以上4个步骤，可以获得Lg波真振幅和信噪比。一旦设定信噪比的门限值，很容易选择用于LgQ层析成像的Lg谱数据体，如图1-5f。

从单台数据中提取双台和双事件数据：

(1)单台数据

给定一个频率f，对于事件k在台站i观测到的Lg波频谱振幅A可以表示为：

A_ki＝S_kG_kiΓ_kiP_iR_ki (3-2b)

S_k是源项，G_ki＝(Δ₀Δ_ki)^-1/2，Δ_ki≥100km是几何扩展因子,Δ_ki是在事件k和台站i的震中距离，并且Δ₀是固定在100km的基准距离。衰减项可以被表示为：

V是Lg波的群速度，是事件k到台站i沿大路径的积分，Q(x,y,f)是Lg波质量因子，是频率的函数，陆地位置(x,y)。P_i是台站响应，由于当地地球结构，用一个乘法系数控制台站放大率。R_ki是在Lg波在地震k和台站i之间传播的随机效应，通常忽略计算。

(2)双台数据

两站记录同样的地震，他们在或大约在一条线上，双台振幅可通过计算标量谱比。理想情况下，两个记录站i和j与源k对齐，双台振幅可以准确的通过下式计算：

A_kj是记录事件k在j台站的振幅谱。实际情况下，两台站也许没有对齐完美,并存在一定的公差。使用两台站的方位角差，设置最大15°的公差。并设置一个参考点l在射线路径上kj让距离kl等于距离ki，并且需要在点i和l之间的距离在基于频谱分辨率上的半个网格空间。当台站符合这一标准，我们大致计算出它们站间振幅：

基于假设，有相似的传播路径和台站，所以在台站i和参考点l观测的振幅是相同的。

(3)双事件数据

类似双台振幅，当台站记录两个地震，并且他们在或大概在一条线上，双事件振幅可通过计算标量频谱比。理想条件下，台站j与源k和h对齐，双事件振幅可以表示为：

A_hj是台站j对于事件h记录的振幅频谱。实际情况下，两个事件并不完全对齐。我们设置一个参考点q在射线路径kj上，让距离qj等于距离hj。当事件h和参考点q之间的距离不超过半格空间，双事件振幅可以通过如下近似计算：

基于假设，考虑到相同的源项和相似的传播路径，事件h和参考点q在台站j观察到的振幅是相同的。

三、Lg波Q值层析成像

将用于速度结构成像的标准的走时层析技术扩展到LgQ值结构重建，具体方案包括：

(1)假定Lg波在地壳内沿震源至台站的大圆路径传播，将LgQ看成位置与频率的函数，Q(x,y,f)，其中x和y分别表示离散区域的经度和纬度，f为参考频率；

(2)根据先验地质信息给定初始Q值模型，任意的参考频率点上的震源函数的初始值设定为1.0；

(3)基于给定的Q模型和震源函数，合成Lg波振幅谱。将合成的Lg波振幅谱与观测值进行比较并计算残差，以此建立方程组求取Q值与源函数的修正量；

(4)对Q值模型和震源函数进行修正；

(5)进入下一步线性化并返回步骤(3)，直到满足某一收敛条件后终止迭代。

下面给出Lg波Q值层析成像的详细过程。

地震Lg波振幅谱可以表示为

A(f,Δ)＝S(f)R(θ)G(Δ)Γ(f,Δ)X(f)r(f), (3-2)

其中f为频率，Δ是震中距，A(f,Δ)为Lg波位移谱，源函数S(f)可表示为

其中S′(f)是归一化的震源谱，M₀是地震矩，ρ和v_s是源区的平均密度值和剪切波的平均速度。

G(Δ)＝(Δ₀Δ)^-1/2 (3-4)

是几何扩展因子，Δ₀是参考距离，通常设定为100km，R(θ)是源辐射花样，随方位角θ变化而变化。X(f)和r(f)分别为与频率相关的台基响应函数和波传播过程中的随机效应。通常情况下，宏观世界们的影响很难从震源和衰减影响中区分开。折衷地考虑，通过假定X(f)＝r(f)＝1，忽略这些因素的影响，而这些影响可以认为是震源或路径衰减测量中的误差。Lg波通常对方向的影响很弱，在实际计算中令R(θ)＝1。Γ(f,Δ)为衰减的函数，可以表示为

其中V是Lg波的群速度，B(f,Δ)是Lg波Q值沿大圆路径的积分

其中Q(x,y,f)为Lg波Q值，它是频率与地表位置(x,y)的函数。在衰减成像过程中，积分(3-6)式可以离散为

其中k指第k个网格单元内的射线路径，为该单元内的积分。

对方程(3-2)式两端取自然对数，可得

令

ln[S(f)]＝ln[S⁰(f)]+δln[S(f)], (3-10)

并代入(3-7)和(3-8)式可得

其中上标0表示初始模型或反演过程中某步迭代后所得的源或Q值模型，δB(f)可以表示为

矩形网格化简单易行，是普遍采用的方法之一。在每一个矩形单元内有四个角点上的Q值与通过它的射线相关。单元内任意点上的Q值可以用双曲线性函数表示

Q(x,y)＝a₀+a₁x+a₂y+a₃xy, (3-13)

其中双曲线性系数a_p(p＝0,1,2,3)把单元内Q值与角点Q值联系在一起，可以通过下式计算

其中L＝4，是矩形网格单元的节点数。Q_l是第l个节点上的Q值，Δ_lp是按Cramer法则解4阶线性方程组的4阶行列式Δ的代数余子式。按微分的定义，可得

把方程(3-13)代入(3-15)式，可得

其中

把(3-16)式代入(3-12)式，得

对于多台多事件，由方程(3-11)和(3-18)能够建立线性系统

H＝A·δQ+E·δU, (3-19)

其中

是观测的Lg波振幅谱与合成数据的残差，以及

是微分系数，i，j，k分别指统一编号网格单元节点的序号，射线序号，源函数的序号。Q⁰是Q⁰(x,y,f)的简化表达。δq_i指第i个节点上的Q值扰动量。E是一个(0,1)矩阵，当第j条射线与第k个源函数有关时，矩阵E的元素e_kj＝1。δu_k＝δ[ln(S_k(f)]为源的扰动量。用LSQR算法求解线性方程(3-19)式，一般迭代4～6次即可收敛。

LgQ层析成像的迭代收敛条件通常采用如下形式：

其中E(k)指第k次迭代后的系统误差，M是射线方程总数，A_j(f)_obs和分别指第j条射线方程中观测的Lg谱和第k次迭代后的合成谱。

Lg波源激励模型

本发明介绍的Lg波Q值层析成像是区域Q值与震源函数联合反演的方法。因此，有必要讨论一下Lg波源激励模型。作为成像结果之一，不同频率的震源函数可反映事件源的特性。这一估计值是否准确，也是评价方法可信性的重要因素。

Lg波源模型，通常对天然地震事件采用ω²模型，对爆炸事件采用简化的Mueller-Murphy模型

其中M₀和f_c是源的标量地震矩和拐角频率，ρ和v_s分别为地壳的平均密度值和平均的剪切波速度，一般在中国东北地区取值分别为2.7g/cm3and3.5km/s。β表示爆炸源的过冲效应，与源区介质有关系，通常对Poisson介质取值为0.75。对于某一事件，得到一系列离散频率点上的源函数后，可以利用上述源模型反推源参数(M₀和f_c)。地震矩M₀，不仅与拐角频率f_c有关系，而且与应力降Δσ有关系

许多关于震源的研究推断，在中亚和西美地区，源参数M₀与f_c的关系为-4次方关系，而不是-3次方的关系。M₀与f_c的-3次方的关系仍然是正确的，而偏差的存在是因为非常数应力降引起的。利用Kanamori的公式通过标量地震矩可以计算矩震级，

M_w＝[log₁₀(M₀)-9.1]/1.5, (3-26)

其中M₀的单位是Nm。

成像分辨率分析

对于小规模的反演计算，可以采用分辨率矩阵，信息密度矩阵，协方差矩阵来评价利用数据进行模型估计的分辨率。然而，对于大规模的反演计算来说，因为它们的母矩阵，广义逆矩阵本身难以获得，所以不能利用这些信息进行分辨率估计。在本研究中，我们采用“检测板”技术来检验Lg波Q值成像的分辨率。检测板技术通常在旅行时层析成像、波形成像及区域震相传播的效果成像中。至今未见在Lg波Q值成像中采用检测板技术的报道。首先，给定一个平均的Q值模型，在此基础上叠加行列正、负交替的Q值异常值来构建检测板。Q值扰动量是实际Q值的一个百分比。基于检测板的Q值模型，在给定任意震源函数后可以利用方程(3-2)合成Lg波振幅谱。震源与台站的几何关系与实际资料一致。采用与处理实际资料一致的方法，手段，程序和区域网格化，利用合成的Lg波振幅谱进行反演计算。数据对模型空间内任意一个节点上Q值的分辨力可以由重建的Q值异常与给定的Q值扰动量的比值来确定。震源函数的重建能力也是评价层析成像方法的重要因素，可以通过比较重建的和给定的源函数来确定。

针对中国东北及邻近地区的研究区域，在0.05至10.0Hz之间的任一频率平面上，震源与台站的几何关系，Lg波真振幅和相应信噪比都是已知的。通过设定信噪比的阀值(2.0)来构建用于成像的Lg波谱振幅数据体。对于某一事件，如果只有一个Lg波谱振幅数据可用，那么该事件将被删除以保证源函数估计的精度。这样，在某一参考频率上，具有有效Lg波谱振幅数据通过其事件与台站几何关系可以确定成像的观测系统。然而，不同频率上的可用数据是不同的，所以射线覆盖不同，区域分辨率也不同。下面，我们解释如何获得地壳Q值和源函数重建的分辨率，以1Hz Lg波谱振幅为例。

图1-6a显示了由可用数据确定的1370条射线路径，其中十字丝表示事件，三角表示台站。很明显，中国东北地区的射线覆盖比其周边地区密度大。一般来说，射线覆盖越大，分辨率越高。图1-6b显示了检测板模型，网格剖分为1°×1°，Q值扰动量为±30，背景Q值为420。图1-6c绘出了125个1Hz源函数，该函数可由全球震级估算。通常可以假定全球震级与矩震级相当，因此可以通过方程(3-26)得到地震矩；然后固定应力降Δσ为1.0，可以利用方程(3-25)得到拐角频率；将M₀和f_c代入方程(3-23)或(3-24)，即可得到任意频率的源函数。这样，我们可以基于图1-6a给定的观测系统进行正演计算，合成Lg谱振幅。接下来，将合成的Lg波谱振幅输入到成像方程(3-19)式，反演区域Q值分布和震源函数。反演结果的源函数如图1-6d所示，与给定的源函数相比较，可以发现相关系数为0.99，标准差为0.03。这意味着我们的成像方法具有较强的源重建能力。反演结果的Q值示于图1-6e，可以看出，中间区域的分辨率较周边区域高，与射线覆盖密度是一致的。

将输入的检测板模型与重建的模型进行比较，可以评价成像反演的分辨能力。可以定义一个节点分辨率值

其中ΔQ_i和分别为真实的和重建的Q值扰动量，M是测试次数。当R值小于0时，固定R＝0。以测区内的分辨率测量R，作为地壳Q值模型内某一点上扰动量重建的的分辨能力评价值。图1-7显示了不同网格剖分尺度的Q值重建分辨率。从a至d，网格剖分分别为0.7°×0.7°,1°×1°,1.5°×1.5°,2°×2°。一般地，高密度的剖分网格有利于准确拾取路径Q值，从而更加容易地应用到衰减结构分析及震级测量中。然而，在固定的观测系统中，剖分网格密度太大则使Q值重建不准确。例如，图1-7a显示了非均匀的重建结果。随网格密度减小，Q值重建效果更佳。换言之，因为在Q值成像效果与网格密度之间存在较强的依赖关系，寻找最优网格单元至关重要。由图1-7b可以看出，1Hz Lg波谱数据在中国东北的中间区域具有1°的分辨率。

总之，在某一频率平面上对层析成像进行分辨率分析，主要包括以下几步：(a)采用均匀Q值模型为初始模型，一般采用Q值的频率依赖模型Q(f)＝Q₀(f)·f^η计算Q值，其中Q₀固定为420，η设定为0.15；(b)在给定初始模型的基础上施以±7％的Q值扰动量，构建目标模型；(c)考虑与实际观测的Lg波谱一致的观测系统；(d)利用由不同台网管理机构给定的地震震级估计源函数；(e)生成合成的Lg谱；(f)对合成的Lg谱施以5％的随机扰动后，作为观测谱进行成像反演，重建区域Q值分布及震源函数；(g)计算区域Q值异常与实际扰动量的比值，重建的源函数与给定值的比值，比值为1时分辨率最高。

在0.05至10.0Hz之间的58个频率点上进行LgQ值成像的分辨率分析。由于在不同频率平面上的射线覆盖与观测系统的不一致，空间分辨率是不同的。图1-8a显示了每个频率的射线数。也就是说，采用统一的网格进行区域剖分是很困难的。幸运的是，从图1-7a我们可以发现，即使射线稀疏，对某一给定的网格节点，仍然有一些节点上的Q值能够被很好地重建。这样，我们就可以采用统一的1°×1°网格剖分，并获取节点Q值及相应的分辨率。根据节点Q值重建分辨率的大小进行Q值取舍。在本项研究中，节点Q值分辨率低于0.25则舍弃该节点Q值。然而，这样的作法会引起一些节点上的分辨率过低而没有数值。可以采用双曲线性插值的方法进行Q值补充。图1-8b显示了1HzQ值重建的迭代过程。

四、成像结果

中国东北及邻近地区地震Lg波Q值层析成像的结果包括在58个参考频率点上的区域Lg波Q值分布和125个Lg波震源谱函数估计。

区域Lg波Q值分布

中国东北及邻近地区是西太平洋构造带、北美大陆板块、中亚造山带、华北板块和扬子板块的交汇部位。中-新生代以来经历了多次构造事件的叠加，诸如南北大陆的汇聚碰撞、东部地体拼贴、晚中生代的大陆岩石圈减薄的伸展作用、大型走滑断裂的错移、东部陆缘的俯冲增生、新生代太平洋板块运动的转向及东亚西太平洋裂谷系的形成等，构造面貌比较复杂。谢鸣谦用“中国东北拼贴板块”解释微板块的广泛分布。这些微板块有着相同的基底，在演化过程中保持长期稳定的洲际板块特征。地槽和褶皱带位于微板块之间，且附属于微板块运动。图1-9显示了该地区的区域大地构造分布图。

图1-10显示了三个频率点上的含有主要断裂、地体缝合线和盆地边界的Q值分布图像，以及相应的射线路径分布与分辨率。随着频率的增加，Q值增大。在三个频率平面上，Q值的分布特征是一致的。从全部58个Q值图像，可以看出主要的区域特征包括：(1)盆地表现为低Q值，相对的高Q值对应于火山岩地区；(2)区域断裂与Q值的剧烈变化有关系，也就是说，活动断裂通常位于Q值梯度大的地方。例如，剡庐断裂恰位于较强的Q值变化带上，在其东边是长白山脉(CMs)，具有较高的Q值分布，而在其西侧是松辽盆地(B2)，低Q值；(3)在某些独立的微板块上存在着相对的Q值变化。典型地，对于微板块II8，与其板块均值比较，中间为高Q值，而两端为低Q值。

为了更加准确地描绘区域Lg波Q值特征，对不同的地质单元进行全频率的Q值统计，包括整个探测区域，大型板块，盆地，山脉，微板块及褶皱带等26个不同的区域。图1-11展示了宽频带的Q值的统计结果，每个小的框图为一个独立的地质单元，单元的编号标示于右下角。图中浅灰色的十字丝表示Q值，而统计的Q值用黑色实心圈表示，并给出了误差棒。统计的Q值随频率的增大而增大。对Q值进行统计计算，首先选定参考频率点，在频带内计算算术平均值。在每个框图内还给出了Q0值(在1Hz的统计计算的Q值)，这一数值用来刻划该地质单元对Lg波的衰减。图1-12展示了全部26个地质单元的Q值随频率变化。

图1-10典型的Lg波Q值分布及其分辨率。图中标出了参考频率。Q值分布图像采用统一的色标，列于最下端。

中国东北及邻近地区的平均Q0值为429。Q0值的变化范围较大，从150至800。在对数域低频Q值较高频Q值更为离散。在0.1至10Hz范围内，平均Q值随频率的增加而增大。当频率超过1Hz，Q值增加较快。在1Hz附近出现了一个LgQ值的拐点。现有技术中调查了欧亚大陆的Lg波Q值并给出横向的Q0值图像，对于我们的探测区域，他们得到的Q0均值约为450，与我们的结果很接近。

按大型板块划分，研究区可分为西伯利亚板块(I)，中国东北拼贴板块(II)，和华北克拉通(III)。统计的平均Q0值分别为544±174，459±176，和377±175，具有相对较大的标准偏差，说明Q值数据较为离散，LgQ值不能有效地约束一阶地质板块单元。

从地形地貌的角度，有6个沉积盆地(B1～6)，大兴安岭(DMs)和长白山脉(CMs)分布于研究区内。由于沉积层对Lg波的强吸收影响，盆地表现为低Q值。渤海湾盆地(B5)的Q值最低，Q0值为261±131。1Hz拐点非常明显，高频Q值快速增大。其它盆地也表现出清晰的1Hz拐点，以及Q值与频率在对数域的非线性关系。这一现象表明，传统的频率依赖的Q值模型Q(f)＝Q₀f^η，在一个较宽的频带范围内是不合适的。这与传统Q值模型在一个较窄的频带与频率呈线性关系并不矛盾。B5位于华北克拉通的东部。因此，要解释低Q值现象，除了地表沉积层的影响外，还要考虑华北岩石圈的减薄机制。大兴安岭DMs和长白山CMs的Q值较高，统计的Q0值分别为649±114and603±190。

中国东北拼贴板块由14个地块构成，包括7个褶皱带和7个微板块。除了那些包含或部分包含盆地的地质块体外，独立的拼贴地块上统计的Q0值在500左右。两个微板块之间Q的统计值差异不大，支持现有技术中得出的微板块具有相似的基底结构的结论。一些微板块上的Q值存在明显的横向差异，可能是微板块受不同期次的造山运动影响所致。例如，从图1-10可以看出，微板块II8的两端Q值较低，而中间Q值较高。这可能与加里东造山带(II4)和华力西造山带(II2,II3,II6)有关系。在横向的Q值的数值上，不能区分微板块与褶皱带，这一结果可支持现有技术中的结论，微板块之间的褶皱带是微板块运动的表现形式，从属并受制于微板块的运动。

现有技术中调查了欧亚大陆的横向地壳Lg波Q值变化。如果只考虑中国东北及邻近地区，他们得到的Q0值分布范围是250至800。平均Q0值在450左右，与我们的结果非常接近。然而，现有技术中得到频率依赖的Q值存在一个偏差。现有技术中检测了594个η测量值的分布，得到η值的范围是-0.2至0.5，平均值约为0.18。与现有技术中得到的平均值为0.4～0.5的结果存在一个系统的偏差，现有技术中认为造成这一系统偏差的原因可能是他们采用的测量方法不同。从我们的结果上看，导致这一系统偏差的真实原因是他们采用数据不同。由图1-12可以明显看出，低频段Q值变化平缓，高频段Q值变化较快。现有技术中采用的是Lg尾波Q数据，其频率成分较高，估计在0.8～5.0Hz范围内。而现有技术中所用的数据频率范围是0.1～2.0。那么，这一偏差就显而易见了。现有技术中调查了南墨西哥Lg波传播的横向变化。他们的发现之一，区域Q值模型具有一个1Hz的吸收峰，也同样支持频率依赖的Q值模型在一个较宽的频带内的非常数指数属性。

Lg波源谱函数

作为成像结果的一部分，得到了125个事件的Lg波源谱函数。对于一个独立的事件而言，由于数据选择方案的限制，如信噪比门限值，单记录事件丢弃等，不一定能够得到全部58个参考频率点上的Lg波源谱函数。但是，利用ω²源模型(方程3-23和3-24)仍然能够拟合有限的源函数得到标量地震矩M₀和拐角频率f_c。在拟合过程中，采用Bootstrap方法估计标准差。图1-13是全部125个事件的Lg波源谱函数及拟合曲线，其中十字丝表示源函数，浅灰色的曲线为由Bootstrap方法得到源模型。首先，对成像反演得到的Lg波源谱函数施以20％的随机扰动；然后，采用快速模拟退火算法获得最优的源模型，其中优化方法的目标函数为增加扰动的源谱函数与理论值的残差之和；最后，经过50次的重复操作后得到源参数及其标准差。作为结果的地震矩和拐角频率标注于每一个事件框图内。由图1-13可以看出，成像反演的Lg波源谱函数对ω²源模型参数提供很好的约束。

图1-14把全部125个事件的源函数绘于一张图上，成像反演结果的源函数与最佳拟合的源模型进行比较，其中同一事件的源函数曲线在颜色上是相同的。合成的源模型与成像反演结果非常一致，特别是在0.05至5.0Hz的频带内。从图1-14中可以看出，对于一个特定事件，随频率增加源函数表现为稳定的减小趋势。文中的成像方法为频率独立的Lg波Q值成像，在进行源函数估计时，任意两个频率之间没有施以任何约束。由此可以推断，成像结果是准确可靠的。反演结果列于表3-2。

为了验证成像结果的准确性，用方程(3-25)讨论标量地震矩M₀和拐角频率f_c之间的关系，同时检测M₀和全球震级，包括体波震级m_b，面波震级M_s和近震震级m_L。如果假定应力降Δσ是一个常数，通过结性回归，可以得到M₀～f_c的关系式

log₁₀(M₀)＝15.01(±0.05)-2.79(±0.17)log₁₀(f_c), (3-28)

其中地震矩M₀的单位是Nm。由(3-28)式可以看出，中国东北及邻近地区地震事件的源参数M₀和f_c的指数率为-2.79，与理论上预测值-3非常接近。Taylor et al.[2002]建议采用一个非常数的应力降来解释指数率的偏差。如果令v_s＝2.7km/s，方程(3-25)可以变形为：

log₁₀(M₀)＝9.71+log₁₀(Δσ)-3log₁₀(f_c), (3-29)

其中Δσ＝1.89bar，为这一地区的平均应力降。如图1-15所示，应力降不同，但不依赖于地震矩。

检测标量地震矩M₀与全球震级，体波震级m_b，面波震级M_s和近震震级m_L的关系，通常对M₀取对数后，分别采用非常系数和常系数进行线性拟合，得：

log₁₀(M₀)＝10.28(±0.23)+1.17(±0.05)m_b,sd＝0.40,r＝0.85 (3-30a)

log₁₀(M₀)＝11.02(±0.27)+1.00(±0.06)m_b,sd＝0.47,r＝0.79 (3-30b)

log₁₀(M₀)＝9.01(±0.01)+1.57(±0.00)M_s,sd＝0.02,r＝0.99 (3-31a)

log₁₀(M₀)＝9.29(±0.11)+1.50(±0.02)M_s,sd＝0.18,r＝0.98 (3-31b)

log₁₀(M₀)＝11.12(±0.47)+0.95(±0.11)m_L,sd＝0.63,r＝0.67 (3-32a)

log₁₀(M₀)＝10.89(±0.44)+1.00(±0.11)m_L,sd＝0.61,r＝0.69 (3-32b)

其中±之后为相应的标准差，sd表示线性拟合的标准偏差，r为相似系数。图1-16显示了标量地震矩与全球震级的线性关系。很明显，log₁₀(M₀)和M_s之间存在着良好的线性关系。这一结果意味着成像反演的Q值模型可以用于这一地区的地震尺度的测量，而成像反演的标量地震矩可以作为该区域地震事件直接的量度。

讨论与推论

基于大量的区域地震资料，建立了中国东北地区分辨率为1°×1°的地壳Lg波Q值成像模型。研究区平均的Q0值为429，由南向北呈增大的趋势。Q0值分布的相对低值与盆地沉积层相对应。活动构造恰好处于Q0值梯度较大的位置。Q值的分布与构造的活动性具有良好的对应关系。文中给出的频率独立的Lg波层析成像是对区域Q值和Lg波源谱函数同时反演的方法。虽然在任意两个频率间没有施以任何约束条件，对于同一事件，所得源函数表现出一致的随频率的变化规律。用ω²源模型拟合成像反演的源函数，得到标量地震矩M₀和拐角频率f_c。M₀与f_c，与各个全球震级之间存在良好的线性关系可以说明文中给出的LgQ成像方案是可信的，在中国东北及邻近地区的成像结果是可靠的。

在中国东北及邻近地区，火山岩分布广泛，典型的区域为大兴安岭和长白山，均具有较高的Q值；在盆地区域则表现为低Q值。早在中生代的火山岩喷发后冷却，以及沉积层能够对Q值变化提供直接的地质解释。However,Lg波的强烈衰减还与深部构活动性相联系。这对研究中国东部岩石圈减薄具有重要意义。西太平洋板块的俯冲与若干块体拼贴可能是形成火山岩带的构造环境.板块俯冲作用使岩石圈底部发生弯曲，向下弯曲的地方由于地幔对流作用发生拆沉作用，而后软流圈上涌使岩石圈拉张而形成盆地系。考虑到地幔对流可能存在不均匀性，使拆沉部分下沉进入软流圈而导致又一个软流圈上涌高峰，从而形成断陷盆地。无论是岩石圈物质沉入软流圈，还是软流圈物质上升与壳-幔物质作用，均会产生大量的热并上涌。这样对盆地底部的加热，使深部构造的活动性加强，从而导致低的Lg Q值。现有技术中观测到东北亚地区的低速且强烈衰减的Pn，Sn。这表明壳幔之间强烈的相互作用，暗示了下地壳与上地幔顶部的物质可能处于半熔融状态。深部构造发育，构成了热液上升通道，在中地壳表现为存在低速高导层，在地表则表现为大地热流值较高，温泉发育。

在我们的成像方法中，忽略了震源辐射花样、台基效应和传播过程中产生的随机误差，都可能限制Q值估计的分辨率。因此，我们的成像结果的Q值对火山通道，热点机制等没有直接的显示。通过对不同尺度的地质块体进行统计分析，发现Q值与频率的关系在宽频带范围内不再是对数域的线性关系。但在某一窄的频段内，如0.05-1.5Hz,0.5-2.Hz,1.6-8.0Hz,Q值与频率仍然表现为对数域的线性关系。我们已经确定的中国东北及邻近地区的Lg波Q值分布，有助于提高该地区m_b(Lg)测定和核爆当量估计精度。

Claims (1)

1.基于单台、双台和双事件数据二维Lg波Q值层析成像方法，其特征在于，包括如下步骤：

(a)地震资料预处理，获得Lg波真振幅和信噪比，设定信噪比的门限值，选择用于Lg波Q值层析成像的Lg波数据体；

地震资料预处理分为4个步骤：

(a1)分离Lg波与P波波前噪声；

为了识别Lg波并选择有效的Lg波时间窗，使用递归的零相位四阶Butterworth滤波器对地震记录进行滤波处理，其拐角频率设定为1Hz，从而得到增强的Lg波震相以便于确定其群速度窗口；或者在宽频带地震图上卷积一个全球范围标准地震台网的短周期仪器响应，并通过地震图可视化来选择Lg波群速度窗口；因为全球范围标准地震台网仪器的短周期仪器响应的高通属性，可识别的离散的Lg波在全球范围标准地震台网地震图上清晰可见；通过对全部地震记录的可视化检测，可以确定Lg波的群速度窗口为3.6-3.0或3.5-2.9km/s；噪声在P波初至前与Lg波等时间窗内提取；

(a2)计算Lg波和噪声的Fourier振幅谱；

确定Lg波和噪声的时间窗，在宽频带地震图上获得有效Lg波的持续时间；在有效Lg波两端各增加时间窗长度的10％用于余弦镶边消减；对镶边的信号进行快速Fourier变换，得到Lg波和噪声的振幅谱；

(a3)拾取振幅谱；

在0.05至10.0Hz的频带范围内选取参考频率，并使之在对数域中为等间隔，即，Δlog₁₀f_ref＝0.04；对于每个参考频率f_ref，与其对应的Lg波振幅谱A(f_ref)采用均方根的方法来求取，即，其中f_t∈[f_ref-Δf,f_ref+Δf]，这里的Δf可由Δlog₁₀(Δf)＝0.02来确定，M是采样点的总数；

(a4)去噪处理；

假定Lg波时间窗内的波列为Lg波与噪声的能量叠加，并认为Lg波与噪声的振幅谱是不相关的，利用下式估计每个频率上的Lg波振幅谱：

$A_{sig}^2\left(f_{ref}\right)=A_{obs}^2\left(f_{ref}\right)-A_{noi}^2\left(f_{ref}\right)---(3-1)$

其中A表示振幅谱，下标sig、obs和noi分别指真振幅、观测值以及噪声；

(b)从单台数据中提取双台和双事件数据；

(b1)单台数据：

给定一个频率f，对于事件k在台站i观测到的Lg波振幅谱A_ki可以表示为：

A_ki＝S_kG_kiΓ_kiP_iR_ki (3-2b)

S_k是事件k的震源函数，G_ki＝(Δ₀Δ_ki)^-1/2，Δ_ki≥100km是几何扩展因子,Δ_ki是事件k和台站i的震中距离，并且Δ₀是固定在100km的基准距离；衰减项Γ_ki可以被表示为：

${\mathrm{\Γ}}_{ki}=\exp \[-\frac{\mathrm{\π}f}{V}\·{\∫}_k^i\frac{ds}{Q(x,y,f)}\]---(3-3b)$

V是Lg波的群速度，是事件k到台站i沿大路径的积分，Q(x,y,f)是Lg波质量因子，它是频率与陆地位置(x,y)的函数；P_i是当地地球结构控制台站i响应放大率的乘法系数；R_ki是Lg波在事件k和台站i之间传播的随机效应，忽略不计；

(b2)双台数据：

两台站记录同样的地震，它们在一条线上，双台站振幅谱可通过计算标量谱比获得；理想情况下，两台站i和j与事件k对齐，双台站振幅谱可以准确的通过下式计算：

$A_{ij}=\frac{A_{kj}}{A_{ki}}={\left(\frac{{\mathrm{\Δ}}_{kj}}{{\mathrm{\Δ}}_{ki}}\right)}^{-1/2}\·\exp \[-\frac{\mathrm{\π}f}{V}\·\left({\∫}_i^j\frac{ds}{Q(x,y,f)}\right)\]\·\frac{P_j}{P_i}---(3-4b)$

A_kj是记录事件k在j台站的振幅谱，Δ_kj是事件k和台站j的震中距离，P_j是当地地球结构控制台站j响应放大率的乘法系数；实际情况下，两台站也许没有对齐完美，并存在一定的公差；使用两台站的方位角差，设置最大15°的公差；并设置一个参考点l在射线路径kj上，让距离kl等于距离ki，并且需要使台站i和参考点l之间的距离在基于频谱分辨率上的半个网格空间之内；当台站符合这一标准，计算出台站间振幅谱：

$A_{lj}\≈\frac{A_{kj}}{A_{ki}}={\left(\frac{{\mathrm{\Δ}}_{kj}}{{\mathrm{\Δ}}_{kl}}\right)}^{-1/2}\·\exp \[-\frac{\mathrm{\π}f}{V}\·\left({\∫}_l^j\frac{ds}{Q(x,y,f)}\right)\]\·\frac{P_j}{P_i}---(3-5b)$

Δ_kl是事件k和参考点l之间的距离，基于假设，有相似的传播路径和台站，所以在台站i和参考点l观测的振幅谱是相同的；

(b3)双事件数据：

类似双台振幅，当台站记录两个地震，并且它们在一条线上，双事件振幅谱可通过计算标量频谱比获得；台站j与事件k和事件h对齐，双事件振幅谱可以表示为：

$A_{hk}=\frac{A_{kj}}{A_{hj}}=\frac{S_k}{S_h}\·{\left(\frac{{\mathrm{\Δ}}_{kj}}{{\mathrm{\Δ}}_{hj}}\right)}^{-1/2}\·\exp \[-\frac{\mathrm{\π}f}{V}\·\left({\∫}_h^k\frac{ds}{Q(x,y,f)}\right)\]---(3-6b)$

A_hj是台站j对于事件h记录的振幅谱，S_h是事件h的震源函数，Δ_hj是事件h和台站j的震中距离；实际情况下，两个事件并不完全对齐；设置一个参考点q在射线路径kj上，让距离qj等于距离hj；当事件h和参考点q之间的距离不超过半个网格空间，双事件振幅谱可以通过如下近似计算：

$A_{qk}\≈\frac{A_{kj}}{A_{hj}}=\frac{S_k}{S_h}\·{\left(\frac{{\mathrm{\Δ}}_{kj}}{{\mathrm{\Δ}}_{hj}}\right)}^{-1/2}\·\exp \[-\frac{\mathrm{\π}f}{V}\·\left({\∫}_q^k\frac{ds}{Q(x,y,f)}\right)\]---(3-7b)$

基于假设，考虑到相同的震源函数和相似的传播路径，事件h和参考点q在台站j观察到的振幅谱是相同的；

(c)基于振幅谱的Lg波Q值层析成像；

将用于速度结构成像的标准的走时层析技术扩展到Lg波Q值结构重建，具体方案包括：

(c1)假定Lg波在地壳内沿震源至台站的大圆路径传播，将Lg波Q值看成位置与频率的函数，Q(x,y,f)，其中x和y分别表示离散区域的经度和纬度，f为给定频率；

(c2)根据先验地质信息给定初始Q值模型，任意的给定频率的震源函数的初始值设定为1.0；

(c3)基于给定的初始Q值模型和震源函数，合成Lg波振幅谱；将合成的Lg波振幅谱与观测的Lg波振幅谱进行比较并计算残差，以此建立方程组求取Q值与震源函数的修正量；

(c4)对Q值模型和震源函数进行修正；

(c5)进入下一步线性化并返回步骤(c3)，直到满足某一收敛条件后终止迭代；

Lg波Q值层析成像的详细过程包括：

地震Lg波振幅谱可以表示为：

A(f,Δ)＝S(f)R(θ)G(Δ)Γ(f,Δ)X(f)r(f) (3-2)

其中f为频率，Δ是震中距，A(f,Δ)为Lg波振幅谱，震源函数S(f)可表示为：

$S\left(f\right)=\frac{M_0}{4{\mathrm{\π\ρv}}_s^3}{S}^{\′}\left(f\right)---(3-3)$

其中S′(f)是归一化的震源函数，M₀是地震矩，ρ和v_s是震源区的平均密度值和剪切波的平均速度；

G(Δ)＝(Δ₀Δ)^-1/2 (3-4)

是几何扩展因子，Δ₀是参考距离，设定为100km，R(θ)是震源辐射花样，随方位角θ变化而变化；X(f)和r(f)分别为与频率相关的台站基响应函数和Lg波传播过程中的随机效应；假定X(f)＝r(f)＝1，Lg波对方向的影响很弱，在实际计算中令R(θ)＝1；Γ(f,Δ)为衰减的函数，可以表示为：

$\mathrm{\Γ}(f,\mathrm{\Δ})=\exp \[-\frac{\mathrm{\π}f}{V}\·B(f,\mathrm{\Δ})\]---(3-5)$

其中V是Lg波的群速度，B(f,Δ)是Lg波Q值沿大圆路径的积分，

$B(f,\mathrm{\Δ})={\∫}_{ray}\frac{ds}{Q(x,y,f)}---(3-6)$

其中Q(x,y,f)为Lg波Q值，它是频率与陆地位置(x,y)的函数；在衰减成像过程中，积分(3-6)式可以离散为：

$B\left(f\right)=\underset{z=1}{\overset{Z}{\Σ}}{\∫}_z\frac{ds}{Q(x,y,f)}---(3-7)$

其中z指第z个网格单元内的射线路径，Z是射线经过网格单元的总数，∫_zds为该单元内的积分；

对方程(3-2)式两端取自然对数，可得：

$ln\[A(f,\mathrm{\Δ})\]-ln\[G\left(\mathrm{\Δ}\right)\]=ln\[S\left(f\right)\]-\frac{\mathrm{\π}f}{V}\·B(f,\mathrm{\Δ})---(3-8)$

令

$\frac{1}{Q(x,y,f)}\≈\frac{1}{Q^0(x,y,f)}-\frac{\mathrm{\δ}Q(x,y,f)}{{\[Q^0(x,y,f)\]}^2}---(3-9)$

ln[S(f)]＝ln[S⁰(f)]+δln[S(f)] (3-10)

并代入(3-7)和(3-8)式可得：

$ln\[A(f,\mathrm{\Δ})\]-ln\[G\left(\mathrm{\Δ}\right)\]-ln\[S^0\left(f\right)\]+\frac{\mathrm{\π}f}{V}\·B^0\left(f\right)=\mathrm{\δ}ln\[S\left(f\right)\]-\frac{\mathrm{\π}f}{V}\·\mathrm{\δ}B\left(f\right)---(3-11)$

其中Q⁰表示初始Q值模型，S⁰表示反演过程中迭代后所得的震源函数，B⁰是将Q⁰代入式(3-6)中得到的中间函数，δB(f)可以表示为：

$\mathrm{\δ}B\left(f\right)=\underset{z=1}{\overset{Z}{\Σ}}{\∫}_z\frac{\mathrm{\δ}Q(x,y)}{{\[Q^0(x,y)\]}^2}ds---(3-12)$

其中z指第z个网格单元内的射线路径，Z是射线经过网格单元的总数，采用矩形网格化的方法，在每一个矩形单元内有四个角点上的Q值与通过它的射线相关；矩形单元内任意点上的Q值可以用双曲线函数表示：

Q(x,y)＝a₀+a₁x+a₂y+a₃xy (3-13)

其中双曲线系数a_p-1，p＝1,2,3,4把单元内Q值与角点Q值联系在一起，可以通过下式计算：

$a_{p-1}=\underset{w=1}{\overset{W}{\Σ}}{(-1)}^{p+1}\frac{{\mathrm{\Δ}}_{wp}^{\′}}{{\mathrm{\Δ}}^{\′}}{Q}_w,p=1,2,3,4---(3-14)$

其中W＝4，是矩形网格单元的节点数；Q_w是第w个节点上的Q值，Δ′_wp是按Cramer法则解4阶线性方程组的4阶行列式Δ′的代数余子式；按微分的定义，可得：

$\mathrm{\δ}Q(x,y)=\underset{w=1}{\overset{W}{\Σ}}\frac{\∂Q(x,y)}{\∂Q_w}{\mathrm{\δQ}}_w---(3-15)$

把方程(3-13)代入(3-15)式，可得：

$\mathrm{\δ}Q(x,y)=\underset{w=1}{\overset{W}{\Σ}}(\frac{\∂a_0}{\∂Q_w}+\frac{\∂a_1}{\∂Q_w}x+\frac{\∂a_2}{\∂Q_w}y+\frac{\∂a_3}{\∂Q_w}xy){\mathrm{\δQ}}_w---(3-16)$

其中

$\frac{\∂a_{p-1}}{\∂Q_w}={(-1)}^{p+1}\frac{{\mathrm{\Δ}}_{wp}^{\′}}{{\mathrm{\Δ}}^{\′}},p=1,2,3,4---(3-17)$

把(3-16)式代入(3-12)式，得：

$\begin{array}{c}\mathrm{\δ}B\left(f\right)=\underset{z=1}{\overset{Z}{\Σ}}\underset{w=1}{\overset{W}{\Σ}}(\frac{\∂a_0}{\∂Q_w}{\∫}_z\frac{1}{{\[Q^0(x,y,f)\]}^2}ds+\frac{\∂a_1}{\∂Q_w}{\∫}_z\frac{x}{{\[Q^0(x,y,f)\]}^2}ds\\ +\frac{\∂a_2}{\∂Q_w}{\∫}_z\frac{y}{{\[Q^0(x,y,f)\]}^2}ds+\frac{\∂a_3}{\∂Q_w}{\∫}_z\frac{xy}{{\[Q^0(x,y,f)\]}^2}ds){\mathrm{\δQ}}_w\end{array}---(3-18)$

对于多台多事件，由方程(3-11)和(3-18)能够建立线性系统：

H＝A·δQ+E·δU (3-19)

基于单台、双台、双事件可建立如下线性系统：

$\left[\begin{array}{c}{H}_{1s}\\ H_{2s}\\ H_{2e}\end{array}\right]=\left[\begin{array}{c}{E}_{1s}\\ 0\\ E_{2e}\end{array}\right]\·\mathrm{\δ}U+\left[\begin{array}{c}{A}_{1s}\\ A_{2s}\\ A_{2e}\end{array}\right]\·\mathrm{\δ}Q+\left[\begin{array}{c}{F}_{1s}\\ F_{2s}\\ 0\end{array}\right]\·\mathrm{\δ}P---(3-19a)$

H是一个由观测的Lg波振幅谱和合成Lg波振幅谱间的残差组成的向量，下标1s,2s,2e分别代表单台站，双台站和双事件数据集，δU是一个有震源项扰动组成的向量，基数E是代表震源扰动和观测Lg波振幅谱的关系，δQ是Q值模型扰动组成的向量，基数A代表Q值模型扰动和观测Lg波振幅谱间的关系，δP是台站响应扰动组成的向量，基数F代表台站响应扰动和观测Lg波振幅谱间的关系；采用两步法去解决反问题；这样，台项在等式(3-19a)中可以表达为：

$\left[\begin{array}{c}{F}_{1s}\\ F_{2s}\\ 0\end{array}\right]\·\mathrm{\δ}P=0---(3-19b)$

其中H的元素h_u为：

$h_u=ln\[A_u(f,\mathrm{\Δ})\]-ln\[G_u\left(\mathrm{\Δ}\right)\]-ln\[S_z^0\left(f\right)\]+\frac{\mathrm{\π}f}{V}\·B_u^0\left(f\right)---(3-20)$

是观测的Lg波振幅谱与合成Lg波振幅谱的残差，公式(3-20)中的下标u和z分别表示第u条射线和第z个地震事件，

A的元素a_tu为：

$a_{tu}=-\frac{\mathrm{\π}f}{V}\·\underset{n=1}{\overset{N_{tu}}{\Σ}}(\frac{\∂a_0}{\∂Q_t}{\∫}_n\frac{1}{{\[Q^0\]}^2}ds+\frac{\∂a_1}{\∂Q_t}{\∫}_n\frac{x}{{\[Q^0\]}^2}ds+\frac{\∂a_2}{\∂Q_t}{\∫}_n\frac{y}{{\[Q^0\]}^2}ds+\frac{\∂a_3}{\∂Q_t}{\∫}_n\frac{xy}{{\[Q^0\]}^2}ds)---(3-21)$

a_tu是微分系数，公式(3-21)中t，u，z分别指统一编号网格单元节点的序号，射线序号，震源函数的序号，n指第n个网格单元内的射线路径，N是射线经过网格单元的总数；Q⁰是Q⁰(x,y,f)的简化表达；δQ_t指第t个节点上的Q值扰动量；E是一个(0,1)矩阵，当第u条射线与第z个震源函数有关时，矩阵E的元素e_zu＝1；δu_z＝δ[ln(S_z(f)]为震源的扰动量；用LSQR算法求解线性方程(3-19a)和(3-19b)式，迭代4～6次即可收敛；

Lg波Q值层析成像的迭代收敛条件采用如下形式：

$E\left(iter\right)=\frac{1}{M_1}\underset{u=1}{\overset{M_1}{\Σ}}|A_u{\left(f\right)}_{obs}-A_u{\left(f\right)}_{syn}^{iter}|\→min---(3-22)$

其中E(iter)指第iter次迭代后的系统误差，M₁是射线方程总数，A_u(f)_obs和分别指第u条射线方程中观测的Lg波振幅谱和第iter次迭代后的合成Lg波振幅谱。