池 越1,2，丁 木1,2，周亚同1,2，白 阳1,2
(1.河北工业大学 电子信息工程学院，天津 300401; 2.河北工业大学 天津市电子材料与器件重点实验室,天津 300401)
摘 要:地震信号重建具有重要意义，其中傅里叶重建算法受到广泛关注。这些算法具有原理简单、纯数据驱动、假设前提少等优势，但也存在着重建速度慢等不足。以FRMN为基础，构建了一种新的Landweber迭代傅里叶快速重建算法(Fourier Reconstruction with accelerated Landweber Iteration，FRLI);该算法在经典的FRMN基础上，利用快速Landweber迭代求解傅里叶系数，再经过傅里叶反变换重建地震信号。通过不同维数的地震信号重建实验，结果表明:FRLI与经典的FRMN和FRSI重建算法相比，大幅提高了重建速度，且重建效果良好。因此采用FRLI算法重建地震信号可以解决传统算法的缺陷。
关键词:地震信号;傅里叶重建;快速Landweber迭代;重建速度
中图分类号:P631.4
文献标志码:A
文章编号:0253-9993(2018)09-2562-08
移动阅读
池越，丁木，周亚同，等.地震信号的Landweber迭代傅里叶快速重建[J].煤炭学报,2018,43(9):2562-2569.doi:10.13225/j.cnki.jccs.2017.1482
CHI Yue,DING Mu,ZHOU Yatong,et al.Fourier reconstruction of seismic signal with accelerated Landweber iteration[J].Journal of China Coal Society,2018,43(9):2562-2569.doi:10.13225/j.cnki.jccs.2017.1482
收稿日期:2017-10-29
修回日期:2018-07-27
责任编辑：韩晋平
基金项目:国家自然科学基金资助项目(61401307);河北省科学技术研究与发展基金资助项目(11213565);河北省引进留学人员联合资助项目(CL201707)
作者简介:池 越(1977—)，男，河北唐山人，副教授，博士。E-mail:chiyueliuxin@126.com
通讯作者:周亚同(1973—)，男，湖北荆州人，教授，博士。E-mail:zyt@hebut.edu.cn
CHI Yue1,2,DING Mu1,2,ZHOU Yatong1,2,BAI Yang1,2
(1.School of Electronics and Information Engineering,Hebei University of Technology,Tianjin 300401,China; 2.Tianjin Key Laboratory of Electronic Materials and Devices,Hebei University of Technology,Tianjin 300401,China)
Abstract:The reconstruction of seismic signals is of great significance,and the Fourier reconstruction algorithm has received extensive attention.These algorithms have some advantages of simple principle,pure data drive and little assumption.However,there are also the shortcoming of slow computation speed.Therefore,on the basis of FRMN,the new algorithm of Fourier Reconstruction with accelerated Landweber Iteration (FRLI) was put forward.Based on the classical FRMN,the algorithm uses accelerated Landweber iteration to solve the Fourier coefficients.Then the seismic signals were reconstructed through the inverse Fourier transform.Through the reconstruction experiments of seismic signals with different dimensions,the results show that the FRLI greatly improves the reconstruction speed and the reconstruction effect is better than the classic FRMN and FRSI reconstruction algorithms.Therefore,the FRLI algorithm for reconstructing seismic signals can solve the defects of traditional algorithms.
Key words:seismic signal;fourier reconstruction;accelerated landweber iteration;reconstruction speed
通过地震信号可以探查地下地质结构，以便准确寻找煤炭[1]、石油和天然气。但是在野外地震勘探中，由于各种原因有时导致采集的地震信号并不完整，因此需要将缺失的地震信号重建成完整的地震信号。常见的地震信号重建算法有3类:基于数学变换的重建、基于预测滤波的重建和基于波动方程的重建，其中基于数学变换的重建包括傅里叶变换、小波变换[2]、K-L变换及Radon变换重建等。近年来，随着稀疏信号约束及压缩感知理论的发展，上述理论已被应用于地震信号重建[3-4]。
当前地震信号重建算法中，傅里叶重建受到广泛关注。经典算法包括最小范数傅里叶重建(Fourier Reconstruction with Minimum Norm，FRMN)[5-6]和稀疏反演傅里叶重建(Fourier Reconstruction with Sparse Inversion，FRSI)[7-8]等。FRMN和FRSI的基本思路是:时域中非均匀采样信号经过傅里叶变换到频域，得到均匀采样的傅里叶系数，再由傅里叶系数经过傅里叶反变换得到均匀采样的时域信号。针对大道距地震信号重建问题，文献[9]对地震信号的大道距位置添加预测误差算子的滑动窗口，滑动窗口采用最小二乘加权方式去填补空白地震道;文献[10]中提出用正则化功率谱约束代替稀疏约束，使得快速傅里叶变换式自适应修改迭代约束项，能有效避免空间假频;文献[11]提出辅助插值FRSI算法，先在大间距处插值作为辅助地震道，再由FRSI重建;文献[12]引入稀疏约束矩阵，提出可变稀疏度傅里叶重建算法，利用贝叶斯推理解释并验证其合理性。
因此，对傅里叶重建算法的深入研究具有重要意义。Landweber迭代是一种优秀的不定方程组求解算法，具有较快的计算速度，已经被用于图像处理[13]、信号重建、计算机视觉等多个领域。例如阈值Landweber迭代已被应用于磁感应成像(MIT)图像重建[14]，文献[15]将小波域加速Landweber迭代法用于图像恢复。Landweber预迭代和改进粒子群算法已联合起来用于电阻层析成像图像重建[16-17]。广义逆预迭代快速投影Landweber已被应用于电容层析成像(Electrical Capacitance Tomography)图像重建[18]。将FFT和Landweber迭代相结合，可实现带限信号非均匀釆样重建[19-20]。另外，Landweber在语音重构中也有应用[21]，还可利用Landweber迭代分析滤波反向投影特征[22]，而且在数学领域应用更为广泛[23-24]。
Landweber迭代法具有优势，可以替代傅里叶重建算法中的迭代，分析常见傅里叶重建算法，将Landweber迭代与其结合。FRMN的重建精度低，FRSI的重建精度高但重建速度慢，而Landweber迭代具有较快的计算速度。因此，笔者将FRMN和Landweber迭代二者结合起来，提出一种新的基于Landweber迭代的地震信号傅里叶快速重建算法。该算法在FRMN的基础上，利用Landweber迭代求解傅里叶系数，再重建规则化等间距矩阵，最后通过傅里叶反变换重建地震信号，应用效果良好。
若信号采样点是非均匀的，设非均匀采样点位置为[x0,x1,x2,…,xN-1]，对应的非均匀离散傅里叶变换(Non-uniform Discrete Fourier Transform，NDFT)为
Δ
Δxn
(1)
Δkx]e-jmΔkxxn
(2)
其中,Δxn=(xn+1-xn-1)/2表示采样间隔。频域带宽限制在[-MΔkx,MΔkx]，其中Δkx=2π/(xN-x0)。
FRMN建立在NDFT的基础上，首先将式(2)表示为矢量的形式
(3)
式(3)中傅里叶系数
是未知的，需要由矢量y估计出
在实际应用中信号并没有精确的带宽限制，带宽外也存在部分频率信息，带宽外的信号作为噪声处理，n表示噪声，常将式(3)写成:
(4)
文献[5]通过最大后验概率估计(MAP)推出的计算式，设噪声为高斯白噪声n=N(0,Cn)(Cn为高斯白噪声协方差矩阵)，先验系数
,
，由MAP估计出傅里叶系数,即
(5)
其中,
可以看作是的协方差矩阵;
，
，I为单位阵;
和
分别为信号和噪声的方差;与
的平方最大值成一定比例。Cn由对角阵c2W-1代替，且Wii=Δxi，则式(5)可以简化为
(6)
/
(7)
其中,Mp=2M+1代表频域上采样点数，也是傅里叶系数总个数;F代表预期信噪比;N为非规则地震数据的总道数。式(6)给出了FRMN的非均匀采样信号频谱。
FRSI和FRMN都是经典频率域非均匀采样信号重建算法，其中FRSI是FRMN的一种改进。FRSI使下面目标函数最小化以求解傅里叶系数
ρ
(8)
式(8)右边第1项为数据，第2项为罚函数。Zwartjes通过Z变换[8]，将式(8)转化为
′
(9)
其中,
，A′
通过贝叶斯推理使式(9)的估计量为
′HWA′
(10)
然后，将代入式(10)得
(11)
其中，
即为FRSI的非均匀采样信号频谱。
其中
，为稀疏平滑因子，a和
用来调节稀疏度，a选取不同对应不同的罚函数模型。使用最小二乘加权迭代法(IRLS)求解矩阵，并在迭代中采用共轭梯度法(CGNE)，通过CGNE反复计算将非线性系统线性化。
在求出傅里叶系数后，通过傅里叶反变换实现信号重建
(12)
其中,Ar为规则化的等间距Δx重建矩阵，其定义与矩阵A的定义类似。
通过理论分析可知FRMN的优势在于理论简单，FRSI理论复杂但重建精度高，因此需要对傅里叶重建算法进行深入研究，以便找到一种理论简单且重建精度较高的算法。
对于式(13)所列方程组，其中A:x→y为有界线性矩阵
(13)
已知矩阵A和y，求解向量当A可逆时，也就是A为非奇异阵，则
，可以通过矩阵计算直接求出。但很多情况下A不可逆，或A根本就不是一个方阵，这时式(13)求解是一个不适定问题，Landweber迭代是求解不适定问题的有效手段[25]，其迭代式为
(14)
其中0<b<1/‖A‖2。迭代终止条件为
/‖y‖≤δ1，
/
≤δ2，其中δ1和δ2为任意小数。
张军以级数形式分解迭代公式，对Landweber迭代进行了改进[26]，从而大幅减少了迭代次数，提高了收敛速度。Landweber迭代给定初始条件为p0=bA*y，A0=bA*，将式(14)以级数形式表示
(15)
式(15)中的
由A的迭代式来表示
(16)
(17)
这样
的迭代转化为Al的迭代，式(16)和(17)与式(14)是等价的。当终止条件相同时，迭代次数却有很大不同，由式(14)和(17)比较可以得出:Al迭代l次相当于迭代2l-1次，从而大幅减少迭代次数，有效提高了收敛速度。
针对FRMN算法中的式(3)，y为n×1的向量，A为n×m矩阵，为m×1的向量，并且已知n>m，则式(3)为超定方程，所以可以利用快速Landweber迭代求解该方程组，最终得出傅里叶系数序列
地震信号的Landweber迭代傅里叶快速重建(FRLI)步骤如下:
(1)给定地震信号
;
(2)设置频域带宽参数M，最大迭代次数L，迭代终止条件δ1和δ2;
(3)FRMN的非均匀采样信号频谱如式(6)所示，而FRLI的频谱为
，这是两种算法的不同之处。因此只需先求出矩阵A，然后通过Landweber迭代计算，即可求出FRLI的频谱;
(4)给定频谱的初始值:
，A0=bA*，其中0<b<1/‖A‖2;
(5)计算矩阵Al+1=Al(2I-AAl);
(6)得出矩阵后计算向量
，如果满足迭代终止条件
/‖y‖≤δ1和/≤δ2则迭代结束;否则l=l+1，转至步骤(5);
(7)计算重建矩阵Ar。A由非均匀采样信息[x0,x1,x2,…,xN-1]及它们的间距Δxi计算出来，而重建矩阵Ar的计算只需把这些非均匀分布的采样点置换成均匀等间隔网格，即可得重建矩阵;
(8)由
进行傅里叶反变换，即可从频率域中恢复到空间域得到地震重建信号。
Landweber迭代傅里叶快速重建算法中，影响重建效果的主要参数有:频域带宽参数M，最大迭代次数L，迭代终止条件δ1,δ2。L通常选取1～100次就可达到收敛效果。
重建实验在一台联想Think Pad笔记本上进行，处理器为i5-3320M，CPU为2.60 GHZ，内存为4 GB，64位Windows操作系统，安装有MATLAB-R2014b。每个重建算法都是通过多组实验，调整参数反复比较，选取一组参数使得重建效果最佳。通过重建均方误差(MSE)、重建误差极值(EE)及相对重建误差极值(REE)，还有重建时间等量化指标，比较FRMN、FRSI和FRLI算法的性能。
(18)
EE=max[max(|yr-y|)]
(19)
(20)
该叠后地震信号道间距10 m，共计200道，时间轴上采样率为4 ms，每道含700个均匀采样点，记录长度2 800 ms，原始地震信号如图1(a)所示;图1(b)为人工随机删除50道后的信号，缺失率为25%。现采用3种重建算法对之重建，相关参数选择见表1。
FRMN重建结果如图1(c)所示，图1(d)为FRMN重建结果与原始地震信号的差值;FRSI重建结果如图1(e)所示，图1(f)为FRSI重建结果与原始地震信号的差值;FRLI重建结果如图1(g)所示，图1(h)为FRLI重建结果与原始地震信号的差值。
如图1所示，3种重建算法基本上都能完成重建任务。但是在地震信号的第一个缺口处，FRLI重建结果显得比较“干净”;比较3种算法的重建误差图1(d),(f),(h)，发现FRSI和FRLI的重建误差比FRMN要小一些，3种重建效果的具体差异已通过量化比较。表2列出了叠后复杂地震信号的重建误差以及重建时间。
从表2可看出:FRLI重建时间为16.8 s耗时最少，比FRMN节省15.4 s，即耗时减少47.83%，这意味着FRLI的重建时间比FRMN少了将近一半，明显提高了重建速度。3种重建算法中，FRSI的重建误差极值和相对重建误差极值最小，FRLI重建均方误差相对FRSI稍大。上述结果表明:FRLI重建效果介于FRMN和FRSI之间，但优势在于计算简便，大幅减少了重建时间。
图1 叠后复杂地震信号及重建结果
Fig.1 Post-stack complex seismic signal and reconstruction results
表1 3种重建算法所选取的相关参数
Table 1 Selected parameters of the three reconstruction algorithms
表2 叠后复杂地震信号重建误差及重建时间
Table 2 Post-stack complex seismic signal reconstruction error and reconstruction time
现选取三维地震信号进行重建，地震信号共有55×55=3 025道，单道含128个采样点。图2为该地震信号。地震信号采取横纵切片分别重建，横切片选取60 ms处切片，纵切片为14和15号切片。为了使得重建后的地震信号效果更优，因此对比了几种参数的重建效果，表3列出了不同频带参数对应的重建误差以及重建时间。
从表3看出，频域带宽参数Mkx和Mky选择10，重建效果好。3种算法的相关参数选择如下:FRMN的频域带宽参数Mkx为10，Mky为10，预期信噪比为F=0.1;FRSI的频域带宽参数Mkx为10，Mky为10，预期信噪比为F=0.15，比例为P2=0.5，稀疏度a=1，最大迭代次数为L=6;FRLI的频域带宽参数Mkx为10，Mky为10，最大迭代次数为L=50，迭代截止条件δ1=10-4，δ2=10-5。
由图3可看出，3种算法都能还原地震信号，具体效果还要通过重建误差等指标量化分析。
图2 三维地震信号
Fig.2 3D seismic signal
表3 不同频带参数对应的重建误差及重建时间
Table 3 Reconstruction error and reconstruction time for different frequency band parameters
图3 合成三维地震信号及重建结果
Fig.3 Synthesis of 3D seismic signal and reconstruction results
表4列出了三维地震信号的重建误差以及重建时间。
表4 三维地震信号的重建误差及重建时间
Table 4 Reconstruction error and Reconstruction Time of 3D Seismic Signal
从表4可看出:FRLI对三维地震信号重建效果优于FRMN和FRSI，在重建时间上FRLI也是最快的，有效的缩短了重建时间，降低了均方误差，在三维地震信号重建中具有优势。
(1)二维地震信号的FRMN、FRSI和FRLI重建可知:FRSI的重建精度最高，但重建时间最长;FRLI的重建精度介于FRMN和FRSI之间，但由于其计算复杂度低，重建时间最短。因此，采用FRLI重建二维地震信号具有优势。
(2)三维地震信号的FRMN、FRSI和FRLI重建中，在重建精度方面FRLI最高，FRSI次之，FRMN重建精度最低;在重建速度上FRLI最快，FRMN次之，FRSI最慢。FRLI在不损耗较长时间的前提下，有效地提高了重建精度。
参考文献(References) :
[1] 金朝娣,能昌信,刘玉强,等.小波域相关算法在随钻地震信号处理中的应用[J].煤炭学报,2012,37(4):621-626.
JIN Chaodi,NENG Changxin,LIU Yuqiang,et al.Application of wavelet domain correlation algorithm in seismic signal processing while drilling[J].Journal of China Coal Society,2012,37(4):621-626.
[2] 周亚同,滕琳琳,李玲玲.基于高阶扩展快速行进法的缺失地震数据重建[J].石油地球物理勘探,2015,50(5):873-880，803.
ZHOU Yatong,TENG Linlin,LI Lingling.Reconstruction of missing seismic data based on high-order extended rapid marching method[J].Petroleum Geophysical Exploration,2015,50(5):873-880，803.
[3] 周亚同,王丽莉,蒲青山.压缩感知框架下基于K-奇异值分解字典学习的地震数据重建[J].石油地球物理勘探,2014,49(4):652-660，2.
ZHOU Yatong,WANG Lili,PU Qingshan.Reconstruction of seismic data based on K-singular value decomposition dictionary learning in compressed sensing[J].Petroleum Geophysical Exploration,2014,49(4):652-660，2.
[4] 周亚同,刘志峰,张志伟.形态分量分析框架下基于DCT与曲波字典组合的地震信号重建[J].石油物探,2015,54(5):560-568,581.
ZHOU Yatong,LIU Zhifeng,ZHANG Zhiwei.Reconstruction of seismic signals based on combination of DCT and Curvelet dictionary under morphological component analysis framework[J].Petroleum Geophysical Prospecting,2015,54(5):560-568,581.
[5] DUIJNDAM A J W,SCHONEWILLE M A,HINDRIKS C O H.Reconstruction of band-limited signals,irregularly sampled along one spatial direction[J].Geophysics,1999,64(2):524-538.
[6] HINDRIKS K,DUIJNDAM A J W.Reconstruction of 3D seismic signals irregularly sampled along two spatial coordinates[J].Geophysics,2000,65:253-263.
[7] ZWARTJES P M,GISOLF A.Fourier reconstruction with sparse inversion[J].Geophysical Prospecting,2007,55:199-221.
[8] ZWARTJES P M,SACCHI M D.Fourier reconstruction of non-uniformly sampled,aliased seismic data[J].Geophysics,2007,72(1):21-32.
[9] 孟小红,刘国峰,周建军.大间距地震数据重建方法研究[J].地球物理学进展,2006,21(3):687-691.
MENG Xiaohong,LIU Guofeng,ZHOU Jianjun.Research on reconstruction methods of large spacing seismic data[J].Progress in Geophysics,2006,21(3):687-691.
[10] 孟小红,郭良辉,张致付,等.基于非均匀快速傅里叶变换的最小二乘反演地震数据重建[J].地球物理学报,2008,51(1):235-241.
MENG Xiaohong,GUO Lianghui,ZHANG Zhifu,et al.Reconstruction of seismic data based on least squares inversion of non-uniform fast Fourier transform[J].Chinese Journal of Geophysics,2008,51(1):235-241.
[11] 周亚同,穆向阳,姜焕.大间距地震记录的插值辅助稀疏反演傅里叶重建算法[J].西安石油大学学报(自然科学版),2013,28(4):17-22,6.
ZHOU Yatong,MU Xiangyang,JIANG Huan.Interpolation-assisted Fourier reconstruction sparse inversion algorithm for large interval seismic records[J].Journal of Xi’an Shiyou University(Natural Science),2013,28(4):17-22,6.
[12] 周亚同,张伟,杨瑞霞.一维非均匀采样信号可变稀疏度傅里叶重建算法研究[J].微电子学与计算机,2012,29(7):188-191.
ZHOU Yatong,ZHANG Wei,YANG Ruixia.Research on variable sparsity Fourier reconstruction algorithm for one-dimensional non-uniform sampling signals[J].Microelectronics & Computer,2012,29(7):188-191.
[13] 王化祥,王超,陈磊.基于Landweber迭代的图像重建算法[J].信号处理,2000,16(4):354-356,328.
WANG Huaxiang,WANG Chao,CHEN Lei.Image reconstruction algorithm based on Landweber iteration[J].Signal Processing,2000,16(4):354-356,328.
[14] 王旭,王静文,王柯元.阈值Landweber在MIT图像重建中的应用[J].东北大学学报(自然科学版),2016,37(4):477-480.
WANG Xu,WANG Jingwen,WANG Keyuan.Application of Threshold Landweber in MIT image reconstruction[J].Journal of Northeastern University(Natural Science),2016,37(4):477-480.
[15] 张伟斌,冯象初,王卫卫.图像恢复的小波域加速Landweber迭代阈值方法[J].电子与信息学报,2011,33(2):342-346.
ZHANG Weibin,FENG Xiangchu,WANG Weiwei.An algorithm for accelerating Landweber iteration threshold in wavelet domain for image restoration[J].Journal of Electronics and Information Technology,2011,33(2):342-346.
[16] 肖理庆,王化祥,厉丹.改进Landweber预迭代ERT图像重建算法[J].中国电机工程学报,2013,33(23):118-125,1.
XIAO Liqing,WANG Huaxiang,LI Dan.Improved Landweber pre-iteration ERT image reconstruction algorithm[J].Proceedings of the CSEE,2013,33(23):118-125,1.
[17] 张立峰,王化祥.一种修正的电阻层析成像Landweber迭代算法[J].计量学报,2016,37(3):271-274.
ZHANG Lifeng,WANG Huaxiang.A modified Landweber iterative algorithm for Electrical Resistance Tomography[J].Acta Metrologica Sinica,2016,37(3):271-274.
[18] 董向元,陈琪,李惊涛.基于快速投影Landweber法的电容层析成像图像重建算法研究[J].中国电机工程学报,2005,25(14):89-92.
DONG Xiangyuan,CHEN Qi,LI Jingtao.Research on image reconstruction algorithm of electrical capacitance tomography based on rapid projection Landweber method[J].Proceedings of the CSEE,2005,25(14):89-92.
[19] 张春利.基于Landweber迭代格式的带限信号非均匀采样重建算法研究[D].北京:北京交通大学,2014.
ZHANG Chunli.Research on non-uniform sampling reconstruction algorithm for band-limited signals based on Landweber iterative format[D].Beijing:Beijing Jiaotong University,2014.
[20] 钟瑞淇.基于FFT-Landweber迭代的带限函数外推算法研究[D].北京:北京交通大学,2015.
ZHONG Ruiqi.Research on extrapolation algorithm with band-limit function based on FFT-Landweber iteration[D].Beijing:Beijing Jiaotong University,2015.
[21] 杨真真,杨震,李雷.语音重构的DCT域加速Landweber迭代硬阈值算法[J].信号处理,2012,28(2):172-178.
YANG Zhenzhen,YANG Zhen,LI Lei.Accelerated Landweber iterative hard threshold algorithm in DCT domain for speech reconstruction[J].Signal Processing,2012,28(2):172-178.
[22] ZENG G L.A filtered backprojection algorithm with characteristics of the iterative Landweber algorithm[J].Medical Physics,2012,39(2):603-607.
[23] HEIN T,KAZIMIERSKI K S.Accelerated Landweber iteration in banach spaces[J].Inverse Problems,2010,26(26):1037-1050.
[24] BRIANZI P,BENEDETTO F D,ESTATICO C.Improvement of space-invariant image deblurring by preconditioned Landweber iterations[J].Siam Journal on Scientific Computing,2008,30(3):1430-1458.
[25] LANDWEBER L.An iteration formula for fredholm integral equations of the first kind[J].American Journal of Mathematics,1951,73(3):615-624.
[26] 张军.求解不适定问题的快速Landweber迭代法[J].数学杂志,2005,25(3):333-335.
ZHANG Jun.Fast Landweber iterative method for solving ill-posed problems[J].Mathematic Journal,2005,25(3):333-335.