[发明专利]一种大地电磁信噪辨识及分离方法

权利要求书

1.一种大地电磁信噪辨识及分离方法，其特征在于，包括以下步骤：

步骤1、初始化参数，将大地电磁信号时间序列等间距划分为J段，并对每个大地电磁信号段分别进行去均值处理；

步骤2、对去均值处理后的每个大地电磁信号段x(t)分别进行优化固有时间尺度分解，得到N个旋转分量PR_n(t)和一个呈单调变化趋势的残余分量R_N(t)，其中n＝1,2,…,N；

步骤3、分别提取各个大地电磁信号段经优化固有时间尺度分解得到的N个旋转分量的样本熵和模糊熵；将第j个大地电磁信号段经优化固有时间尺度分解得到的N个旋转分量的样本熵和模糊熵分别表示为SampEn_j和FuzzyEn_j；

步骤4、将SampEn_j和FuzzyEn_j作为第j个大地电磁信号段的联合特征参数，j＝1,2,…,J，对各个大地电磁信号段进行模糊C均值聚类，将其自动划分为有用信号和强干扰信号两大类，实现信噪辨识；

步骤5、将步骤4中辨识为强干扰信号的大地电磁信号段运用小波阈值法进行去噪处理；

步骤6、将步骤5中去噪后的大地电磁信号段与步骤4中辨识为有用信号的大地电磁信号段进行叠加，并加上步骤1中的均值，得到重构后的大地电磁信号；

所述步骤2中，优化固有时间尺度分解的步骤为：

1.1)确定原信号x(t)的极值点序列x_k及其对应的时刻t_k，k＝1,2,…,K；

1.2)利用式计算出极值点序列中每个极值点对应的基线控制点L_k，k＝1,2,…,K，其中0<α<1，t₀＝0；

1.3)对基线控制点L_k进行样条插值拟合，获取基线信号L₁(t)；将原信号x(t)减去基线信号L₁(t)得到第一阶旋转分量PR₁(t)，对应的基线信号L₁(t)称为第一阶残余分量；

1.4)把基线信号作为原信号，重复上述步骤1.1)～1.3)，直到基线信号满足单调条件，终止分解，记此时得到的旋转分量的个数为N，此时的基线信号，即第N阶残余分量为R_N(t)；至此，原信号x(t)被分解为N个旋转分量和呈单调变化趋势的残余分量R_N(t)之和，即

所述步骤3中，N个旋转分量的样本熵的计算方法如下：

2.1)预先给定嵌入维数m和相似容限r，依据N个旋转分量构造一组m维向量，其中第i个向量记为X_m(i)，X_m(i)＝{PR_i(t),PR_i+1(t),...,PR_i+m-1(t)}，i＝1,2,...,N-m+1；

2.2)对每个向量X_m(i)，计算它与其余向量X_m(l),l＝1,2,…,N-m+1且l≠i的距离d_il，并统计d_il小于r的数目及此数目与距离总数N-m的比值，记作其中X_m(i)与X_m(l)之间的距离d_il为两者对应元素差值绝对值的最大值，即d_il＝d[X_m(i),X_m(l)]＝max[|X_m(i+f)-X_m(l+f)|],f＝0,1,…,m-1；

2.3)求的平均值：

B^m(r)表示向量X_m(i)与X_m(l)在相似容限匹配m个点的概率；

2.4)令嵌入维数m＝m+1，重复上述步骤2.1)～2.3)，得到新的B^m(r)，将其记为A^m(r)，A^m(r)表示向量X_m(i)与X_m(l)在相似容限匹配m+1个点的概率；

2.5)根据以下公式计算N个旋转分量的样本熵：

所述步骤3中，N个旋转分量的模糊熵的计算方法如下：

3.1)预先给定嵌入维数m、相似容限r以及相似容限边界梯度v；依据N个旋转分量构造一组向量，其中第i个向量记为s_m(i)，s_m(i)＝{PR_i(t),PR_i+1(t),...,PR_N(t)}-PR_i0，i＝1,2,...,N-m+1，

3.2)计算向量s_m(i)与其余向量s_m(h),h＝1,2,…,N-m+1且h≠i的距离s_m(i)与s_m(h)之间的距离为两者对应元素差值绝对值的最大值，即

3.3)通过模糊函数计算向量s_m(i)与s_m(h)的相似度：

3.4)计算

同时计算

3.5)令嵌入维数m＝m+1，重复上述步骤3.1)～3.4)，得到新的将其记为φ^m(v,r)；

3.6)通过以下公式计算N个旋转分量的模糊熵：

所述步骤4中，模糊C均值聚类的方法为：

4.1)用值在[0，1]间的随机数初始化隶属矩阵U_C×J，使其元素u_cj满足式以下约束条件：

其中，C为聚类个数，u_cj为U_C×J中第c行第j列的元素，表示第j个大地电磁信号段的联合特征参数Q_j属于第c个聚类的隶属度，Q_j＝[SampEn_j,FuzzyEn_j]；

4.2)计算c个聚类中心V_c：

其中，z∈(1,+∞)为加权指数；

4.3)计算目标函数：

其中，d_cj代表Q_j与V_c的欧式距离；

如果目标函数值小于某个确定的阈值，或它相对上次目标函数值的改变量小于某个阈值，则算法结束；否则，进入步骤4.4)；

4.4)根据以下公式更新U_C×J中的元素u_cj，得到新的隶属矩阵U_C×J，并返回步骤4.2)；

其中，d_bj代表Q_j与V_b的欧式距离。