2021, Vol. 51
2. 山东省海洋环境地质工程重点实验室，山东 青岛 266100
地球表面时刻都在进行着微弱的振动，它区别于地震或其它有特定震源的振动，位移振幅一般为几微米，这种微弱振动在工程地震领域中被称为微动。微动信号中含有丰富的地层内部结构信息，因此在浅层工程地震勘探中备受关注。微动由体波和面波组成，其中70%以上的微动能量都来自于面波[1]。早在1887年，英国物理学家Rayleigh在其研究中预言了这种波[2]，后来的科研工作中也证实了面波的存在，并将其称为瑞利面波[3]。20世纪50年代初，研究者们又发现瑞利面波在层状介质中具有频散特性，即在均匀水平层状介质中，其速度会随着频率变化。在此基础上Haskell首次采用矩形台阵的方法[4]得到了层状介质中瑞利面波的频散曲线，他的研究成果成为利用瑞利面波探测地层结构的理论基础。Aki和Capon分别用空间自相关法(SPAC)[5]和频率-波数(F-K)法[6]，从微动信号中提取出了面波，进而求出面波的频散曲线。
剪切波只能在土骨架中传播，可以体现土的骨架性质。正是由于剪切波所具有的这种独特性质，使其成为反映场地土动力特性以及作为勘察地质灾害的一项重要指标。因此通过微动探测获得的面波频散曲线反演土层剪切波速成为了工程物探中的一种新方法—微动探测技术。早期微动探测技术主要用于地幔结构[7]、地基土类型划分等研究，近年来，随着城市化进程的发展，传统的工程物探法由于受到噪声源多以及破坏环境等方面的制约，不易在城市及邻近区域开展。微动探测法则由于其无需震源、经济、安全、环保等优点被科研工作者尝试用到地质灾害[8]、地下资源[9]、城市物探[10]等领域，但这些方面的研究成果仍然较少, 应用于地下空间的效果仍不明朗。
目前常用的瑞利面波频散曲线反演方法有遗传算法、下山单纯形法和模拟退火法等。遗传算法[11]是一种模拟自然进化过程搜索最优解的方法，其优越性在于不依赖初始模型的选择、不容易陷入局部极小等，但其收敛较慢、运行时间较长、局部搜索能力弱；下山单纯形法(DHSM)[12]是一种并行处理的局部优化算法，具有简单、效率高等优点，但其全局搜索能力差；模拟退火法(SA)[13]是一种全局优化算法，可以获得全局最优解，但是其退火速度缓慢、计算时间长。可见单一的反演方法已经无法满足当前的需要，将不同反演方法相结合进行混合优化反演受到了广泛研究[14]。
本文采用SPAC法对试验区进行微动探测工作，通过分析数据资料，并与其它物探方法所得结果对比，综合分析了组合反演方法的效率以及微动技术应用到地质灾害测量中的可行性。
1 微动探测技术 1.1 微动观测方式微动探测通常采用点测的方式，单个测点常用的布设方式为多重圆观测台阵(见图 1)。为了高效率完成野外施工，多重圆观测台阵之间通常采用等间距布设方式。对于某些特殊的工区，也可采用不等间距方式。圆形台阵的半径可依据探测深度确定。
|
图 1 多重圆观测台阵 Fig. 1 Multicircular observatory array |
SPAC法是微动技术中发展最早应用最广泛的主流方法之一，假定随机振动在时空上具有稳定性，微动波场主要由传播在地表面的频散性波动组成，并且震源是随机分布的。
根据Aki的理论[5]，空间自相关系数可以用零阶贝塞尔函数表示为：
| $ \begin{array}{l} \rho \left( {f, r} \right) = \frac{1}{{2{\rm{ \mathsf{ π} }}}}\int_0^{2{\rm{ \mathsf{ π} }}} {\frac{{S\left( {f, r, \theta } \right)}}{{\sqrt {{S_0}\left( {f, 0} \right){S_r}\left( {f, r} \right)} }}} {\rm{d}}\theta = \\ {J_0}\left( {\frac{{2{\rm{ \mathsf{ π} }}f}}{{c\left( f \right)}}r} \right) = {J_0}\left( {{x_0}} \right)。\end{array} $ | (1) |
式中：S(f, r, θ) 为圆心处观测点与圆周观测点信号的交叉谱；S0(f, 0) 与Sr(f, r) 分别为圆心处和圆周上观测点信号的功率谱；J0为第一类零阶贝塞尔函数；x0=2πfr/c(f) 为J0的宗量；θ为波的入射角度，c(f) 为波的传播速度；ρ(f, r) 为空间自相关系数。
处理数据时，先将实测记录分成若干数据段，剔除干扰明显的记录段，然后将各测点之间进行互相关计算，经过方向平均后求得空间自相关系数ρ(f, r)。由ρ(f, r)=J0(x0) 求出宗量x0，再由x0=2πfr/c(f) 求出相速度c(f)，最终获得相速度频散曲线。
下山单纯形法[15]是一种并行处理的局部优化算法，这种算法首先在N维空间中构造一个包含N+1个顶点的凸多面体，求出各顶点的函数值，并确定其中的最大值、次小值和最小值，然后通过反射、扩张、内缩、缩边等方法求出一个较好解，用之取代最大(差)点，从而构成新的多面体，这样多次迭代则可逼近一个性能较好的极小点。该种方法对目标函数的解析性没有要求，收敛速度快，但是对于目标函数具有多个极小值时，容易陷入一个局部极小值处，全局搜索能力较差。
模拟退火法[16]是一种模拟熔融状态下物体逐渐冷却达到结晶状态的一种反演方法，就是将待反演的模型的参数看作是熔化物体的每一个分子，将目标函数看作是熔化物体的能量函数，通过缓慢减小一个模型温度的控制参数T0来进行迭代反演，并且有一定的概率跳出局部极小值，使目标函数最终达到全局极值点。模拟退火法的主要缺点是收敛速度较慢，为得到一个好的近似最优解，需要反复进行迭代运算，从一个充分大的控制参数初始值T0开始缓慢减小，在变量较多、目标函数较复杂、搜索空间比较大时，完成计算所需时间过于冗长。
本次试验的反演方法结合了下山单纯形法(DHSM)与快速模拟退火法(VFSA)两种方法的优点，联合求解目标函数的极小值[17-18]，该组合反演方法的基本思路为：对于一个给定的初始值，首先用DHSM快速求得一个极小值点，然后改用VFSA随机搜索，跳离该局部最小值，一旦找到一个比该极小值点更小的点时，立即以该点为初始值调用DHSM直接搜索该点附近的另一个极小值点，如此交叉进行，直到满足精度条件，算法结束，得到全局最小值。
2 实验结果 2.1 试验场地概况试验选取莱芜某工程场地进行, 场地的岩性主要以石灰岩为主，溶蚀发育。工作区依据地下水的赋存条件、水力性质等特征将地下水分为松散岩类孔隙水和基岩裂隙水两类。
2.2 数据采集本次微动探测采用重庆奔腾数控研究所生产的WZG-6A工程综合探测仪。探测中主要技术参数如表 1所示。
|
表 1 探测仪主要技术参数 Table 1 Main technical parameters of detector |
根据现场条件，本次微动探测工作共布设两条测线，其中，Ⅰ测线沿山前道路布设，多为硬化路面，局部有基岩出露；Ⅱ测线沿山根布设，多为基岩出露区。每条测线上布设数目不同的测点，观测台阵的布设方式采用1.1中的多重圆观测台阵。
2.3 数据处理采用1.2中的SPAC法反演剪切波速，其基本流程见图 2。
|
图 2 数据处理流程图 Fig. 2 Data processing flow chart |
其中，质量控制主要是通过在数据中引入的一个限制范围来筛选观测曲线与理论曲线间吻合范围。反演方法采用DMSH和VFSA组合反演方法。
2.4 结果分析与讨论本次试验结合微动探测结果及试验区场地概况，对各测线所得探测剖面进行了解释及成因推测。并通过与该试验区高密度电法探测结果对比，对微动探测应用于地质灾害的效果进行了综合推断。
2.4.1 Ⅰ测线该测线共布设11个测点，测点间距5.0 m。图 3为该测线起始点处采集的来源于不同方向随机分布噪声的典型原始信号波形图(0~1.4 s)。由图可见每个通道波形一致性好，形态相似、能量稳定，无持续的强干扰源信号，表明所采集的微动信号在时间域和空间域是稳定的，为后续数据处理过程中获得准确的结果奠定了基础。
|
图 3 Ⅰ测线起始点处原始信号波形图(0~1.4 s) Fig. 3 The original signal waveform figure of lineⅠat initial point |
对该测点求取频散曲线并分别采用三种反演剪切波速，反演时初始模型设置为10层，各层初始参数见如表 2，最大厚度设置为20 m，剪切波速度搜索范围为100~800 m/s，反演误差限为0.002。测点频散曲线与反演层速度见图 4，图中圆圈表示频散曲线；实线表示为组合反演出的剪切波速，共迭代了1 000次；虚线表示为DHSM反演出的剪切波速图。由图可知，组合反演方法反演出三处明显的低速异常，与已知地质资料吻合较好；DHSM反演结果仅在浅层20 m以上与实际地质条件吻合；VFSA反演方法迭代23 000次后得到与组合反演方法基本相同的结果，反演时间过于冗长。
|
|
表 2 初始模型各层参数 Table 2 Initial model parameters of each layer |
|
图 4 Ⅰ测线起始点处频散曲线与反演层速度图 Fig. 4 The dispersion curve and inversed velocity of line Ⅰ at initial point |
对Ⅰ测线上各测点采用与起始处测点相同的反演参数进行反演，将各测点反演速度采用克里金方法插值，得到Ⅰ测线的剪切波速度断面图如图 5所示。结合试验区面波信息和地质资料综合分析：0~10 m为低速层，Vs < 300 m/s，推测为第四系覆盖层或强风化层；10~70 m为中等速度层，形态结构较复杂，平均速度约为480 m/s，15 m处出现明显的低速异常区，推测为溶蚀裂隙发育区，局部发育为溶洞。70~90 m存在小规模的不连续异常区，推测为完整基岩顶部的溶蚀发育。90 m以下Vs>550 m/s，推测为完整岩体。
|
图 5 Ⅰ测线剪切波速断面图 Fig. 5 The shear-wave velocity profile of line Ⅰ |
图 6为Ⅰ测线所进行的高密度电法视电阻率断面图。由图可知，该断面电阻率以高阻为主要特征，反映靠近山前基岩埋深变浅的趋势。高阻异常以不连续的团块异常为主，局部存在低阻凹陷及明显的低阻区，主要发育在0~40 m。其中，15 m处低阻异常与剪切波波速断面图中15 m处低速异常区位置吻合很好。此外，在探测工作结束后施工方根据微动方法确定的异常进行了钻探(见图 5)验证，在15 m位置处钻遇了溶洞，再次证明了微动探测法的准确性。
|
图 6 Ⅰ测线反演视电阻率断面图 Fig. 6 The inversion of apparent resistivity profile on line Ⅰ |
对比两种方法的探测结果，还可发现电阻率剖面反映的异常体范围要比微动探测获得的异常范围要大，这是因为电阻率法是地下不同地质体电性差异的综合反映，容易受到异常体的形态、产状、物性差异等方面的影响，从而产生附加的“假异常”，目前的反演算法尚不足以完全消除这种影响。而微动方法是一种基于波动理论的探测方法，同时考虑了波在介质中传播的时间因素和空间因素，能够更好地反映异常体的空间范围，使得微动探测法对地质灾害定位更加准确、空间分辨率高。
2.4.2 Ⅱ测线该测线共布设了5个观测点，测点间距5.0 m，各测点原始数据具有与Ⅰ测线相同的特征。利用前述相同的处理流程和相同的反演参数，分别对每个测点求取频散曲线并反演了剪切波速。再将各测点反演速度采用克里金插值法做空间插值，得到Ⅱ线的剪切波速断面图(见图 7)。
|
图 7 Ⅱ测线剪切波速断面图 Fig. 7 The shear-wave velocity profile of line Ⅱ |
由图可知：浅部地层主要显示为低速带，推测为第四系及强风化层所致。10 m以下为中高速层，局部出现低速异常区，推测为溶蚀裂隙发育区。在测线位置15 m、深度30~40 m处，存在一处明显的规模较小的低速异常，推测为溶洞。
图 8为Ⅱ测线所进行的高密度电法视电阻率断面图。由图可知，该断面电阻率整体上以高阻为主要特征，局部存在小范围低阻异常，推测为浅部发育的裂隙充水所致。测线0~15 m处有明显低阻异常，该异常区域与微动探测法测线15 m处的波速异常位置吻合很好。
|
图 8 Ⅱ测线反演视电阻率断面图 Fig. 8 The inversion of apparent resistivity profile on line Ⅱ |
本文分析了微动探测技术在地质灾害中的适用性，得出如下结论：
(1) 下山单纯形法(DMSH)和快速模拟退火法(VFSA)组合的反演方法弥补了传统反演方法的不足，具有效率高、方法简便的优点。
(2) 微动探测技术在反映地下空间结构、勘探地质灾害等方面具有很好的可行性和适用性。
随着城市化进程的发展，传统的工程物探方法受到越来越多的制约与限制，微动探测法由于其安全、环保、便捷等特点，在未来的城市及其邻近区域工程物探中具有潜在的应用前景。
| [1] |
林孝城, 刘宏岳, 陈宗刚. 微动探测方法在桩基工程不良地质体探测中的应用--以福建某拟建厂房项目为例[J]. 福建建筑, 2020(8): 90-93+147. Ling X C, Liu H Y, Chen Z G. The Appli-cation of Microtremor detection method in detecting unfavorable Geological body in pile foundation engineering-Take a proposed factory building project in Fujian as an example[J]. Fujian Architecture & Construction, 2020(8): 90-93+147. (  0) |
| [2] |
Rayleigh L. On waves propagated along the plane surface of an elastic solid[J]. Proceeding of the London Mathematic Society, 1887, 17: 4-11.
( 0) |
| [3] |
徐义贤, 张煜. 瑞雷面波在弹性楔状体中的传播特性研究[J]. 工程地球物理学报, 2009, 6(2): 143-149. Xu Y X, Zhang Y. Study of transmission characteristics of Rayleigh waves on an elastic wedge[J]. Chinese Journal of Engineering Geophysics, 2009, 6(2): 143-149. DOI:10.3969/j.issn.1672-7940.2009.02.003 ( 0) |
| [4] |
Haskell N A. The dispersion of surface waves on multilayered media[J]. Bulletin of the Seismological Society of America, 1953, 1(43): 17-34.
( 0) |
| [5] |
Aki K. Space and time spectra of stationary stochastic waves, with special reference to microtremors[J]. Bulletin of the Earthquake Research Institute, 1957(35): 415-456.
( 0) |
| [6] |
Capon J. Applications of detection and estimation theory to large array seismology[J]. Proceedings of the IEEE, 1970, 58(5): 760-770. DOI:10.1109/PROC.1970.7730
( 0) |
| [7] |
Okada H, Matsushima T, Moriya T, et al. An exploration technique using long period microtremors for determination of deep geological structures under urbanized areas[J]. Butsuri-Tansa, 1985, 43: 402-417.
( 0) |
| [8] |
Tavanaei T, Hassanpour J, Memarian H. Urban noises and earthquakes effects on dynamic slope stability-a case study: Arash-Esfandiar tunnel[J]. International Journal of Geo-technical Engineering, 2020, 14(4): 420-427. DOI:10.1080/19386362.2018.1433347
( 0) |
| [9] |
田宝卿, 庞忠和, 孔彦龙, 等. 基于微动技术的浅层地热资源勘探与储量精细评价方法研究[J]. 科技促进发展, 2020, 16(Z1): 367-374. Tian B Q, Pang Z H, Kong Y L, et al. Exploration and amount fine evaluation of geothermal resources based on Microtremor Survey Method[J]. Science & Technology for Development, 2020, 16(Z1): 367-374. ( 0) |
| [10] |
Ku T, Palanidoss S, Zhang Y H, et al. Practical configured microtremor array measurements (MAMs) for the geological investigation of underground space[J]. Underground Space, 2020, 1(4): 1-12.
( 0) |
| [11] |
Holland J H. Adaptation in Natural and Artificial System[M]. Michigan: University of Michigan Press, 1975.
( 0) |
| [12] |
Nelder J, Mead R. A simplex method for function minimization[J]. Comput J, 1965, 7(4): 308-313. DOI:10.1093/comjnl/7.4.308
( 0) |
| [13] |
Kirkpatrick S, Gelatt C D, Vecchi M P. Optimization by simulated annealing[J]. Science, 1983, 220: 671-680. DOI:10.1126/science.220.4598.671
( 0) |
| [14] |
董高玉, 赵晓飞. 地球物理非线性反演方法综述[J]. 科技资讯, 2011(6): 249. Dong G Y, Zhao X F. Overview of Geophysical nonlinear inversion methods[J]. Science & Technology Information, 2011(6): 249. DOI:10.3969/j.issn.1672-3791.2011.06.214 ( 0) |
| [15] |
张志新, 张明廉. 基于并行混沌和单纯形法的混合全局优化算法[J]. 系统仿真学报, 2004(1): 35-37. Zhang Z X, Zhang M L. A hybrid global optimization algorithm based on parallel chaos optimization and simplex search[J]. Journal of System Simulation, 2004(1): 35-37. DOI:10.3969/j.issn.1004-731X.2004.01.012 ( 0) |
| [16] |
师学明, 王家映. 地球物理资料非线性反演方法讲座(三) 模拟退火法[J]. 工程地球物理学报, 2007(3): 165-174. Shi X M, Wang J Y. Simulated annealing method[J]. Chinese Journal of Engineering Geophysics, 2007(3): 165-174. DOI:10.3969/j.issn.1672-7940.2007.03.001 ( 0) |
| [17] |
贾陆锋, 龚俊, 贾善坡. 基于单纯形-模拟退火算法的岩土力学参数反演[J]. 人民长江, 2019, 50(10): 187-192. Jia J F, Gong J, Jia S P. Inversion of geotechnical mechanical parameters based on simplex-simulated annealing algorithm[J]. Yangtze River, 2019, 50(10): 187-192. ( 0) |
| [18] |
邹士新, 杨坤德. 用单纯形模拟退火混合算法反演地声参数[J]. 电声技术, 2006(6): 4-8. Zou S X, Yang K D. Geoacoustic inversion via hybrid simplex simulated annealing algorithm[J]. Elementary Electro-Acoustics, 2006(6): 4-8. ( 0) |
2. Key Laboratory of Marine Environmental Geological Engineering, Shandong Provincial, Qingdao 266100, China