0 引言

场地条件是影响震害分布的重要因素，而场地浅层波速是描述工程场地条件最基本的定量指标.常用的确定场地浅层速度结构的方法有直接的钻孔测井法和间接的主动源地震勘探方法.钻孔测井法费用昂贵，探测深度有限；地震勘探方法易受场地条件限制而难以展开测线，分辨能力也有限.需要开发更简便经济的浅层波速结构探测手段.

利用水平与竖向谱比(Horizontal-to-Vertical Spectral Ratio, HVSR)方法反演场地速度结构是当前国际上正迅速兴起的研究领域，国内外许多学者在这方面开展了研究工作，但已有的研究中多采用噪声作为数据源进行HVSR反演.例如Herak(2008)基于环境噪声HVSR反演方法，发展了一套Matlab程序(ModelHVSR)，并将其用于反演场地浅层速度结构.张立等(2009)认为地脉动以基阶瑞利面波为主，提出了单个三分量检波器地脉动信号的HVSR反演剪切波速结构的遗传算法.王伟君等(2009, 2011)的研究表明噪声HVSR方法不仅可以给出土层卓越频率等场地响应特征，还可以探测近地表土层的界面起伏，有助于识别活动断层的近地表位置及其活动特征，但该方法容易受到风、地形、附近瞬态震源等因素影响.Bignardi等(2016)进一步扩展了Herak的工作，将噪声HVSR方法用于二维和三维场地速度模型的联合反演.以上这些研究对Nakamura(1989)的原始HVSR表达式进行了简单推广应用.近十余年来，HVSR方法的理论和应用得到了极大的发展.Sánchez-Sesma等(2011)和Kawase等(2011)的研究表明了HVSR与场地的S波和P波传递函数及其基岩S波和P波波速密切相关，给出了基于散射场理论的HVSR计算公式，夯实了HVSR反演方法的理论基础.场地速度HVSR反演不局限于地震噪声，已扩展至地震观测数据(Kawase et al., 2011).与噪声数据类似，地震观测数据(如速度、加速度记录)也是一种可以广泛利用的数据类型.Nagashima等(2014)应用散射场理论对多个K-NET台站场地速度结构进行了反演，反演结果均能很好地重现观测地震动，在较宽的频带范围内很好地解释了不同强度地震作用下观测HVSR曲线的差异.Rong等(2016)利用中国西部强震观测网记录到的小震HVSR反演了若干个台站场地的一维等效模型，该等效模型揭示的场地效应与强震观测具有很好的一致性.但总的来说，利用单台地震记录的HVSR反演场地浅层速度结构仍存在两个方面的问题：(1)不同学者针对不同数据源给出的HVSR理论公式存在一个系数差别，采用何种公式、何种类型的地震波进行反演最为有效仍需观测记录的验证；(2) HVSR反演是一个土层场地模型空间搜索的全局优化问题，目前已有的蒙特卡罗搜索算法等计算效率较低，需要发展高效的反演方法.

针对第一方面的问题, 笔者已另撰文予以讨论(荣棉水等，2016；Rong et al., 2017)，本文仅对HVSR方法反演场地速度结构的理论依据作简要介绍，主要针对第二方面的问题开展研究工作.当模型搜索空间的复杂程度增大，蒙特卡罗方法收敛速度很慢，误差概率性等问题突出.相比较而言，遗传算法(GA)是一种更智能的自适应全局优化方法，其概率搜索高效，在波速结构反演中得到广泛应用(Yamanaka and Ishida, 1996；师黎静，2007).但大量的实际应用表明遗传算法存在早熟和进化停滞等缺陷.引入其他优化方法对遗传算法进行改进是非常有效的技术途径，如崔建文(2004)引入模拟退火(SA)发展了一种混合SAGA算法，SA和GA均是模拟自然法则实现优化计算的广义Monte Carlo法，由于引入了优胜劣汰机制，且系统总是趋向于平衡和稳定，这两种算法的计算效率要高于传统的Mente Carlo法.一般而言，GA算法具有隐并行性，全局搜索能力强，但容易收敛到局部最优解，即收敛早熟；SA算法局部搜索能力强，但全局搜索能力却较差，二者的结合能取长补短极大提高反演的效率与准确性.许多研究者将混合的遗传模拟退火算法应用于各个学科领域，如用于电路优化设计(刘志刚等，2004)、结构破坏模式优化(刘文政和叶继红，2013)、军事火力打击(吴坤鸿和詹世贤，2016)、通信(崔小勇和林宁，2016)等，但鲜有研究将该混合方法应用于S波HVSR反演.为此，本文构造了分层场地S波HVSR反演的数学模型，交替进行GA和SA操作，实现了一种结合遗传和模拟退火算法优点的HVSR反演混合全局优化算法，较好地解决了蒙特卡罗算法搜索效率低的问题，将这一算法应用到单台加速度记录的HVSR反演中获得了很好的反演效果.

考虑到HVSR方法采用的是水平与竖向分量傅氏谱的比值，使用加速度记录和速度记录在理论上能得到一致的结果，这就为充分利用强震观测的加速度记录进行反演提供了新的途径.随着强震观测技术的发展，我国已经积累了极为丰富的强震观测资料，国际上也已拥有海量强震观测数据，但利用地震加速度记录探测浅地层速度结构的研究仍相当少，本文在研究中采用加速度记录作为数据源.

1 S波HVSR反演场地浅层速度结构的理论依据

根据Nakamura(1989)的定义，S波HVSR_S定义为地表水平分量傅氏谱H_S(f)与竖向分量傅氏谱V_S(f) (如图 1)的谱比：

(1)

式(1)可改写为如下形式：

(2)

式中H_B(f)与V_B(f)分别为基岩处水平分量和竖向分量的傅氏谱(图 1).假设：

(3)

(4)

(5)

ETF_H、ETF_V分别为土层水平方向和竖向的经验传递函数，HVSR_B为基岩处的水平与竖向谱比，在场地下行波效应可忽略的情况下，ETF_H、ETF_V可分别写为场地的水平和竖向传递函数TFH、TFV(Rong et al., 2017)，则式(2)可改写为如下形式，即HVSR_S等于HVSR_B乘以水平与竖向传递函数的比值：

(6)

若假定地震动为垂直入射的体波，则水平向传递函数TF_H可视为S波传递函数，竖向传递函数TF_V可视为P波传递函数，垂直入射假定看似具有很大局限性，实际上在土层较深层与最浅层具有较大波速差时，是一个很好的近似(Herak，2008).此外，笔者也曾利用美国、日本多个竖向观测台阵的强震观测记录证实基岩层位的水平与竖向谱比HVSR_B可近似看作一个常数.因此，式(6)可近似表示为：

(7)

其中C为与HVSR_B有关的常数.基于散射场理论，Kawase等(2011)给出了体波HVSR的具体表达式：

(8)

其中，TF_S(f)和TF_P(f)与式(7)的含义相同，分别表示S波和P波传递函数，α、β分别为基岩处P、S波波速.根据以上论述可知，无论是直接从HVSR定义式的推导，还是基于散射场理论推导，均能获得形式一致的公式，这也侧面证实了公式(8)的合理性.该公式得到了强震观测记录数据的验证(Kawase et al., 2011)，并应用于场地浅层速度结构反演(Nagashima et al., 2014).值得注意的是，公式(8)并无地震波垂直入射的限制，对于地震动中的S波和尾波均适用，笔者对美国和日本多个竖向观测台阵的强震观测记录的研究证实，对于地震动中的S波，基于公式(8)的理论HVSR与观测HVSR符合较好.此外，Nakamura(1989)在早期研究中认为，HVSR具有类传递函数效应，而S波传递函数一般能较好体现场地的放大效应，具有清晰的物理意义.因此，本文所指的HVSR为利用地震观测中的S波成分计算得到.式(8)中的参数仅与场地本身的特征有关，可将该式用于反演场地速度结构.

2 遗传模拟退火反演方法

HVSR曲线与场地特征密切相关，可构建计算场地速度结构的反演问题.Herak(2008)采用模型空间蒙特卡罗搜索反演了一个四层理论模型，张立等(2009)采用遗传算法反演了成都地区某较浅场地的剪切波速结构.蒙特卡罗搜索方法是利用随机数改变与土层模型相关的控制参数(如剪切波速、泊松比、密度、厚度等)，每个控制参数都是多解的，对不同的控制参数进行组合，就可以获得一个完整的模型空间.对每个模型进行正演计算，利用正演结果对目标HVSR曲线进行拟合，反演最佳模型.但当土层层数增多时，模型空间的复杂程度迅速增大，蒙特卡罗方法收敛速度慢、误差概率性等问题突出.本文构造了分层场地S波HVSR反演的数学模型，交替进行GA和SA操作，实现了一种混合SAGA全局优化算法，算法步骤如下：

(1) 对所需要解决的问题进行数学建模，确定目标函数和计算精度，确定适应度函数.设置优化参数，包括种群大小(假定含M个个体)、变量的维数、二进制位数、模拟退火初始温度.

(2) 初始化种群，计算初始种群的目标函数值.种群中每个个体是对模型参数采用二进制编码并将码串按顺序连接形成的以0和1表示的一串代码.

(3) 计算初始种群中每个个体的适应度值，进入GA的选择、重组、变异、早熟判断操作，对这些操作步骤简要介绍如下：

a，用轮盘赌方法进行选择操作

吸收文献(崔建文，2004)中对适应度函数的定义，用相对适应度表征群体中模型的相对优劣，表达式为：

(9)

其中归一化因子.定义平均适应度为：

(10)

轮盘赌方法即用相对适应度形成累计概率A_i：

(11)

随机转动轮盘，在[0，1]区间中产生一随机数ζ，若A_i-1 < ζ≤A_i，则将第i个个体选入种群，重复操作M次将产生经选择后的新种群.由于个体进入下一代的概率等于其相对适应度，相对适应度越高的个体被选中的可能性就越大，相对适应度较小者将被淘汰，体现了算法在模型空间的优胜劣汰作用.

b，重组操作：将群体中个体的二进制码串随机两两配对并随机确定一个码串的交叉位置，互换配对的两个个体在交叉位置之后的二进制串，重组成新的群体.

c，变异操作：随机确定新群体中个体的变异位置，按给定的概率将变异位置上的二进制码进行0、1互换.

d, 早熟判断：利用ψ_max-ψ_min≤η判断群体是否早熟，η可根据经验确定，本文中一般可取为0.05.ψ_max、ψ_min分别为群体中个体适应度的最大、最小值，若满足早熟条件，则进入SA操作.

(4) SA操作：对适应度满足ψ_j < ψ的个体加一随机扰动生成新个体，计算两个个体目标函数的差异ΔΦ，采用Metropolis接受准则(Metropolis et al., 1953)判断是否接受该个体.若ΔΦ < 0，接受该个体；若ΔΦ≥0，则以Boltzman接受概率p=exp(-ΔΦ/T_k)接受该个体.

本文混合优化方法可用以下流程图 2来表示，基于该算法流程，编制了matlab计算程序.

3 反演实例 3.1 理论模型的构造

为验证本文反演算法，将崔建文(2004)中的理论模型从3层增加至10层，构造一新的理论计算模型如表 1所示.由于土介质的剪切波速一般在100~500 m·s^-1范围内变化，且工程上一般将剪切波速大于500 m·s^-1的层位定义为软基岩，因此表 1中理论模型的各层剪波速均取为100~500 m·s^-1之间的值，底层剪波速取为500 m·s^-1.

采用式(8)计算一维分层模型的理论HVSR，式(8)中的理论传递函数则采用Tsai(1970)提出，并经Herak(2008)改进后的快速递归方法进行计算.该方法将土层近似为均匀半空间上覆的多个黏弹性土层，每层土的特性参数包括层厚(h)、体波传播速度(V_P、V_S)、密度(ρ)、频率相关的Q值(Q_P、Q_S)，Q值控制了土层的滞弹特性.由于剪切波速是我们特别关注的反演参数，我们将Q_S值依照Wang等(1994)给出的经验公式，取为与V_S相关的参数，本文中取Q_S=0.08×V_S，Q_P的值取为2倍Q_S，再引入泊松比，假定各层泊松比(μ)为一个常数，按照式(12)可计算出V_P.

(12)

至此，每一土层独立的土体特性参数变为(h、ρ、V_S、μ).图 3给出了初始计算模型的理论HVSR曲线.若假定土层结构的厚度、密度、泊松比已知，则土体的理论HVSR为仅与土层剪切波速度结构有关的函数.

3.2 目标函数和计算精度的确定

定义目标函数为如下形式:

(13)

其中X={x₁，x₂，x₃，…，x_n}是模型的参数向量，x_j={V_Sj, h_j, μ_j, ρ_j}是第j层介质的剪切波速、厚度、泊松比和密度，n为土层的层数，HVSR_the(f_i)为理论计算S波HVSR，HVSR_obs(f_i)为观测S波HVSR，f_i为频率，HVSR_avg(f_i)为HVSR_obs(f_i)的均值，I为计算频率点的个数.当n较大时，Φ(X)为多参数、多极值的非线性函数.将表 1中V_S向量任意增加一个代表V_S扰动的很小的随机向量(元素的值小于1)，代入式(13)得到目标函数值为5.19×10^-5，将该值作为判别反演是否结束的计算精度指标.

3.3 反演结果

假定土层厚度、密度、泊松比为已知参数，待反演的为土层剪切波速V_S，按照前文论述，Q_S、Q_P、V_P均可依据V_S求得，则待反演的V_S构成一个10维的模型空间.利用第2节所述遗传模拟退火方法进行优化反演，计算开始时的初始种群为10维速度向量，向量中每个速度值为100到500 m·s^-1之间的随机值，每代种群中个体数目为400个，模拟退火初始温度设为10 ℃，最终温度为0.5 ℃，采用比例退火策略，比例系数为0.99.图 4给出了利用不同反演方法各进行5次随机试验获得的收敛曲线，在整个反演过程中种群目标函数的均值和最小值都随着种群代数的增长快速减小.当采用传统遗传算法，遗传代数达到500代时5次随机试验中有4次未达到收敛(表 2)，即使不中止计算继续遗传迭代也无法使得种群目标函数的最小值继续减小，说明遗传算法已收敛至局部极值且无法继续向总体极值靠拢，遗传进入早熟状态.当采用本文混合方法时，5次随机试验均在遗传至约300代时其种群目标函数最小值迅速降低到小于预设计算精度值，且计算停止，计算总模型数约12万个.这说明本文混合方法相比遗传算法已有明显改进，能有效地避免遗传早熟且能快速向总体极值收敛.第一次混合反演试验计算终止时最小目标函数对应的10维剪切波速向量为[150.0978, 199.4129, 300.3914, 400.5871, 309.0020, 471.0372, 490.6067, 279.2564, 330.1370, 500.0000]，与表 1中V_S近乎完全一致，根据该反演模型得到的理论HVSR曲线示于图 3中，与反演的目标曲线(初始理论计算模型)的HVSR曲线几乎完全一致.图 5给出了不同方法得到的理论模型与反演模型的剪切波速度结构，由该图可知混合方法反演模型与理论模型已经非常吻合而传统遗传法反演模型与理论模型差别明显，其获得较精确解的概率是比较小的.此外，从图 5a可知，本文混合方法反演得到的速度模型在浅部更接近于理论模型，深部的收敛效果不如浅部.这主要是因为场地卓越频率主要受浅部土层特性的影响，而本文反演方法的实质就是寻找最佳的土层模型，使其理论HVSR曲线与观测HVSR曲线吻合最好.当给定的目标函数值足够小，反演朝设定目标函数值趋近的过程中，由于卓越频率附近的HVSR值较大，对目标函数的影响最为显著，只有卓越频率附近拟合较好的结果才有可能保留，因此浅部土层速度值相比深部更易收敛到理论值附近.

3.4 利用强震加速度记录的反演

作为方法应用于强震加速度记录反演的一个示例，本文选择了已有大量加速度记录且有台站下方钻孔剖面的日本IWTH25台站，该台隶属于KIK-NET台网.由于地表PGA能间接反映场地非线性程度，地表PGA较小时可认为场地近似处于弹性状态，本文选择了PGA较小的22次地震事件(信息示于表 3)，其PGA范围为4.0~4.9 Gal.图 6给出了一次典型地震事件地表记录的加速度及其S波窗.这些地震事件的S波HVSR具有很好的一致性，其观测平均HVSR示于图 7a中.表 4给出了台站场地初始模型信息，场地速度结构为需要反演的未知量，则整个问题化为6维速度模型空间的反演问题.将目标函数取为式(14)形式，介质的剪切波速取为100~2000 m·s^-1范围内变化的随机值，利用本文提出的混合优化方法，反演至100多代时即收敛到稳定的最小值.此时目标函数最小值对应的反演速度模型如图 7b所示.由图 7可知，利用观测得到的地表HVSR曲线及初始模型信息，能通过反演方法获得场地浅层的速度结构，其速度结构与实测初始速度模型较为接近，反演模型的理论HVSR曲线与实测初始模型的理论HVSR曲线在卓越频率和幅值上均较为接近，表明反演速度模型能较好体现场地的基本自振周期和放大效应.需要说明的是，本文用于反演的目标曲线为观测HVSR曲线，其与初始实测模型的理论HVSR仍存在一定差异(图 6a)，这就使得反演速度结构与初始实测模型必然存在差异(图 6b)，但实际上，观测HVSR曲线由于体现了场地实际观测环境，如场地地形效应等，其反演模型可视为一维等效模型，在研究不同地震动输入作用下场地的放大效应时，可能比实测模型更接近于实际情况，这一点已有研究予以揭示(Rong et al., 2016).

对于以上两个算例，我们在反演过程中均引入了几点假设，包括场地土层层数已知，土层Q值与V_S存在经验关系，土层的特性参数h、ρ、V_S、μ中只有V_S为未知参数，而实际情况下往往难以预知场地的分层特征和各土层的特性参数，因此本文反演方法离广泛实用仍有一定差距.由于本文目标为实现一种HVSR反演混合全局优化算法并验证该方法的合理性和适用性，我们主要选取了场地分层数较少的理论和实际场地.为进一步提高本文方法的实用性，我们可以采取以下技术途径：(1)分析不同参数对反演结果不确定性的影响，以对模型空间进行合理约束；(2)进一步优化HVSR反演混合全局优化算法，从设置固定土层厚度进行速度反演或同时反演土层厚度和波速两方面实现高维度模型空间的反演.这些工作留待未来进一步研究.

4 结论

本文简述了S波HVSR反演场地浅层速度结构的理论依据，实现了一种结合遗传和模拟退火算法优点的HVSR反演混合全局优化算法，较好地解决了蒙特卡罗方法收敛速度慢、遗传算法收敛早熟和模拟退火算法搜索效率低的问题.利用一理论模型的反演验证了本文方法的合理性和准确性，随后，将本文反演方法用于实际强震加速度记录观测的S波反演，也获得了较好的反演效果.本文的研究为基于单个地震台的地震观测记录反演浅层粗略速度结构提供了一种高效且较为准确的反演方法.

致谢

谨此祝贺姚振兴先生从事地球物理教学科研工作60周年.本文所使用的IWTH25台站及其强震数据资料均来源于日本防灾科学技术研究所强震观测网(http://www.kyoshin.bosai.go.jp/kyoshin/)，在此表示感谢！