基于SVD-ICA的TSP强工业电干扰压制方法与流程

本发明涉及一种基于svd-ica的tsp工业电干扰压制方法，该方法用于隧道超前地质预报中tsp强工业电干扰的压制。

背景技术：

svd-ica(奇异值分解-独立成成分分析)的tsp(tunnelseismicprediction)工业电干扰压制方法。

tsp是采用多点激发、单点接收的多波多分量高分辨率地震反射波探测技术，在隧道超前地质预报中极其重要。然而tsp采集数据时不可避免会受到50hz的工业电干扰，因此得到的地震数据质量变差，进而严重影响tsp方法速度估计的准确性、反射目标定位可靠性及精度甚至是岩性分析结果，为了解决上述问题，有时只能进行多次采集，这样不仅影响tsp施工效率，造成工期延长甚至由于预报不可靠或者无法预报引起大量不必要的经济损失。目前其他领域对强工业电干扰抑制的方法主要有频域法、时域法、奇异值分解法(svd)和盲源分离法四大类。频域法如陷波法，小波滤波法等都在频率域进行工频压制，但由于tsp有效信号频率和工业电干扰频率混叠，这类方法会造成有用信号能量损失，影响资料信噪比，同时也会令数据处理结果产生一定的相移，影响不良地质体定位精度。时域法如正余弦逼近法，工频回归相减法，自适应滤波法等是将工业电干扰表示为以振幅、频率、相位为变量的函数，再通过正余弦函数逼近法来估计工业电干扰，但这类方法需要对强工业电干扰频率进行精确估计，不适用于频率不稳定的工业电干扰及谐波压制。svd或者主成分分析(pca)法，要求被处理数据为多道数据，而tsp受工业电干扰的数据通常只有一道。盲源分离法如独立成分分析法(ica)，需要满足正定条件，否则无解，且该方法无法确独立成分的次序和极性。可见，上述方法在tsp采集数据中难以应用于工业电干扰压制或压制效果不理想。

技术实现要素：

本发明的目的就在于针对上述现有技术的不足，提供一种基于svd-ica的tsp强工业电干扰压制方法。

本发明的思想是：当tsp采集数据的电缆线附近存在施工用电输电线路时，整个地震波记录上或部分记录道上经常存在50hz的正弦干扰波，有时其强度是有效波信号的数倍，严重影响后续的数据处理及地质预报精度，会给后续隧道施工带来不必要的财产损失。本发明首先对强工业电干扰的单道数据构造hankel矩阵，使其变成二维矩阵形式，然后对构造的hankel矩阵进行svd分解并重构工业电干扰信号，并分析工业电干扰频谱得出其准确频率，最后根据其频率构造对应的正弦和余弦信号做为fastica输入信号，并通过互相关处理fastica的输出分量得到有效信号，该方法能够抑制工业电基频及谐频干扰。

本发明是通过以下技术方案实现的：

基于svd-ica的tsp强工业电干扰的压制方法，包括以下步骤：

a、tsp野外施工采用共检波器多炮点激发的工作形式，采集数据包含x、y、z三分量，并独立存储，将所有炮点记录记为g'，对g'的各道数据逐一进行频谱分析，找出所有受到强工业电干扰的地震道，依次记为p1,p2,...,pg，构建干扰道集，记为

d＝[p1,p2,...,pg](1)

其中g为受强工业电干扰的总道数，ps表示第s个受强工业电干扰的地震道信号，s＝1,2,...,g；

b、从d中取第一道记录p1，令x(k)^t＝p1，则x(k)为行向量，记为

x(k)＝[x(1),x(2),...,x(n)]；(2)

其中k为采样序列，k＝1,2,...,n，n为最大采样点，t为矩阵的转置运算；

c、对x(k)构造hankel矩阵，记为

该矩阵的行数记为m，列数记为n，若n为偶数，则令m＝n/2，n＝n/2+1，若n为奇数，则令m＝(n+1)/2，n＝(n+1)/2；

d、对h进行svd分解，可得

h＝u·s·v^t(4)

其中“·”表示矩阵乘法，u和v都是单位正交矩阵，v^t为v的转置矩阵，s为特征值由大到小排列的对角矩阵，记为s＝diag[σ1,σ2,...,σn]，这里diag代表对角阵，且σ1≥σ2≥...≥σn≥0；

e、若对u和v用列向量集合形式表示，可记为u＝[u1,u2,...,um]，v＝[v1,v2,...,vn]，则

其中ui∈r^m×1，vi∈r^n×1，i＝1,2,...,q，q为m与n中的最小值；

f、利用特征值构造矩阵f和e，即令

表示为

表示为

g、将矩阵f的第一行与最后一列首尾相连合成信号a，记为

a＝[f1,1,f1,2,...,f1,n,f2,n,...,fm,n](10)

将矩阵z的第一行与最后一列首尾相连合成信号b，记为

b＝[e1,1,e1,2,...,e1,n,e2,n,...,em,n](11)

h、对信号a进行傅里叶变换得到对应的幅频曲线，找到该曲线中工频附近峰值频率记为f，若工频谐波干扰较强，根据具体数据处理质量要求，同时找到该曲线中工频谐波附近多个峰值频率记为谐频f2,f3,f4,……；

i、根据f构造对应频率的正弦和余弦信号i2和i3，分别为

i2＝sin(2π·f·k/c)，i3＝cos(2π·f·k/c)(12)

其中c为采样率；

j、令i1＝x(k)，将i1,i2,i3作为fastica的输入信号，记为

k、根据fastica的实现步骤，对矩阵i进行解混，得到解混矩阵w，定义矩阵o＝w·i，形式上记为

l、分别计算o1,o2,o3分量与信号b的相关系数r(b,oj)，即

其中cov表示协方差运算，从三个互相关系数r(b,oj)中选择出绝对值最大的相关系数对应的分量oj0，若r(b,oj0)＞0，令p1'＝oj0，若r(b,oj0)＜0，令p1'＝-oj0，p1'是对p1去工业电基波干扰后的数据；

m、如果p1'中工业电压制效果达到数据处理质量要求，执行步骤n，否则令f＝f2,f3,f4,……，按照步骤i到步骤l去除工业电主要谐波干扰；

n、类似地，对d中p2到pg的所有记录重复步骤c到步骤m，最终得到去除工业电干扰的p2',...,p'g，定义d'＝[p1',p2',...,p'g]，d'为去除工业电干扰后的干扰道集；

o、用d'中的各地震道p1',p2',...,p'g替换g'中对应的p1,p2,...,pg，得到g”,g”即为去除工业电干扰的共检波点道集。

有益效果：tsp是隧道超前探测中探测距离最远，而且探测结果比较有效的方法，其对前方灾害预报将会极大减少施工过程中财产及人员损失。但是tsp采集数据过程中受到50hz的工业电干扰，甚至这种干扰是有效波信号的数倍时，tsp方法难于进行超前地质预报或预报精度下降。经验证，本发明公开的基于svd-ica的tsp工业电干扰压制方法能够实现对工业电基频及谐波干扰的有效压制，特别是对强工业电干扰的压制，与传统的工业电干扰压制方法相比，噪声压制修改好，且对有效信号损害小，处理过程无相位偏移，不影响目标定位精度。该方法处理数据快，可减少由于工业电干扰过大时重复采集的费用，从而提高施工效率和节约施工成本，在不改变采集施工条件下，提高数据质量，对改善地质预报精度具有重要意义。

附图说明：

图1为含强工业电干扰的tsp记录；

图2去除工业电噪声后的tsp记录与理想不含强工业电干扰记录对比

(a)时域波形细节对比；

(b)频域特征细节对比。

具体实施方式：

下面结合附图和实施例做进一步的详细说明：

在本实施例中使用三分量记录中x分量的一道数据进行处理，记录时间为382ms，采样率1000hz。

1、基于svd-ica的tsp强工业电干扰的压制方法，包括以下步骤：

a、tsp野外施工采用共检波器多炮点激发的工作形式，采集数据包含x、y、z三分量，并独立存储，将所有炮点记录记为g'，对g'的各道数据逐一进行频谱分析，找出所有受到强工业电干扰的地震道，依次记为p1,p2,...,pg，构建干扰道集，记为

d＝[p1,p2,...,pg](1)

其中g为受强工业电干扰的总道数，ps表示第s个受强工业电干扰的地震道信号，s＝1,2,...,g；

b、从d中取第一道记录p1，令x(k)^t＝p1，则x(k)为行向量，记为

x(k)＝[x(1),x(2),...,x(n)]；(2)

其中k为采样序列，k＝1,2,...,n，n为最大采样点，t为矩阵的转置运算，本例中p1为x分量中的1道数据，n＝382，有效信号频谱范围为20～100hz；

c、对x(k)构造hankel矩阵，记为

该矩阵的行数记为m，列数记为n，若n为偶数，则令m＝n/2，n＝n/2+1，若n为奇数，则令m＝(n+1)/2，n＝(n+1)/2，本例中m＝192，n＝191；

d、对h进行svd分解，可得

h＝u·s·v^t(4)

其中“·”表示矩阵乘法，u和v都是单位正交矩阵，v^t为v的转置矩阵，s为特征值由大到小排列的对角矩阵，记为s＝diag[σ1,σ2,...,σn]，这里diag代表对角阵，且σ1≥σ2≥...≥σn≥0；

e、若对u和v用列向量集合形式表示，可记为u＝[u1,u2,...,um]，v＝[v1,v2,...,vn]，则

其中ui∈r^m×1，vi∈r^n×1，i＝1,2,...,q，q为m与n中的最小值；

f、利用特征值构造矩阵f和e，即令

可表示为

可表示为

g、将矩阵f的第一行与最后一列首尾相连合成信号a，记为

a＝[f1,1,f1,2,...,f1,n,f2,n,...,fm,n](10)

将矩阵z的第一行与最后一列首尾相连合成信号b，记为

b＝[e1,1,e1,2,...,e1,n,e2,n,...,em,n](11)

h、对信号a进行傅里叶变换得到对应的幅频曲线，找到该曲线中工频附近峰值频率记为f，若工频谐波干扰较强，根据具体数据处理质量要求，同时找到该曲线中工频谐波附近多个峰值频率记为谐频f2,f3,f4,……，本例中f＝50.01hz；

i、根据f构造对应频率的正弦和余弦信号i2和i3，分别为

i2＝sin(2π·f·k/c)，i3＝cos(2π·f·k/c)(12)

其中c为采样率；

j、令i1＝x(k)，将i1,i2,i3作为fastica的输入信号，记为

k、根据fastica的实现步骤，对矩阵i进行解混，得到解混矩阵w，定义矩阵o＝w·i，形式上记为

l、分别计算o1,o2,o3分量与信号b的相关系数r(b,oj)，即

其中cov表示协方差运算，从三个互相关系数r(b,oj)中选择出绝对值最大的相关系数对应的分量oj0，若r(b,oj0)＞0，令p1'＝oj0，若r(b,oj0)＜0，令p1'＝-oj0，p1'是对p1去工业电基波干扰后的数据，本例中r(b,o1)＝0.9196，r(b,o2)＝-0.0416，r(b,o3)＝-0.2411，本例中选择p1'＝o1，p1'是对p1去工业电干扰后的数据；

m、如果p1'中工业电压制效果达到数据处理质量要求，执行步骤n，否则令f＝f2,f3,f4,……，按照步骤i到步骤l去除工业电主要谐波干扰，本例中谐波干扰弱，因此只进行了工业电基频干扰压制；

n、类似地，对d中p2到pg的所有记录重复步骤c到步骤m，最终得到去除工业电干扰的p2',...,p'g，定义d'＝[p1',p2',...,p'g]，d'为去除工业电干扰后的干扰道集；

o、用d'中的各地震道p1',p2',...,p'g替换g'中对应的p1,p2,...,pg，得到g”，g”即为去除工业电干扰的共检波点道集。