爆炸近场高精度P波角度提取方法

孟 铭，杨 明，李 剑，韩 焱

(中北大学信息探测与处理山西省重点实验室，山西 太原 030051)

0 引言

基于偏振角度信息的DOA(Direct of Angle)定位方法是实现地下浅层震源近场定位的主要手段[1]，该方法采用三分量数据来计算P波的偏振方向，通过多个检波器的P波偏振方向求取交汇点，从而实现震源定位。由于该方法无需求取初至时间，原则上仅需要两组检波器数据即可实现震源定位，因此是地下浅层空间定位领域研究的热点。然而在爆炸点近场定位过程中，由于群波混叠严重、多径干扰效应大，导致P波、S波偏振信息混叠，P波偏振角度 “提不准”，最终造成DOA定位模型“建不准”。

在直达P波角度提取算法的研究中，目前国内外学者主要采用MUSIC(Multiple Signal Classification)算法、ESPRIT(Estimating Signal Parameter via Rotational Invariance Techniques)算法[2-3]、协方差矩阵重构法和极化分析法等[4]。其中MUSIC算法、ESPRIT算法和协方差矩阵重构法均是基于均匀线性阵列构型的远场角度估计方法，极化分析法最重要的工作就是时窗长度的选取，在信号强混叠的情况下，极化分析法难以实现时窗长度的准确选取[5-6]。

因此上述方法只适用于远场大区域、大深度条件下P波偏振角度的提取，在近场信号强混叠情况下无法适用。针对这个问题，本文利用地下波场在Radon域的聚焦差异和P波的偏振特性，开展了基于高分辨率抛物Radon(High Resolution Parabolic Radon Transform，HRP-Radon)联合ACM的直达P波偏振角度提取方法的研究

1 直达P波偏振角度提取原理

由地震波理论可知,同一类型的波引起空间质点的运动,其轨迹近似于椭球体，如图1所示[7]。按照概率统计学方法，通过统计直达P波一段时间内主运动事件，即可得到该波主运动方向。

图1 笛卡尔坐标系下的空间质点运动椭球示意图Fig.1 Schematic diagram of spatial particle motion ellipsoid in Cartesian coordinate system

目前主要采用SCM(标准协方差极化滤波)提取直达P波角度信息。设直达P波对应三分量震动数据ui(t),(i=x,y,z)，t∈[t1～t2]。求取各分量数据的平均值μi(ξ)：

(1)

(2)

(3)

(4)

(5)

式中，θ0(t)为瞬时方位角，δ0(t)为瞬时倾角。基于SCM原理的偏振角度提取算法简单，但是对时窗长度的选取较为敏感，并且在给定长度的时窗内求得的极化参数不具有时变特性。同时由于爆炸近场存在波形混叠严重、多径干扰效应大等问题，导致无法有效和准确的提取时变信号的偏振角度信息。

2 爆炸近场直达P波偏振角度提取方法的改进

针对爆炸近场群波混叠条件下，SCM算法无法提取有效直达P波角度信息的难题，本文利用传感器阵列的波场信息，通过HRP-Radon变换和远程自然分离的P波振相特征在Radon域提取混叠信号中的P波信号。利用震动波的时变特性，自适应调整SCM算法中时窗长度，提取直达P波的瞬时极化特性，最终获取精准的偏振角度信息。

2.1 基于HRP-Radon直达P波信号提取

本文提出的HRP-Radon变换结合P波波形振相特征的方法。HRP-Radon变换可根据横纵波的特征参数(慢度、曲率等)差异将P波和S波映射为Radon域内不同位置的峰值点，实现纵、横波的分离[9-11]。该方法首先提取P波的首振相、主振相和尾振相曲线，通过HRP-Radon变换得到这些振相曲线对应的Radon域投影，并将其作为Radon域投影区域P波提取边界条件。然后参照信号特征边界条件保留混叠波形中P波在Radon域有效区域，并根据P波有效投影区域实现P波反演重建。P波分离流程如图2所示。

图2 P波分离流程图Fig.2 P-wave separation flow chart

首先对三分量混叠信号从时间域数据沿抛物线路径进行积分变换到Radon域，并完成其逆变换(从Radon域变换到时间域)，具体公式如下：

(6)

(7)

对式(6)和式(7)两端进行傅里叶变换，然后将公式转变为矩阵形式：

m=LHd

(8)

d=Lm

(9)

采用高分辨率抛物线Radon变换，可以较大程度地降低由于空间采样范围有限导致的截断效应对Radon域内信号分辨率的影响[12-13]。式(8)中的函数可优化为:

(10)

式(10)中的加权矩阵Qm与m的关系表达式如下：

(11)

共轭梯度算法占用内存小而且计算效率高，在式(10)中，使用共轭梯度算法完成矩阵(LHL+Qm(mk))的求逆过程。高分辨率抛物线Radon变换的数据重建误差小，在变换过程中能量损耗特别低，变换后数据可以得到较完整的恢复，有效的实现P波的分离和反演重建。

在实现强混叠条件下直达P波有效分离的基础上，还需对分离P波进行ACM极化分析，以提取P波偏振角度等极化参数，为构建高精度DOA定位模型提供精确的方向参量信息。

2.2 基于ACM的直达P波极化角度的提取

本文在SCM算法的基础上，利用爆炸波时变特性瞬时极化，使时窗长度自适应调整到时窗内信号的最小周期，构建了ACM(自适应协方差极化滤波)偏振角度提取模型。

(12)

(13)

根据地震波理论，极化度T是评价质点在地下空间运动偏振特性的一个指标，如图1所示。当极化度接近于1时，该质点在空间中呈线性偏振，由主特征值对应的特征向量得到的偏振角度信息可以最大程度反应该质点在地下空间的运动方向。当极化度接近于0时，该质点在空间中呈无序圆形运动，得到的偏振角度信息无法有效的衡量该质点在地下空间的运动方向。因此可以用极化度T来评价P波在地下空间中的偏振特性，即偏振角度信息的获取精度。具体公式如下所示：

(14)

式(14)中，ρ表示主椭圆率极化参量，ρ1表示次椭圆率极化参量。

3 算法仿真

地下浅层爆炸近场的仿真模型如下：假设在地下介质为均匀水平层状介质。建立XYZ三轴坐标系，地平面为XY平面，炸点位置为(0，30，-200)。以原点为中心，在X轴上布置21个传感器，间距为31 m，P波传播速度为5 500 m/s，S波传播速度为3 800 m/s。采用雷克子波激励方式产生相应信号，设信号采样率为20 kHz，采样时间为0.16 s。爆破震动波的初至时刻为0.01 s，频率为160 Hz，衰减因子为e-i/160(i=0,1,…,n)，横波在0.02 s产生，频率为80 Hz，衰减因子为e-i/160(i=0,1,…,n)，噪声源为高斯白噪声，信噪比为35 dB,并抽取第5道数据为例进行算法仿真验证。

图3 仿真模型示意图Fig.3 Schematic diagram of simulation model

通过矢量计算生成传感器阵列的三轴数据，具体P波偏振角度信息提取方法如图4所示。

图4 混叠波形及Radon域投影Fig.4 Aliasing waveform and Radon domain projection

图4(a)是模拟三分量传感器接收到的地下爆炸产生的横纵波混叠信号，图4(b)为混叠信号经过HRP-Radon变换后在Radon域的投影图。

图5(a)是P波振相特征曲线，图5(b)为P波振相位特征曲线在Radon域的投影区域，可作为P波Radon域有效区域提取的依据。

由图4(b)和图5(b)结合可以得出P波合理保留区域，如图6(a)所示。图6(b)为反变换重建后的P波。

图5 P波振相特征曲线及Radon域投影Fig.5 P-wave vibration phase characteristic curve and Radon domain projection

图7、图8、图9分别为混叠波形和分离P波在x，y，z轴的瞬时夹角对比图。其中信号有效区域由信号初至波到时时刻开始到设定时间0.08 s结束，信号无效区域表示激励产生的高斯白噪声。实际数据处理时以信号初至波到时作为有效区域起始位置,以极化度瞬时变化时刻作为有效区域结束时间的划分节点。通过图10的极化度曲线对比图可以看出分离P波偏振角度的极化度曲线[0.95，1]之间，并由表1可得，HRP-Radon变换联合ACM算法能够得到高精度直达P波的角度信息。

图6 P波Radon域分离及时域重建Fig.6 P-wave Radon domain separation and time domain reconstruction

图7 混叠波形与分离P波x轴夹角Fig.7 Inclusion angle between aliasing waveform and separated P-wave x-axis

图8 混叠波形与分离P波y轴夹角Fig.8 Inclusion angle between aliasing waveform and separated P-wave y-axis

图9 混叠波形与分离P波z轴夹角Fig.9 Inclusion angle between aliasing waveform and separated P-wave z-axis

图10 混叠波形与分离P波极化度曲线Fig.10 Mixed Waveform and Separated P-Wave Polarity Curve

表1P波的偏振信息表
Tab.1 The polarization information table of P wave

理论值实际值绝对误差提取精度vx/rad1.02081.0210.00020.9998vy/rad1.69761.6980.00040.9997vz/rad0.56780.56790.00010.9998

4 结论

本文提出了爆炸近场高精度P波角度提取方法。该方法首先利用P波、S波在Radon聚焦特性的不同，以HRP-Radon分离混叠波群中的P波信号。其次，利用爆炸近场震动信号的时变特性，自适应调整波束角度信息检测的窗长，在SCM算法基础上进行改进，构建基于ACM的偏振角度提取模型，以实现有效P波角度的提取。最后，利用极化度对波束角度信息的提取精度进行评价。数值仿真结果表明,上述方法能够有效获取爆炸近场混叠信号中P波角度信息，本文研究成果将为爆炸近场波场特征参数高精度提取提供有力的技术支撑，对地下空间近场震源定位研究具有一定的参考价值。但由于本算法利用远场P波振相特征作为有效区域划分边界条件，因此在保证直达P波角度提取精度的前提下，如何优化近、远场传感器阵列布设是下一步研究的重点。