一种面向预警的地震信号压缩感知采集传输方法
技术领域
本发明涉及地震信号检测与预警技术领域，尤其涉及一种面向预警的地震信号压缩感知采集传输方法。
背景技术
地震时严重威胁人类生命财产安全的自然灾害之一，常常造成严重的建筑物破坏与人员伤亡，并引发火灾、水灾、管道泄漏、细菌和放射性物质扩散，以及山体滑坡、海啸、崩塌以及地裂缝等次生灾害。大震发生后在很短的时间内释放出大量的能量，几十秒甚至是几秒之内就造成毁灭性的破坏。随着全球都市化的发展，城市逐渐成为政治、经济、文化、社会活动的中心，并且城市人口密度大、建筑物多，突发性的地震给城市造成了严重的威胁。地震的发生与否不受人类的控制，这是地震和许多其它自然灾害的共同特征。由于地震孕育和发生的过程极度复杂，目前不能获得可靠的短期地震预报。世界各国已经把地震减灾措施的重点放在了提高建筑抗震能力、地震预警和地震应急救援等方面。
在地震发生后，如果在断层破裂初期就能利用P波识别地震的最终规模，那么在破坏性的S波到达前，将地震预警信息通过电信号发送到即将遭受地震袭击的预警目标区。这样接收到预警信息的地区就可以预先采取减灾措施。地震预警可在地震发生后，破坏性地震波到达前，提供数秒至数十秒的预警时间，并采取相应应急处置措施，是近年来新发展起来的减轻地震损失、降低地震次生灾害、减少人员伤亡的有效手段。目前在日本、美国、墨西哥、土耳其、中国等都建立了较成熟的地震预警系统。
地震预警系统其主要流程是通过数据采集、数据传输、数据处理及分析判断等步骤实现，数据采集环节观测到的地震数据输送到数据处理环节进行加工，削弱噪声的干扰，提高地震信号的信噪比和分辨率，降低地震信号失真度，以便得出较为精确的预警结果。由于地震预警系统观测台网的数量不断增加，采集的数据呈几何级数增长，而地震预警对于数据的时效性要求较高，采集的数据与数据传输速度之间存在着矛盾。常见的地震台网数据采集基于经典的奈奎斯特采样定理；从数据处理流程上来看，现有的地震观测方法是将实时的采样数据经过完整采集记录后再进行传输，采集的数据中存在冗余度，采集过程会消耗大量系统资源，而且对采集的数据进行压缩和传输，更是消耗了相当的系统资源，为数据的记录和实时传输带来了较大压力。如在采样的同时进行信号的压缩，则可以消耗较低的系统资源实现信号的采集、传输和重构。
发明内容
有鉴于此，本发明提出了一种在采集地震信号的同时对地震信号进行压缩，降低采集的数据中的冗余度，并可靠还原该地震信号的面向预警的地震信号压缩感知采集传输方法。
本发明的技术方案是这样实现的：本发明提供了一种面向预警的地震信号压缩感知采集传输方法，包括如下步骤：
S1：地震预警系统对采集的信号中含有噪音的原始强震信号进行分帧处理，在兼顾计算效率的同时还提升信号恢复的效果，相邻的帧之间无重叠，得到分帧处理后的强震信号x；
S2：对分帧处理后的强震信号x进行压缩采样，构建基于压缩感知的强震信号x的模型：令Ψ为稀疏基；Ψ'是稀疏基Ψ的转置矩阵；x是含噪音的分帧原始强震信号，x不是稀疏的；a是x在稀疏基Ψ的稀疏表示向量，a是稀疏的，a＝Ψ'x；Φ为与稀疏基Ψ不相关的观测矩阵，则有：
y＝Φa＝ΦΨ'x；y为x在观测矩阵Φ下的低维观测向量；
S3：完成压缩采样后得到低维观测向量y后，低维观测向量y送入地震预警系统的通信发送缓冲中，并被在接收端存入接收缓冲；每隔K帧，地震预警系统的发送缓冲还将一帧分帧处理后的原始强震信号的发送到系统的接收缓冲中；
S4：信号接收压缩采样后的强震信号和原始强震信号，并对压缩采样后的强震信号进行重构；重构算法为：
根据观测矩阵Φ、Ψ'是稀疏基Ψ的转置矩阵、含噪音的分帧原始强震信号x和压缩后的观测向量y；将ΦΨ'作为恢复矩阵，采用正交匹配追踪算法对线性采样的压缩后的M维的观测向量y进行重构，得到分帧处理后的强震信号x的逼近信号
S5：将重构后的强震信号x的逼近信号
S6：每隔K帧，对重构后的强震信号的逼近信号
在以上技术方案的基础上，优选的，所述分帧处理是设定每一分帧内的数据点的数量，每一分帧内的数据点的数量为200、400、600、512或者800中的一种，设定每秒采样率为100sps或者200sps。
进一步优选的，所述稀疏基Ψ为DCT、DFT或者小波DWT基。
更进一步优选的，所述观测矩阵Φ为随机高斯矩阵、随机贝努利矩阵、稀疏随机测量矩阵、部分傅立叶矩阵或者部分哈达玛矩阵。
再进一步优选的，所述对重构后的强震信号
其中，N为分帧长度，表示每一分帧中强震数据的个数；测量数M，即压缩采样后的强震信号的观测数，xi为地震预警系统发送的一帧被分帧处理后的强震信号；
在以上技术方案的基础上，优选的，还包括信号捡拾步骤，信号捡拾的方法为：地震预警系统在分帧处理前对采集到的地震P波信号进行0.1—33Hz的带通滤波，采用STA/LTA方法对带通滤波后的地震P波进行粗略捡拾：
STA为短时间窗，LTA为长时间窗；在地震信号P波到达时STA/LTA会有突变，当两者的比值大于设定阈值R时，判断为有效地震信号：
CFP为地震信号P波的捡拾特征函数，CFP＝Y(k)2+[Y(k)-Y(k-1)]2；Y(k)和Y(k-1)分别为k和k-1时刻垂直方向的加速度记录值；i为当前时刻；k1和k2是当前时刻之间某时刻，且k2＜k1＜i；上式中，公式分子部分STA短时间窗的数据点数为50个，公式分子部分LTA长时间窗的数据点数为1000个数据；
应用STA/LTA方法只能粗略识别地震P波震相的到时，为更加精确对地震P波震相的到时进行捡拾，采用VAR—AIC法进行精确捡拾：
AIC(n)＝n·log{var(x[1,n])}+(L-n-1)·log{var(x[n+1,L])}；
n的范围为地震P波STA或者LTA窗口内所有的采样点；var表示方差；窗内AIC的最小值是到时点，然后将采用STA/LTA方法粗略捡拾的P波数据根据采集时间向前推AIC的最小值相应的点数或者相应的采样时间，实现信号精确捡拾。
本发明提供的一种面向预警的地震信号压缩感知采集传输方法，相对于现有技术，具有以下有益效果：
(1)本发明通过设定稀疏基、观测矩阵、压缩采样模型和重构算法，实现对不稀疏信号的稀疏表示，降低采样数据量的冗余度，保证原始信号重构的可靠性和再现性，节省了系统开销，提高了信号处理效率；
(2)信号预先进行分帧处理，一方面可以减小数据处理量，另一方面后续的评价过程能够及时调整后续分帧信号采用的稀疏基和观测矩阵，及时提高信号压缩和重构的精度；
(3)定期对原始强震信号和重构后的强震信号进行信噪比和均方误差的评价，能及时选择最适合的稀疏基和观测矩阵进行信号压缩和重构，提高实时采样的信号的可靠性。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明一种面向预警的地震信号压缩感知采集传输方法的流程图。
具体实施方式
下面将结合本发明实施方式，对本发明实施方式中的技术方案进行清楚、完整地描述，显然，所描述的实施方式仅仅是本发明一部分实施方式，而不是全部的实施方式。基于本发明中的实施方式，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施方式，都属于本发明保护的范围。
如图1所示，本发明提供了一种面向预警的地震信号压缩感知采集传输方法，具体包括如下步骤：
S1：地震预警系统对采集的信号中含有噪音的原始强震信号进行分帧处理，在兼顾计算效率的同时还提升信号恢复的效果，相邻的帧之间无重叠，得到分帧处理后的强震信号x；分帧处理是设定每一分帧内的数据点的数量，每一分帧内的数据点的数量为200、400、600、512或者800中的一种，设定每秒采样率为100sps或者200sps，即每秒采样数量，如采样率为100sps，分帧点数为400，则采样时间为4秒，相邻点数之间的时间间隔为0.01秒；相邻的帧不重叠，可以是连续分帧的，也可以是间隔分帧的；
S2：对分帧处理后的强震信号x进行压缩采样，构建基于压缩感知的强震信号x的模型：令Ψ为稀疏基；Ψ'是稀疏基Ψ的转置矩阵；x是含噪音的分帧原始强震信号，x不是稀疏的；a是x在稀疏基Ψ的稀疏表示向量，a是稀疏的，a＝Ψ'x；Φ为与稀疏基Ψ不相关的观测矩阵，则有：
y＝Φa＝ΦΨ'x；y为x在观测矩阵Φ下的低维观测向量；
其中稀疏基Ψ可以是DCT、DFT或者小波DWT基：
DFT为离散傅里叶变换，它把时域和频域联系起来，将随时间变换的信号转换为频域内的傅里叶变换；
DCT为离散余弦变换，类似于离散傅里叶变换，但是只使用实数。离散余弦变换后在频域中得到的也是实数信号，相比DFT可以减少一半的计算量；DCT具有很强的能量集中作用，其能量主要集中在低频分量上；
小波DWT基克服了传统傅里叶变换无法表述信号局部特征的不足，较好的解决时间和频率分辨率的矛盾，在时—频平面的不同位置具有不同的分辨率；
观测矩阵需要足够稳定，为了确保信号的线性观测能够保持信号的原始结构，观测矩阵必须满足约束等距性条件(RIP)，本发明的观测矩阵Φ为满足约束等距性条件RIP的随机高斯矩阵、随机贝努利矩阵、稀疏随机测量矩阵、部分傅立叶矩阵或者部分哈达玛矩阵：
随机高斯矩阵中的元素在(0，1)区间服从随机独立的正态分布，从而使其与大多数正交基(稀疏基)不相关，满足约束等距条件；
随机贝努利矩阵的每个元素独立地服从对称的贝努利分布，相对于随机高斯测量矩阵，由于随机贝努力测量矩阵的元素为±1，所以在实际应用中更容易实现和存储；
稀疏随机测量矩阵的构造方法如下：先生成一个大小为M×N阶的全零矩阵Φ，且M＜N，然后对于矩阵Φ中的每一列，随机选取K个位置设置为1，K＜M；稀疏随机测量矩阵的每一列只有K个非零元素，结构简单，实际应用中极易构造和保存；
部分傅里叶矩阵属于部分正交矩阵，满足约束等距条件(RIP)，可以作为观测矩阵；其构造方法是，首先生成大小为N×N的傅里叶方阵，然后在生成矩阵中随机地选取M行向量，构成一个M×N的矩阵，最后对列进行归一化；与部分傅里叶矩阵类似，部分哈达玛矩阵亦属于部分正交矩阵，具有较强的非相关性和部分正交性，所以与其他确定性测量矩阵相比，该测量矩阵精确重构所需要的测量次数较少，重建效果好。但是其维数N必须是2的整数倍，一定程度上限制了该矩阵的应用范围及场合；
本发明优选的观测矩阵Φ为随机高斯矩阵、随机贝努利矩阵或者稀疏随机测量矩阵；
S3：完成压缩采样后得到低维观测向量y后，低维观测向量y送入地震预警系统的发送缓冲中并被发送到系统的接收缓冲；每隔K帧，地震预警系统的发送缓冲还将一帧被分帧处理后的强震信号的发送到系统的接收缓冲中；
S4：信号接收压缩采样后的强震信号和原始强震信号，并对压缩采样后的强震信号进行重构；重构算法为：
根据观测矩阵Φ、Ψ'是稀疏基Ψ的转置矩阵、含噪音的分帧原始强震信号x和压缩后的观测向量y；将ΦΨ'作为恢复矩阵，采用正交匹配追踪算法对线性采样的压缩后的M维的观测向量y进行重构，得到分帧处理后的强震信号x的逼近信号
S5：将重构后的强震信号x的逼近信号
S6：每隔K帧，对重构后的强震信号的逼近信号
对重构后的强震信号
其中，N为分帧长度，表示每一分帧中强震数据的个数；测量数M，即压缩采样后的强震信号的观测数，xi为地震预警系统发送的一帧被分帧处理后的强震信号；
为了进一步减小地震预警系统的信号处理量，从检测信号中准确识别强震信号，本发明的方法中还包括信号捡拾步骤，信号捡拾的方法为：地震预警系统在分帧处理前，对采集到的地震P波信号进行0.1—33Hz的带通滤波，采用STA/LTA方法对带通滤波后的地震P波进行粗略捡拾：
STA为短时间窗，LTA为长时间窗；在地震信号P波到达时STA/LTA会有突变，当两者的比值大于设定阈值R时，判断为有效地震信号：
CFP为地震信号P波的捡拾特征函数，CFP＝Y(k)2+[Y(k)-Y(k-1)]2；Y(k)和Y(k-1)分别为k和k-1时刻垂直方向的加速度记录值；i为当前时刻；k1和k2是当前时刻之间某时刻，且k2＜k1＜i；上式中，公式分子部分STA短时间窗的数据点数为50个，公式分子部分LTA长时间窗的数据点数为1000个数据；
应用STA/LTA方法只能粗略识别地震P波震相的到时，为更加精确对地震P波震相的到时进行捡拾，采用VAR—AIC法进行精确捡拾：
AIC(n)＝n·log{var(x[1,n])}+(L-n-1)·log{var(x[n+1,L])}；
n的范围为地震P波STA或者LTA窗口内所有的采样点；var表示方差；窗内AIC的最小值是到时点，然后将采用STA/LTA方法粗略捡拾的P波数据根据采集时间向前推AIC的最小值相应的点数或者相应的采样时间，实现信号精确捡拾。由于有每秒采样率，故采样点数前移对应了相应的采样时间的前移，这样采集的数据的时效性上更加可靠。
以上所述仅为本发明的较佳实施方式而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
