0 引言
根据M. J. Friedel等[1]和Dou Linming等[2]的研究，矿震波速分布与井下应力、强震的分布强相关，因此通过监测目标区域的矿震波速分布可实现矿井灾害预测预报[3]。直接测量井下岩层波速比较困难，通常通过测量与波速相关的物理量，如震动信号的走时[4]或波形本身[5-6]，再通过反演间接获得目标区域的波速分布情况。
现有矿震波速反演大多采用线性反演方法。如基于到时的波速反演，在高频近似假设下，认为震波在均匀介质中沿直线传播，采用线性方程建立波速与到时之间的模型[2，7-8]，进而用SIRT(Simultaneous Iterative Reconstruction,联合代数重建)[9]等方法求解。当震源位置未知时，可通过震源与波速联合反演进行求解[9-11]。但是震波的高频近似假设忽略了波的衍射、衰减、频散等特征，在某些条件下增加反演的多解性。
遗传算法(Genetic Algorithm，GA)[12]、粒子群优化(Particle Swarm Optimization，PSO)算法[19]等非线性反演方法及其衍生算法[13-14]具有天然的并行性，能以较低的代价部署到Spark[15]等分布式集群上，以实现计算速度的提升，是一种更加普适的反演方法。但是非线性反演方法收敛速度较慢、计算量大，特别是处理高维反问题时往往无法在有效时间内获得可用的结果[16]。针对该问题，本文提出一种基于二维快速傅里叶变换(Fast Fourier Transformation，FFT)[17]的矿震波速反演模型空间降维方法。该方法通过二维FFT及高频截断，在损失部分高频信息的前提下降低反演模型空间维数，极大地提高了反演速度。实验结果表明，采用PSO算法对任一50×50目标区域进行反演时，该方法可将反演模型空间维数降低为原来的1/100，并稳定地获得与理论值较为接近的结果。
1 矿震波速反演方法
假设目标区域可能的矿震波速分布模型m=
，vi为第i(i=1,2,…,n,n为网格数)个网格中的波速，
和
为每个网格中波速的上下限。可能的可观测物理量构成观测值d，m到d的映射为G，d到m的映射为G-1，则矿震波速反演的数学模型[8]为
m=G-1(d)
(1)
波速正演的数学模型为
d=G(m)
(2)
显然，若G-1已知，通过直接求解式(1)即可获得m。但由于实际情况复杂，直接获得G-1较为困难，而G可通过理论推导或实验测量等方式较容易地获得，因此，井下波速反演问题一般可等价于以下有约束最优化问题。
已知G和实际观测值dobs，求最优矿震波速分布模型
为第i个网格中的最优波速，使得目标函数f(d)最小。f(d)可选为
(3)
实际情况下G通过射线追踪建立,因此d=G(m)可由以下线性方程组表示：
d=Gm
(4)
以基于射线追踪的二维波速反演为例说明矿震波速反演原理，如图1所示。
图1 矿震波速反演原理
Fig.1 Mine seismic wave velocity inversion principle
反演步骤如下。
(1) 对目标区域进行网格剖分，分成10×30个网格，得到待解模型m。假设每个网格中波速恒定，则待解的300个波速构成整个模型空间，此时模型空间维数为300。
(2) 确定m与观测值d之间的映射关系G，建立正演数学模型d=G(m)。矿震波速反演通常选择震波的初至到时作为观测值，而m和d之间的映射关系通常选用射线追踪[3]来建立，以降低正演数学模型计算的复杂度。
(3) 选择最优化算法求解m。已知实际观测值dobs及m到d的映射关系G，求最优模型mopt，使得目标函数
最小。
2 矿震波速反演模型空间降维方法
2.1 模型空间降维原理
波速反演的目的是获得目标区域各网格的波速，通过波速不同区分目标区域构造或应力分布情况。显然，在保证一定的网格剖分密度情况下，有意义的波速分布图像与无意义的图像有明显区别，如图2所示,其中图像划分为50×50网格,(x,y)为网格坐标。
图2 波速分布图像
Fig.2 Wave velocity distribution image
从图2可得出以下结论。
(1) 有意义的波速分布图像中，各网格的波速不是完全随机的，相邻两格的波速大概率一致或接近。为方便讨论，限定图2(a)和图2(b)中每个网格波速为0～255 m/s且为整数。
(2) 有意义图像的主要能量集中在低频部分，无意义图像的能量在整个频段上近似均匀分布。可见，对于有意义图像，仅保留其低频部分(高频截断)可近似得到原始图像，如图3所示。
图2中幅度谱的中心点坐标为(26，26)且频谱应关于中心点对称[17]。假设频谱中保留点的横纵坐标到中心点横纵坐标的距离为D，则描述一幅图像所需的变量个数为
Nvar=2[(2D+1)D+D+1]-1=
4D2+4D+1
(5)
由式(5)可知，在保证图像近似为原始图像的条 件下，描述一幅图像所需的变量个数由原始图像的50×50=2 500个极大地降低了，如图3(f)只需361个变量即可。
图3 原始图像与高频截断图像
Fig.3 Original image and images with high frequency truncation
2.2 基于模型空间降维的波速反演
在损失部分高频信息的情况下，描述一幅图像所需的变量个数可大幅减少，意味着波速反演模型空间维数可大幅降低，从而保证基于非线性最优化方法的波速反演在有效时间内完成，具体步骤如下。
(1) 确定目标区域的网格剖分及每格波速上下限。假设目标区域被剖分为M行N列，每格波速范围为
在保证一定的波速反演精度条件下，需对波速进行离散化及归一化处理。
(2) 确定拟保留的变量个数及每个变量范围。由式(5)得变量个数Nvar，其中(Nvar+1)/2个为模(下限为0，上限可通过实验确定)，(Nvar+1)/2-1个为幅角(范围为[-π,π])。将目标区域图像的每一格波速均置为
对图像进行二维FFT并取模，则中心点高度可作为模的上限。若对波速进行归一化处理，则模的上限为一固定值，其值取决于归一化的波速上限。
通过步骤(2)，正演数学模型由d=G(m)变为d=G(G′(m))，其中G′为理想低通滤波，此时目标区域的描述由MN维向量m变为Nvar维向量G′(m)，即反演模型空间由空域变为频域。
(3) 构造目标函数，选择最优化方法进行问题求解。目标函数的构造方法不唯一，本文选择
作为寻优目标。考虑到目标函数维数高且峰值较多，选择PSO算法作为最优化方法。
3 实验与分析
3.1 有效性实验与分析
采用仿真实验对矿震波速反演模型空间降维方法的有效性进行验证。
(1) 任意给出某种波速分布mtheo及震源、检波器布置，将mtheo作为真实波速分布，通过正演数学模型获得各检波器理论到时dtheo。为保证计算速度，选择射线追踪建立波速分布模型到观测值之间的映射关系。
(2) 令实际观测值dobs=dtheo，即认为观测结果不含误差，相当于式(2)中d，G已知，m待求。采用2.2节方法求得最优模型mopt。
(3) 对比mtheo与mopt，评价反演结果。
选择图3(a)作为真实波速分布，实际波速范围为3 000～6 000 m/s，通过离散化及线性映射可将波速范围映射为0～255 m/s之间的整数。反演参数见表1，真实波速分布与震源、检波器布置如图4所示。波速反演结果如图5所示。本文保留频域坐标(24，24)～(28，28)的图像，模型空间维数可降为原来的1/100。
表1 反演参数
Table 1 Inversion parameters
图4 真实波速分布与震源、检波器布置
Fig.4 True wave velocity distribution and location of seismic source and seismometer
图5 有效性实验波速反演结果
Fig.5 Wave velocity inversion results in validity experiment
以
作为反演结果与目标接近程度的评价指标，其中
为第i行第j列网格的最优波速，
为第i行第j列网格的真实波速。设中间结果
评价指标为
(6)
选择图5(b)与图5(a)对比(b→a)、图5(c)—图5(f)与图5(a)对比(c—f→a)、图5(c)—图5(f)与图5(b)对比(c—f→b)，所得mcomp_mean如图6所示。
图6 不同迭代次数下反演结果评价指标
Fig.6 Evaluating indicator of inversion results under different iterations
由图5、图6可得出以下结论。
(1) 在PSO算法迭代过程中，mcomp_mean随迭代次数的增加逐渐减小，即反演结果与目标在形状与颜色上越来越接近，说明使用PSO算法进行基于到时的波速反演是有效的。
(2) 反演结果与目标之间的接近程度均大于图5(b)与图5(a)之间的接近程度，说明反演结果与理论模型间仍有较大差距。由于实验选择dobs=dtheo，可排除测量误差对反演结果的影响，又由于正反演模型一致(均由射线追踪建立)，可排除模型误差带来的影响，所以造成反演结果误差较大的主要原因在于观测资料不完整。
为验证该结论，将图4中纵向的震源和检波器取消(取消的震源坐标为(0，100)，(0，200)，(0，300)，(0，400)，(0，500)，取消的检波器坐标为(2 000，0)，(2 000，50)，(2 000，100)，(2 000，150)，(2 000，200)，(2 000，250)，(2 000，300)，(2 000，350)，(2 000，400)，(2 000，450))，再用相同的方法与参数进行反演实验，结果如图7所示。
从图7可看出，由于观测资料不足，反演结果与理论模型间从直观上已无相似性。基于式(6)计算图7与图5(a)、图5(b)的mcomp_mean，分别为5 121.2，4 842.8，远大于图5(f)与图5(a)、图5(f)与图5(b)间的mcomp_mean，说明此时获得的反演结果与真实值误差较大，这与直观感觉一致。
(3) 对比图5(a)与图5(c)—图5(f)可知，反演结果的图形分界线均为曲线，这是由于反演步骤(2)中使用二维FFT对待反演的模型空间进行展开且 做了高频截断。由于模型空间高频成分缺失，所以反演获得的图像边缘呈现正弦状波动，如图5(f)的1与2、2与3交界处，这种曲线分界线会导致mcomp_mean增大。通过增加反演步骤(2)中保留的频谱范围，可以增强模型对高频变化的表现能力，但会导致待反演变量急剧增加，从而增加反演计算难度，使得计算难以收敛。
图7 取消纵向震源和检波器后波速反演结果
Fig.7 Wave velocity inversion results after canceling vertical seismic source and seismometer
3.2 稳定性实验与分析
由于PSO算法为蒙特卡洛类方法，且从理论上并未证明其收敛性，通过该方法获得的结果可能具有随机性，且考虑到波速反演本身的多解性，有必要考查矿震波速反演模型空间降维方法的稳定性。使用3.1节中的反演参数进行100次迭代，共进行3次反演实验，结果如图8所示。
从图8可看出，3次反演的结果各不相同，体现了蒙特卡洛类方法的随机性及波速反演多解性对最终反演结果的影响。从直观上看，图8(c)—图8(e)形状各异，但仍基本反映了图8(a)的特征，即波速分布由上到下逐级降低，其中图8(c)与真实波速分布最接近，图8(e)次之，图8(d)最差。图8(c)—图8(e)对应的目标函数值
分别为0.816，1.168，0.991，与直观感觉一致。为衡量反演结果与真实值的接近程度，分别计算3次反演结果与图8(a)、图8(b)的mcomp_mean及2组mcomp_mean均值，结果如图9所示。
根据图9可看出反演结果与真实值间的接近程度与图9中直观感觉一致。若以图8(c)与图8(a)的mcomp_mean为基准，图8(d)与图8(a)、图8(e)与图8(a)的mcomp_mean及mcomp_mean均值分别比其大43.3%，24.5%，22.6%；若以图8(c)与图8(b)的mcomp_mean为基准，图8(d)与图8(b)、图8(e)与图8(b)的mcomp_mean及mcomp_mean均值分别比其大54.6%，31.4%，28.7%。图7与图8(a)、图7与图8(b)的mcomp_mean分别比图8(c)与图8(a)、图8(c)
(e) 第3次反演结果
图8 稳定性实验波速反演结果
Fig.8 Wave velocity inversion results in stability experiment
与图8(b)的mcomp_mean大109.6%，145%。可见采用本文方法进行波速反演所得结果具有随机性，但与真实值偏离程度不高，仍能给出目标区域的波速分布特征。
图9 反演结果评价指标
Fig.9 Evaluating indicator of inversion results
4 结语
针对矿震波速反演模型空间过大的问题，从正演模型出发，通过二维FFT及高频截断，在损失部分高频信息的前提下将待寻优的模型空间由空域变为频域，极大地降低了模型空间维数(为原来的1/100)，并基于PSO算法实现了无预置初值的波速非线性反演。实验结果表明，该矿震波速反演模型空间降维方法在进行非线性寻优时能够较快地收敛，在观测数据足够的条件下可稳定地获得与理论值较为接近的结果。
致谢：感谢中国矿业大学物联网(感知矿山)研究中心丁恩杰教授在本文研究方向和研究方法上的指导！
参考文献( References) :
[1] FRIEDEL M J,JACKSON M J,SCOTT D F,et al.3-D tomographic imaging of anomalous conditions in a deep silver mine[J].Journal of Applied Geophysics,1995,34(1):1-21.
[2] DOU Linming,CHEN Tongjun,GONG Siyuan,et al.Rockburst hazard determination by using computed tomography technology in deep workface[J].Safety Science,2012,50(4):736-740.
[3] 蔡武,窦林名,李振雷,等.矿震震动波速度层析成像评估冲击危险的验证[J].地球物理学报,2016,59(1):252-262.
CAI Wu,DOU Linming,LI Zhenlei,et al.Verification of passive seismic velocity tomography in rock burst hazard assessment[J].Chinese Journal of Geophysics,2016,59(1):252-262.
[4] AKI K,LEE W.Determination of three-dimensional velocity anomalies under a seismic array using first P arrival times from local earthquakes, 1, a homogeneous initial model[J].Journal of Geophysical Research,1976,81(23):4381-4399.
[5] GAUTHIER O,VIRIEUX J,TARANTOLA A.Two-dimensional nonlinear inversion of seismic waveforms: numerical results[J].Geophysics,1986,51(7):1387-1403.
[6] MORA P.Nonlinear two-dimensional elastic inversion of multioffset seismic data[J].Geophysics,1987,52(52):1211-1228.
[7] 窦林名,蔡武,巩思园,等.冲击危险性动态预测的震动波CT技术研究[J].煤炭学报,2014,39(2):238-244.
DOU Linming,CAI Wu,GONG Siyuan,et al.Dynamic risk assessment of rock burst based on the technology of seismic computed tomography detection[J].Journal of China Coal Society,2014,39(2):238-244.
[8] ASTER R C, BORCHERS B, THURBER C H. Parameter estimation and inverse problems[M].Salt Lake City:Academic Press,2013:411-421.
CHEN Li, DU Wenhua, ZENG Zhiqiang, et al. Automatic separation method of coal and gangue based on wavelet transform[J].Industry and Mine Automation，2018,44(12)：60-64.
[9] GILBERT P.Iterative methods for the three-dimensional reconstruction of an object from projections[J].Journal of Theoretical Biology,1972,36(1):105-117.
[10] 刘福田.震源位置和速度结构的联合反演(Ⅰ)——理论和方法[J].地球物理学报,1984,27(2):167-175.
LIU Futian.Simultaneous inversion of earthquake hypocenters and velocity structure(I)-theory and method[J].Chinese Journal of Geophysics,1984,27(2):167-175.
[11] 黄国娇,白超英.多震相走时联合三参数同时反演成像[J].地球物理学报,2013,56(12):4215-4225.
HUANG Guojiao,BAI Chaoying.Simultaneous inversion of three model parameters with multiple classes of arrival times[J].Chinese Journal of Geophysics,2013,56(12):4215-4225.
[12] GOLDBERG D E, HOLLAND J H. Genetic algorithms and machine learning[J]. Machine Learning,1988,3(2):95-99.
[13] CLAUDE S,GEOFFREY W.Encyclopedia of machine learning[M].New York:Springer,2011:
760-766.
[14] 师黎静,陶夏新,赵纪生,等.一种波速结构的两步优化反演策略[J].地球物理学报,2009,52(8):2105-2112.
SHI Lijing,TAO Xiaxin,ZHAO Jisheng,et al.A two-step optimization strategy for S-wave velocity structure inversion[J].Chinese Journal of Geophysics,2009,52(8):2105-2112.
[15] ZAHARIA M,CHOWDHURY M,FRANKLIN M J,et al.Spark:cluster computing with working sets[C]//Usenix Conference on Hot Topics in Cloud Computing,Berkeley,2010:10.
[16] 王家映.地球物理资料非线性反演方法讲座(一) 地球物理反演问题概述[J].工程地球物理学报,2007,4(1):1-3.
WANG Jiaying.Lecture on non-linear inverse methods in geophysics(I) introduction to geophysical inverse problems[J].Chinese Journal of Engineering Geophysics,2007,4(1):1-3.
[17] 郑君里,应启珩,杨为理.信号与系统[M].2版.北京:高等教育出版社,2000.