一种信号恢复的优化方法 【技术领域】
本发明属于信号处理领域， 尤其适用于应用在地球物理、 通讯工程、 电子工程、 军 事等多个应用领域。背景技术
现有技术中， 基于双谱运算的信号重构， 最重要的一个步骤是对基频率 ( 或基波 数 )f0 的准确估算。现有的方法为 ： 取双谱的对角线值 B(f1， f2)， 使 |B(f， f)| 不为 0 的最 小 f 便是估算的 f′ 0。f′ 0 为 f0 的初次估值， 不一定精确， 特别是当基频 f0 与频率分辨 率 Δf 之比不为整数时。估算的 f′ 0 与真值 f0 的误差范围为 f′ 0 与真值 f0 的误差， 进一步搜索 为减小估算值区间内使 |B(f， f)| 达到最大的 f， 进行频率细化处理。对于实际数据， 由于噪声的影响， |B(f， f)| 不可能为 0， 采用第一个谱峰 来代替不为 0 的 |B(f， f)|。 在进行频率细化处理的搜索计算时， 先设定误差控制的精度 ε。 搜索计算采用迭代算法 ：
(1) 令 B ＝ |B(f′ 0， f′ 0)|，计算 |B(f， f)| ； 若 Δf ＜ ε， 则结(2) 若 B ＜ |B(f， f)|， 则令 B ＝ |B(f， f)|， f′ 0 ＝ f，束， 否则转回 (1) ； 若 B ＞ |B(f， f)|， 则转入 (3)。
(3) 令 B ＝ |B(f′ 0， f′ 0)|，计算 |B(f， f)| ； 若 Δf ＜ ε， 则结束，若 B ＜ |B(f， f)|， 则令 B ＝ |B(f， f)|， f′ 0 ＝ f，否则转回 (1) ； 若 B ＞ |B(f， f)|， 则令
若 Δf ＜ ε， 则结束， 否则转回 (1)。 区间内， 迭代搜索使为减小估算值 f ′ 0 与真值 f0 的误差， 在|B(f， f)| 达到最大的 f， 作为实际计算得到的基频 f0。 现有的技术中存在以下 3 方面的技术问题 ：
(1) 通常情况下， 信号的基频 f0 与频率分辨率 Δf 之比不为整数， 因此现有技术的 方法估算的 f0， 在信噪比较低时是个近似值 ( 如表 1) ； 即使基频 f0 与频率分辨率 Δf 之比 为整数， 当信噪比较低时， 上述方法估算的 f0 也不能确保得到真值 ( 如表 2)。为使基频 f0 估算的精度足够高， 迭代运算的次数、 计算量随之大大增加。
表 1 半基频 f0/2 不为频率分辨率 Δf 的整数倍时不同信噪比信号估算的 f′ 0 及 其双谱值 |B(f′ 0， f′ 0)|
表 2 半基频 f0/2 为频率分辨率 Δf 的整数倍时不同信噪比信号估算的 f′ 0 及其 双谱值 |B(f′ 0， f′ 0)|(2) 现有技术的方法， 频率分辨率 Δf 越高， 则基频 f0 估算的精度越高， 但对于实 际信号， 计算用的信号长度与信号的时间域采样率已经固定， 也即频率分辨率 Δf 已经固 定， 因此， 增加迭代运算的次数是提高估算精度的唯一途径。
(3) 现有技术由于在方法上存在缺陷而带来实际应用上的局限。目前仅用于具有 周期变化的信号的恢复处理， 例如， 仿真试验用的模型信号是具有周期变化的谐波信号， 实 际信号试验用的也是具有周期变化的缝纫机振动信号。
为了克服现有技术中基频 f0 估算不准的缺陷， 增强 “应用双谱运算进行信号恢复” 这一技术的适用性和实用性， 申请人提出了本发明。
发明内容
本发明为了克服上述现有技术中存在的技术问题， 研发了一种信号恢复的优化方 法。本发明使基于双谱运算的信号重构技术适用于任意信号数据的重构， 而且重构信号的 保真度高。
本发明的核心内容是， 通过信号处理方案、 流程的改变， 规避现有技术中方法的缺 陷。 具体地说， 就是在对原始信号做双谱运算前， 先对原始信号做 3 项预处理 ： 去线性背景 ； 零均值化 ； 周期拓展。
其中， 最关键之处是， 对原始信号进行周期拓展处理。因为周期拓展， 使信号的基 频 f0 与频率分辨率 Δf 相等。 这样， 不管原始信号的信噪比如何， 基频 f0 与频率分辨率 Δf 的比值始终是整数 1， 从而， 确保了估算得到的基频 f0 是真值。具体实现时， 通过搜索双谱 值 |B(f， f)| 的第一个谱峰时的频率 f 即可得到基频 f0 的真值。为使后面应用式 (14) 推 算谐波项系数 ai 与相位角 时， 方便、 准确地搜索、 应用双谱振幅谱 By(f0/2， f0/2) 和相位谱 φy(f0/2， f0/2) 的数值， 根据经验将原始信号拓展成 4 个周期比较合适。
具体方法如下，
一种信号恢复的优化方法， 所述方法在对地球物理勘探采集信号做双谱运算前， 先对原始信号做预处理过程 ； 所述预处理过程包括去线性背景、 零均值化、 周期拓展 ； 完成 预处理过程的信号再经过双谱运算过程后， 实现了原始信号的完全恢复 ： 完整描述一个信号的谐波信号基频、 谐波项系数与相位等 3 项参数被准确求取。
所述方法包括如下步骤 ：
步骤 1 读取输入离散化后地球物理采集信号 x(t) 的数据 ： x(i)， i ＝ 1，…， n； 其 中， n 为信号 x(t) 的采样点数 ；
步骤 2 对信号 x(t) 进行预处理过程 ：
(1) 去线性背景过程 ：
利用以下公式计算信号 x(t) 的线性背景斜率 bg ：
式中， x(n) 为信号 x(t) 的最末一个数据， x(1) 为信号 x(t) 的第一个数据， n 为信 号 x(t) 的采样点数 ；
然后利用以下公式消去信号的线性背景 ：
xb(i) ＝ x(i)-x(1)-bg×(i-1)
式中， xb(i)， i ＝ 1，…， n 为信号 x(t) 消去线性背景后的数据， x(i)， i ＝ 1，…， n 为原始信号 x(t) 的数据， x(1) 为原始信号 x(t) 的第一个数据 ；
(2) 零均值化过程 ：
对信号 x(t) 消去线性背景后的数据 xb(i)， i ＝ 1，…， n 的各个样点数据 xb(i)， 求和并除以样点数 n， 得信号去线性背景后数据 xb(i)， i ＝ 1，…， n 的平均值
式中， n 为信号 x(t) 的采样点数 ； 去线性背景后的数据 xb(i)， i ＝ 1，…， n 减去平均值 实现零均值化 ：式中， xe(i)， i ＝ 1，…， n 为信号 x(t) 去线性背景后又零均值化后的数据 ；
(3) 周期拓展过程 ：
原信号是周期数为 1 的信号， 将原信号首尾相连重复放置于后， 重复放置 K-1 次， 即得周期数为 K 的新信号 ；
周期拓展过程的 C 语言计算程序语句如下 ：
for(i ＝ 0 ； i＜k； i++)for(j ＝ 0 ； j＜n； j++)xn[j+i*k] ＝ xe[j] ；
其中， xe[j] 为原始信号 x(t) 去线性背景后又零均值化后数据 xe(i)， i ＝ 1， …， n 的数组， xn[j+i*k] 为周期拓展后信号 xn(t) 的数组， n 为信号 x(t) 的采样点数， k 为周 期拓展后信号的周期数 ；
步骤 3 通过 Fourier 变换， 计算周期拓展后信号 xn(t) 的频谱 X(f) ： jφ(f)
X(f) ＝ A(f)·e
其中， A(f) 和 φ(f) 分别为频谱 X(f) 的幅值和相位 ； 其中， 频谱 X(f) 的长度为 k·n ； 其中， n 为信号 x(t) 的采样点数， k 为周期拓展后信号的周期数 ；
步骤 4 根据下式， 通过三重相关计算信号 xn(t) 的双谱 B(f1， f2) ：
其中， X(f1)、 X(f2)、 X(f1+f2) 分别是信号 xn(t) 在频率变量 f1， f2， f1+f2 上的频谱， X (f1+f2) 是频谱 X(f1+f2) 的共轭， A(f1， f2) 是双谱 B(f1， f2) 的幅值， φ(f1， f2) 是双谱 B(f1， f2) 的相位 ； 其中 f1， f2 是与时间变量 t1， t2 相对应的频率变量， t1， t2 表示两个不同的时间 变量 ；
这一过程的 C 语言计算程序语句如下 ：
for(i ＝ 0 ； i ＜ km ； i++){
for(j ＝ 0 ； j ＜＝ i ； j++){
if(i+j ＞＝ km)break ；
ba[i][j] ＝ a[i]*a[j]*a[i+j] ；
bf[i][j] ＝ f[i]+f[j]-f[i+j] ；
}
}
其中， 数组 a[i]、 f[i] 分别是信号 xn(t) 的频谱 X(f) 的振幅谱和相位谱数组， 数组 ba[i][j]、 bf[i][j] 分别是信号 xn(t) 的双谱 B(f1， f2) 的振幅谱 A(f1， f2) 和相位谱 φ(f1， f2) 数组 ； 其中 km 信号 xn(t) 的频谱 X(f) 的长度的 1/2， 即 km ＝ k·n/2+1 ；
步骤 5 利用对称性计算整个第一象限的双谱 B(f1， f2) ：
B(f2， f1) ＝ B(f1， f2)
其中， f1， f2 的含义同上 ；
这一过程的 C 语言计算程序语句如下 ：
for(i ＝ 0 ； i ＜ km ； i++){
for(j ＝ 0 ； j ＜＝ i ； j++){
if(i+j ＞＝ km)break ；
ba[j][i] ＝ ba[i][j] ；
bf[j][i] ＝ bf[i][j] ；
}
}
其中， 数组 ba[i][j]、 bf[i][j] 分别是信号 xn(t) 的双谱 B(f1， f2) 的振幅谱数组 和相位谱数组 ； 其中 km 信号 xn(t) 的频谱 X(f) 的长度的 1/2， 即 km ＝ k·n/2+1 ；
步骤 6 取信号 xn(t) 的双谱 B(f1， f2) 的模 |B(f， f)| 的第一个峰值时的 f 为估算 的基频 f0 ；
步 骤 7 在 信 号 xn(t) 上 加 一 相 位 为 0 的 半 基 频 余 弦 信 号，即 y(t) ＝ xn(t)+cos(πf0t) ；
步骤 8 同步骤 3， 通过 Fourier 变换， 计算信号 y(t) 的频谱 Y(f) ： jφ(f)
Y(f) ＝ A(f)·e
其中， A(f) 和 φ(f) 分别为频谱 Y(f) 的幅值和相位 ； 频谱 Y(f) 的长度为 k·n ； 其中， n 为信号 x(t) 的采样点数， k 为周期拓展后信号的周期数 ；
步骤 9 同步骤 4， 通过三重相关计算信号 y(t) 的双谱 By(f1， f2) ；
步骤 10 同步骤 5， 利用对称性计算整个第一象限的双谱 By(f1， f2) ；
*步骤 11 设所要恢复的信号 x(t) 可以用谐波信号来描述， 即
其中， f0 是谐波基频， 已在步骤 6 求得 ； ai， i ＝ 1，…， M， 分别是谐波项系数、 相 位角， 可以根据下式由信号 y(t) 的双谱 By(f1， f2) 的振幅谱 Ay(f1， f2) 与相位谱 φy(f1， f2) 求得 ：
(i ＝ 2， 3，…， M)
(i ＝ 2， 3，…， M)其中， M 是谐波项数， 可由下式求得 ：
M ＝ km/nf
其中， km 为信号 xn(t) 的频谱 X(f) 的长度的 1/2， 即 km ＝ k· n/2+1 ； nf 是估算得 到的基频 f0 所对应的频率采样序号 ；
步骤 12 将基频 f0、 谐波项系数 ai、 谐波项相位角 谐波项数 M， 一起代入描述信号 x(t) 的谐波信号表达式 ：
得重构信号 ； 步骤 13 输出信号 x(t) 恢复结果。 本方法在步骤 12 后还包括进行线性背景补偿和均值补偿处理得恢复信号步骤 ； 所述均值补偿处理过程为 ： 对输出信号 x(t)， 加上步骤 2 中零均值化所求得的平均值 即所述线性背景补偿过程为 ：
对均值补偿处理后的信号 x′ (t)， 再加上步骤 2 中去线性背景时所求得的线性背 景斜率 bg ：
x″ (i) ＝ x′ (i)+x(1)+bg×(i-1)
其中， x″ (i)， i ＝ 1，…， n 为最终恢复得到的信号 x″ (t) 的数据， x′ (i)， i＝ 1，…， n 为输出信号 x(t) 作了均值补偿处理后的数据， x(1) 为原始信号的第一个数据， bg 为原始信号的线性背景斜率。
本发明使基频 f0 的估算既简便又准确， 而且使基于双谱运算的信号恢复这一技 术， 能够适用于任意信号的信号恢复处理。 本发明为信号恢复、 图像降噪等处理提供了一个 有用的工具。通过对现有技术文献中相同的谐波信号模型、 地球物理勘探中的重力测量数 据、 地球物理勘探中的地震子波信号、 地球物理勘探中的实际地震数据等多种信号数据， 应用本发明研制的技术进行信号重构计算， 基频 f0 估算准确， 重构信号的保真度高， 方法技术 的适用性和实用性强。 附图说明
图 1 为本发明方法流程框图 ；
图 2 为地球物理勘探中的重力测量信号的恢复处理结果 ；
图 3 为地球物理勘探中的地震子波信号恢复处理的结果 ；
图 4 为地球物理勘探中的原始地震信号剖面 ；
图 5 为应用本发明的方法恢复处理的地震信号剖面 ；
图 6 为放大 100 倍显示的附图 4 与附图 5 的差值。
下面将结合具体实施方式对各幅附图进行说明 具体实施方式
本发明所采用的方法流程详细说明 ：
第一步， 读取离散化后地球物理信号 x(t) 的数据文件。 第二步， 对信号 x(t) 进行去线性背景、 零均值化、 周期拓展等预处理。
第三步， 通过 Fourier 变换， 计算信号 x(t) 的频谱 X(f)。
第四步， 根据式 (4)， 通过三重相关计算信号 x(t) 的双谱 B(f1， f2)。
第五步， 利用对称性计算整个第一象限的双谱 B(f1， f2)。
第六步， 取 |B(f， f)| 的第一个峰值时的 f 为估算的基频 f0。
第 七 步， 在 原 信 号 x(t) 上 加 一 相 位 为 0 的 半 基 频 余 弦 信 号， 即 y(t) ＝ x(t)+cos(πf0t)。
第八步， 通过 Fourier 变换， 计算信号 y(t) 的频谱 Y(f)。
第九步， 根据式 (4)， 通过三重相关计算信号 y(t) 的双谱 By(f1， f2)。
第十步， 利用对称性计算整个第一象限的双谱 By(f1， f2)。
第十一步， 根据式 (14) 计算谐波项系数 ai 与相位角 第十二步， 谐波项系数 ai 与相位角 代入式 (5)， 得重构信号。 第十三步， 进行线性背景补偿和均值补偿处理得恢复信号。 第十四步， 输出信号恢复结果。 在具体各个步骤中， 依次为 (1) 基于双谱运算信号恢复的方法原理 设 x(t) 是一均值为 0 的实信号， 其 Fourier 变换为 ： X(f) ＝ A(f)·ejφ(f) (1) A(f) 和 φ(f) 分别为频谱 X(f) 的幅值和相位。 x(t) 的三阶累积量函数表达式为 ：t1， t2 表示两个不同的时间变量。 C(t1， t2) 的二维 Fourier 变换即为双谱 B(f1， f2)X*(f1+f2) 是频谱 X(f1+f2) 的共轭 ( 注 ： 下同 )。其中 f1， f2 是与时间变量 t1， t2 相 对应的频率变量。
进一步， 可得 x(t) 的双谱表达式 ：
设所要恢复的信号 x(t) 可以用谐波信号来描述， 即
根据 Fourier 变换的特性， 由式 (4) 和 (5) 可以得到双谱振幅谱、 相位谱与谐波项 系数、 相位角的关系式 ： (i ＝ 1， 2，…， M/2， j ＝ i， i+1，…， M-i) (i ＝ 1， 2， …， M/2， j ＝ i， i+1， …， M-i) (6a) (6b)
方程式 (6a)、 (6b) 中的未知数也即式 (5) 中的未知数， 有基频 f0、 谐波项系数 ai 与相位角 (i ＝ 2， 3，…， M)。显然， (6a)、 (6b) 中未知数的个数大于方程数。
实际应用时须先对基频 f0 进行准确地估算， 这是本方法成功实现的关键， 而且，
由式 (5) 可知， 基频 f0 估算的准确度， 决定了恢复信号的保真度。
根据估算得到的基频 f0， 在信号 x(t) 上加一相位为 0 的半基频余弦信号
z(t) ＝ cos(πf0t)， 即
y(t) ＝ x(t)+z(t) ＝ x(t)+cos(πf0t) (7)
由于 y(t) 的 Fourier 变换可以表示为 ：
Y(f) ＝ X(f)+Z(f) (8)
故 y(t) 的双谱可以表示为 ：
By(f1， f2) ＝ Y(f1)·Y(f2)·Y*(f1+f2)
＝ [X(f1)+Z(f1)]·[X(f2)+Z(f2)]·[X*(f1+f2)+Z*(f1+f2)] (9)
对于式 (9)， 如果 f1 ＝ f2 ＝ f0/2， 则根据 Fourier 变换的特性， 有
从而有
如果 f1、 f2 均不为 f0/2， 且 f1+f2 ≠ f0/2， 则有 *
Z(f1) ＝ Z(f2) ＝ Z (f1+f2) ＝ 0 (12)
从而， 由式 (12)、 (9) 可得
By(f1， f2) ＝ X(f1)·X(f2)·X*(f1+f2) ＝ Bx(f1， f2) (13)
即信号 y(t) 的双谱此时就是信号 x(t) 的双谱。从而可得所求谐波信号表达式 x(t) 中的谐波项系数 ai 与相位角 求解的递推公式 ：
(i＝2， 3，…，M)(14a)
(i＝2， 3，…，M)(14b) 由此， 即实现了对信号 x(t) 的恢复。
根据上述技术流程， 发明人在 Windows 操作系统下利用 C++Builder， 依照图 1 的程 序框图， 实现了采用本发明所述的方法做的任意信号的恢复处理， 图中虚框部分为本发明 的核心内容。
(1) 地球物理勘探中的重力测量信号的恢复处理 ( 附图 2)。除了个别的单点异常 外， 原始信号几乎得到了完全的恢复。 这些个别的单点， 则反映了信号采集时的随机噪声的 存在。说明本发明的信号恢复方法， 可以实现对信号中随机噪声的消除。
(2) 地球物理勘探中的地震子波信号的恢复处理 ( 附图 3)。地震子波的主频为 202.4Hz。这一信号近乎脉冲信号， 但本发明的方法也使它实现了信号的完全恢复， 恢复信 号与原始模型信号的均方误差为 0.00076106。
(3) 地球物理勘探中的实际地震记录信号的恢复处理 ( 附图 4 ～附图 6)。从附图 5 与附图 4 的对比来看， 原始地震记录信号得到了完全的恢复， 但从附图 6 上来看， 恢复信号 与原始信号的差值是一些在 “0” 值附近来回震荡的随机噪声。说明应用本发明的方法可以 使图像中的随机噪声有效地消除。