CN104570141B - 一种基于双谱运算的重力异常分离方法 - Google Patents
一种基于双谱运算的重力异常分离方法 Download PDFInfo
- Publication number
- CN104570141B CN104570141B CN201310491933.0A CN201310491933A CN104570141B CN 104570141 B CN104570141 B CN 104570141B CN 201310491933 A CN201310491933 A CN 201310491933A CN 104570141 B CN104570141 B CN 104570141B
- Authority
- CN
- China
- Prior art keywords
- signal
- gravity anomaly
- bispectrum
- reconstructed
- data
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Abstract
本发明提供了一种基于双谱运算的重力异常分离方法，属于应用地球物理领域。本方法在信号重构中，采用一个控制参数，将该控制参数与原完整信号的谐波项数M相乘，得到新的谐波项数M′，利用该新的谐波项数M′对原完整信号进行重构，获得重构重力异常信号，实现重力异常分离；所述控制参数为0.0～1.0之间。本发明用于实际重力异常的分离处理，使线性条带状分布的地质特征得到清楚、正确的显示；本发明为信号恢复、重构处理增加了一个有用的工具。
Description
技术领域
本发明属于应用地球物理领域，具体涉及一种基于双谱运算的重力异常分离方法，同时适用于通讯工程、电子工程、军事、医学等多个应用领域。
背景技术
在石油地球物理勘探等应用中，所获得的重力异常，包含了从地壳深部到地表的所有密度不均匀地质体的影响。不同地质因素引起的重力异常，叠加在一起，给人们在识别、区分和研究上带来相当大的困难，因此针对石油勘探不同阶段不同地质任务的需要，要进行重力异常的分离处理，以便能够对重力异常做出更准确的解释。常规的重力异常分离方法，包括图解法、平均场法、高次导数法、上下延拓法、函数逼近法(也称趋势分析法)以及频率域滤波法等。近年对这些传统的常规分离方法通过引入新的计算技术，延伸出了数值切割法、小波变换法、分数维法等。
现有的这些重力异常分离法，有一个共同的特点，就是呈线性条带状分布的地质构造在应用现有的这些技术进行重力异常分离处理后，特征显示不够明显，给识别与解释带来一定的困难。
发明内容
本发明的目的在于解决上述现有技术中存在的难题，提供一种基于双谱运算的重力异常分离方法，通过引入基于双谱运算的信号重构技术，增强重力异常分离的功能，实现对呈线性条带状分布的地质构造的正确识别与解释。
本发明是通过以下技术方案实现的：
一种基于双谱运算的重力异常分离方法，所述方法在信号重构中，采用一个控制参数，将该控制参数与原完整信号的谐波项数M相乘，得到新的谐波项数M′，利用该新的谐波项数M′对原完整信号进行重构，获得重构重力异常信号，实现重力异常分离；
所述控制参数为0.0～1.0之间。
所述方法包括：
第一步，读取原始重力异常信号x(i),i＝1，…，n；其中，n是采样点数；
第二步，对原始重力异常信号x(i)，i＝1，…，n进行预处理，得到新信号xn(i)，i＝1，…，kn；
第三步，通过傅里叶变换，计算得到新信号xn(i)，i=1，…，kn的频谱X(f)；
第四步，计算信号xn(i)，i=1，…，kn的双谱B(f1，f2)；
第五步，利用对称性计算信号xn(i)，i=1，…，kn在整个第一象限的双谱B(f1，f2)；
第六步，取双谱B(f1，f2)在f1＝f2时的模|B(f，f)|的第一个峰值时的f为估算的基频f0；
第七步，在原信号xn(i)，i=1，…，kn上加一相位为0的半基频余弦信号得到信号y(i)，i=1，…，kn，y(i)=xn(i)+cos(πf0i)，i＝1，…，kn；
第八步，通过傅里叶变换，计算信号y(i)，i=1，…，kn的频谱Y(f)；
第九步，计算信号y(i)，i=1，…，kn的双谱By(f1，f2)；
第十步，利用对称性计算信号y(i)，i=1，…，kn在整个第一象限的双谱By(f1，f2)；
第十一步，计算谐波项系数ai与相位角及谐波项数M；
第十二步，输入控制参数q，将控制参数q与谐波项数M相乘，得重构谐波项数M′；
第十三步，利用谐波项系数ai、相位角与重构谐波项数M′对信号进行重构得到得重构信号数据；
第十四步，对重构信号数据进行线性背景补偿和均值补偿处理得到重构重力异常；
第十五步，将第一步读取的原始重力异常信号x(i)，i＝1，…，n与第十四步得到的重构重力异常对应相减，得剩余重力异常；
第十六步，输出第十五步计算得到的剩余重力异常。
所述第二步的预处理包括去线性背景、零均值化和周期拓展，具体如下：
(1)去线性背景：
利用以下公式计算信号x(i)，i＝1，…，n的线性背景斜率bg：
式中，x(n)为原始重力异常信号的最末一个数据，x(1)为原始重力异常信号的第一个数据，n为原始重力异常信号的采样点数；
然后利用以下公式消去信号的线性背景：
xb(i)=x(i)-x(1)-bg×(i-1)
式中，xb(i)，i=1，…，n为原始重力异常信号消去线性背景后的数据；
(2)零均值化：
对消去线性背景后的数据xb(i)，i=1，…，n的各个样点数据xb(i)，求和并除以样点数n，得消去线性背景后的数据xb(i)，i=1，…，n的平均值
消去线性背景后的数据xb(i)，i=1，…，n减去平均值实现零均值化：
式中，xe(i)，i=1，…，n为原始重力异常信号去线性背景后又零均值化后的信号；
(3)周期拓展：
原始重力异常信号去线性背景后又零均值化后的信号xe(i)，i=1，…，n是周期数为1的信号，将K个该信号首尾相连得到周期数为K的新信号xn(i)，i＝1，…，kn；
所述第三步是通过下式实现的：
X(f)=A(f)·ejφ(f)
其中，A(f)和φ(f)分别为频谱X(f)的幅值和相位；其中，频谱X(f)的长度为k·n。
所述第四步和第九步都是通过公式(4)计算实现的：
B(f1，f2)=X(f1)·X(f2)·X*(f1+f2)
=A(f1)·A(f2)·A(f1+f2)·exp{j[φ(f1)+φ(f2)-φ(f1+f2)]}
=A(f1，f2)·exp[jφ(f1，f2)] (4)
其中，X(f1)，X(f2)分别表示变量f取f1，f2时对应得到的频谱X(f)；X*(f1+f2)是频谱X(f1+f2)的共轭，其中f1，f2是与时间变量t1，t2相对应的频率变量。
所述第十一步中是根据式(14)计算得到谐波项系数ai与相位角
所述谐波项数M是由下式求得：
M=km/nf
其中，km为信号xn(i)，i=1，…，kn的频谱X(f)的长度的1/2，即km=k·n/2+1；nf是估算得到的基频f0所对应的频率采样序号。
所述第十三步是这样实现的：
将第十一步得到的谐波项系数ai、相位角和第十二步得到的重构谐波项数M′代入式(5)，得重构信号数据：
与现有技术相比，本发明的有益效果是：本发明使表征断层构造的台阶模型的重力异常，得到了准确的分离。同时本发明用于实际重力异常的分离处理，使线性条带状分布的地质特征得到清楚、正确的显示；本发明为信号恢复、重构处理增加了一个有用的工具。
附图说明
图1为本发明方法的步骤框图；
图2为某试验区的实际布格重力异常图；
图3为实际布格重力异常应用本发明研制的基于双谱运算的异常分离新方法得到的重构重力异常；
图4为实际布格重力异常减去重构重力异常得到的剩余重力异常。
具体实施方式
下面结合附图对本发明作进一步详细描述：
(1)基于双谱运算的重力异常分离方法原理
设x(t)是一均值为0的连续实信号，其Fourier变换为：
X(f)=A(f)·exp[jφ(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)的双谱表达式：
B(f1，f2)=X(f1)·X(f2)·X*(f1+f2)
=A(f1)·A(f2)·A(f1+f2)·exp{j[φ(f1)+φ(f2)-φ(f1+f2)]}
=A(f1，f2)·exp[jφ(f1，f2)] (4)
设重力异常数据为x(n)是连续实信号x(t)的离散形式，可以用谐波信号表达式来描述，即
式中：M为谐波项数，ai和分别为谐波系数和相位角，f0是基频率。
根据Fourier变换的特性，由式(4)和(5)可以得到x(n)的双谱振幅谱、相位谱与谐波系数、相位角的关系式：
方程式(6a)、(6b)中的未知数也即式(5)中的未知数，有基频f0、谐波系数ai与相位角显然，(6a)、(6b)中未知数的个数大于方程数。实际应用时须先对基频f0进行准确地估算，这是本方法成功实现的关键，而且，由式(5)可知，基频f0估算的准确度，决定了重构信号的保真度(下面第二步中的周期拓展，使信号的基频f0与频率分辨率Δf相等。这样，不管原始信号的信噪比如何，基频f0与频率分辨率Δf的比值始终是整数1，从而，确保了估算得到的基频f0是真值，满足了对基频准确度的要求。这一内容是专利号：201010508001.9的专利的核心内容)。
根据估算得到的基频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)
其中，X(f)为x(t)的频谱，Y(f)为y(t)的频谱，Z(f)为z(t)的频谱。
故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(n)中的谐波系数ai与相位角求解的递推公式：
将式(14)所求得的谐波系数ai和相位角代入式(5)，即可实现对信号x(n)的重构(公式(14)中的A(f0/2，f0/2)，A[f0,9i-1)f0]，φy(f0/2，f0/2)，φ(f0，mf0)等是前面求得的双谱的振幅谱和相位谱)。一个完整的信号x(n)是由M个谐波项组成的，若只取其中某些谐波项进行信号重构，则可以获取某些频率成分的信号，对重力异常数据而言，就可以实现重力异常的分离。
(2)本发明的内容
在应用双谱运算进行重力异常信号重构计算的过程中，在对原始信号做双谱运算前，先对原始信号做去线性背景、零均值化、周期拓展等3项预处理。
本发明的核心内容是，在对预处理后的重力异常信号进行双谱运算后，最后进行信号重构计算时，输入一个控制参数(参数的取值范围为：0.0～1.0之间)，将该参数与谐波项数M相乘，得到一个新的谐波项数M′，将新的谐波项数M′与所求得的谐波系数ai和相位角一起代入式(5)，即得到重构重力异常信号(一个完整的信号x(n)是由M个谐波项组成的，即取控制参数值为1.0时，将使原来的信号得到完全恢复，这是专利号为201010508001.9的专利实现的效果，本发明是取控制参数小于1.0，即只取其中一部分谐波项进行信号重构，实现对原始信号的部分恢复)，将原始重力异常信号减去重构重力异常信号得剩余重力异常信号，从而实现重力异常的分离。
技术流程详细说明：
第一步，读取野外数据采集成果中布格重力异常x(i)，i＝1，…，n的数据文件，其中，n是采样点数。
第二步，对布格重力异常数据x(i)，i＝1，…，n进行去线性背景、零均值化、周期拓展等预处理，得到xn(i)，i＝1，…，kn。预处理的详细过程请参考发明专利“一种信号恢复的优化方法”，专利号：201010508001.9，下面仅给出用到的公式：
(1)去线性背景过程：
利用以下公式计算信号x(i)，i＝1，…，n的线性背景斜率bg：
式中，x(n)为布格重力异常的最末一个数据，x(1)为布格重力异常的第一个数据，n为布格重力异常的采样点数；
然后利用以下公式消去信号的线性背景：
xb(i)=x(i)-x(1)-bg×(i-1)
式中，xb(i)，i=1，…，n为布格重力异常消去线性背景后的数据，x(i)，i=1，…，n为原始布格重力异常的数据，x(1)为原始布格重力异常的第一个数据；
(2)零均值化过程：
对信号x(i)，i＝1，…，n消去线性背景后的数据xb(i)，i=1，…，n的各个样点数据xb(i)，求和并除以样点数n，得信号去线性背景后数据xb(i)，i=1，…，n的平均值
式中，n为布格重力异常的采样点数；
去线性背景后的数据xb(i)，i=1，…，n减去平均值实现零均值化：
式中，xe(i)，i=1，…，n为信号布格重力异常去线性背景后又零均值化后的数据；
(3)周期拓展过程：
原信号是周期数为1的信号，将原信号首尾相连重复放置于后，重复放置K-1次，即得周期数为K的新信号xn(i)，i=1，…，kn；
周期拓展过程的C语言计算程序语句如下：
for(i=0；i<k；i++)for(j=0；j<n；j++)xn[j+i*k]=xe[j]；
其中，xe[j]为原始布格重力异常信号x(i)，i＝1，…，n去线性背景后又零均值化后数据xe(i)，i=1，…，n的数组，xn[j+i*k]为周期拓展后信号xn(i)，i=1，…，kn的数组，n为信号x(i)，i=1，…，n的采样点数，k为周期拓展后信号的周期数；
第三步，通过Fourier变换，计算去线性背景后又零均值化周期拓展后布格重力异常数据xn(i)，i＝1，…，kn的频谱X(f)。
X(f)=A(f)·ejφ(f）
其中，A(f)和φ(f)分别为频谱X(f)的幅值和相位；其中，频谱X(f)的长度为k·n；其中，n为信号布格重力异常信号x(i)，i=1，…，n的采样点数，k为周期拓展后信号的周期数；
第四步，根据式(4)，通过三重相关计算数据xn(i)，i=1，…，kn的双谱B(f1，f2)。
这一过程的C语言计算程序语句如下：
其中，数组a[i]、f[i]分别是信号xn(i)，i＝1，…，kn的频谱X(f)的振幅谱和相位谱数组，数组ba[i][j]、bf[i][j]分别是信号xn(i)，i=1，…，kn的双谱B(f1，f2)的振幅谱A(f1，f2)和相位谱φ(f1，f2)数组；其中km信号xn(i)，i＝1，…，kn的频谱X(f)的长度的1/2，即km＝k·n/2+1；
第五步，由于第四步得到的双谱数据，仅仅在第一象限中只占下半空间，根据双谱性质中的对称性，在第一象限中B(f2，f1)=B(f1，f2)，计算信号xn(i)，i=1，…，kn在整个第一象限的双谱B(f1，f2)；
这一过程的C语言计算程序语句如下：
其中，数组ba[i][j]、bf[i][j]分别是信号xn(i)，i＝1，…，kn的双谱B(f1，f2)的振幅谱数组和相位谱数组；其中km信号xn(i)，i=1，…，kn的频谱X(f)的长度的1/2，即km=k·n/2+1；
第六步，取双谱B(f1，f2)f1＝f2时的模|B(f，f)|(B(f1，f2)表示一个函数有2个变量f1，f2，当2个变量相等时，即f1=f2时，也即相当于坐标系中y=x时，只用一个通用的字母f表示，于是得到B(f，f))的第一个峰值时的f为估算的基频f0。
这一过程的C语言计算程序语句如下：
第七步，在原数据信号xn(i)，i=1，…，kn上加一相位为0的半基频余弦信号，即y(i)=xn(i)+cos(πf0i)，i=1，…，kn。
第八步，通过Fourier变换，计算数据信号y(i)，i=1，…，kn的频谱Y(f)。
第九步，根据式(4)，通过三重相关计算数据信号y(i)，i=1，…，kn的双谱By(f1，f2)。
第十步，利用对称性计算数据信号y(i)，i=1，…，kn整个第一象限的双谱By(f1，f2)。
第十一步，根据式(14)计算谐波项系数ai与相位角及谐波项数M。
谐波项数M，可由下式求得：
M=km/nf
其中，km为信号xn(i)，i＝1，…，kn的频谱X(f)的长度的1/2，即km=k·n/2+1；nf是估算得到的基频f0所对应的频率采样序号；n为信号x(i)，i=1，…，n的采样点数，k为周期拓展后信号的周期数；
第十二步，输入控制参数q(控制参数的数值是人为确定的，一般在计算时，按照0.1，0.2，0.3，…，逐一输入，分别计算得到相应的结果，根据经验和其它资料来判断合适的结果，相应的参数即是合适的参数)，将控制参数q与谐波项数M相乘，得重构谐波项数M′。
第十三步，谐波项系数ai、相位角与重构谐波项数M′代入式(5)，得重构信号数据。
第十四步，对重构信号数据进行线性背景补偿和均值补偿处理得重构重力异常。
所述均值补偿处理过程为：
对输出信号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为原始信号的线性背景斜率。
第十五步，将第一步读取的布格重力异常x(i)，i＝1，…，n与重构得到的重力异常数据对应相减，得剩余重力异常。
第十六步，输出计算得到的剩余重力异常结果。
根据上述技术流程，发明人在Windows操作系统下利用C++Builder，依照图1的程序框图，实现了采用本发明所述的方法做的重力异常分离处理，图中虚框部分为本发明的核心内容。
对某试验工区的实际布格重力异常数据(如图2所示)，应用本发明方法得到的重构重力异常如图3所示，较好地反映了该地区深部地质构造特征；由布格重力异常减去重构重力异常得到的剩余重力异常，由图4可见，一些北北东、北北西及近南北向的线性条带状展布的沉积盖层内地质构造特征得到了清楚、正确的显示，这一特征与该工区的实际地质特征相符。
本发明在应用双谱运算实现信号有效恢复的基础上，在最后进行信号重构计算时，通过输入一个控制参数，取其中部分谐波项进行信号恢复重构计算，得到重构重力异常，将原始重力异常减去重构重力异常得剩余重力异常，从而实现重力异常的分离。
本发明使表征断层构造的台阶模型的重力异常，得到了准确的分离。本发明用于实际重力异常的分离处理，使线性条带状分布的地质特征得到清楚、正确的显示。
上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。
Claims (6)
1.一种基于双谱运算的重力异常分离方法，其特征在于：所述方法在信号重构中，采用一个控制参数，将该控制参数与原完整信号的谐波项数M相乘，得到重构谐波项数M′,利用该重构谐波项数M′对原完整信号进行重构，获得重构重力异常信号，实现重力异常分离；
所述控制参数为0.0～1.0之间；
所述方法包括：
第一步，读取原始重力异常信号x(i),i＝1,…,n；其中，n是采样点数；
第二步，对原始重力异常信号x(i),i＝1,…,n进行预处理，得到新信号xn(i),i＝1,…,kn，k为周期数；
第三步，通过傅里叶变换，计算得到新信号xn(i),i＝1,…,kn的频谱X(f)；
第四步，计算新信号xn(i),i＝1,…,kn的双谱B(f1，f2)；
第五步，利用对称性计算新信号xn(i),i＝1,…,kn在整个第一象限的双谱B(f1,f2)'；
第六步，取双谱B(f1，f2)'在f1＝f2时的模|B(f,f)'|的第一个峰值时的f为估算的基频f0；
第七步，在新信号xn(i),i＝1,…,kn上加一相位为0的半基频余弦信号得到信号y(i),i＝1,…,kn，y(i)＝xn(i)+cos(πf0i),i＝1,…,kn；
第八步，通过傅里叶变换，计算信号y(i),i＝1,…,kn的频谱Y(f)；
第九步，计算信号y(i),i＝1,…,kn的双谱By(f1，f2)；
第十步，利用对称性计算信号y(i),i＝1,…,kn在整个第一象限的双谱By(f1,f2)'；
第十一步，计算谐波项系数ai与相位角及谐波项数M；
第十二步，输入控制参数q，将控制参数q与谐波项数M相乘，得到重构谐波项数M′；
第十三步，利用谐波项系数ai、相位角与重构谐波项数M′对原始重力异常信号进行重构得到重构信号数据；
第十四步，对重构信号数据进行线性背景补偿和均值补偿处理得到重构重力异常信号；
第十五步，将第一步读取的原始重力异常信号x(i),i＝1,…,n与第十四步得到的重构重力异常信号对应相减，得到剩余重力异常；
第十六步，输出第十五步计算得到的剩余重力异常。
2.根据权利要求1所述的基于双谱运算的重力异常分离方法，其特征在于：所述第二步的预处理包括去线性背景、零均值化和周期拓展，具体如下：
(1)去线性背景：
利用以下公式计算原始重力异常信号x(i),i＝1,…,n的线性背景斜率bg：
<mrow>
<mi>b</mi>
<mi>g</mi>
<mo>=</mo>
<mfrac>
<mrow>
<mi>x</mi>
<mrow>
<mo>(</mo>
<mi>n</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mi>x</mi>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</mfrac>
</mrow>
式中，x(n)为原始重力异常信号的最末一个数据，x(1)为原始重力异常信号的第一个数据，n为原始重力异常信号的采样点数；
然后利用以下公式消去信号的线性背景：
xb(i)＝x(i)-x(1)-bg×(i-1)
式中，xb(i),i＝1,…,n为原始重力异常信号消去线性背景后的数据；
(2)零均值化：
对消去线性背景后的数据xb(i),i＝1,…,n的各个样点数据xb(i)，求和并除以样点数n，得消去线性背景后的数据xb(i),i＝1,…,n的平均值
<mrow>
<mover>
<mrow>
<mi>x</mi>
<mi>b</mi>
</mrow>
<mo>&OverBar;</mo>
</mover>
<mo>=</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>n</mi>
</munderover>
<mi>x</mi>
<mi>b</mi>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>)</mo>
</mrow>
<mo>/</mo>
<mi>n</mi>
</mrow>
消去线性背景后的数据xb(i),i＝1,…,n减去平均值实现零均值化：
<mrow>
<mi>x</mi>
<mi>e</mi>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mi>x</mi>
<mi>b</mi>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mover>
<mrow>
<mi>x</mi>
<mi>b</mi>
</mrow>
<mo>&OverBar;</mo>
</mover>
</mrow>
式中，xe(i),i＝1,…,n为原始重力异常信号去线性背景后又零均值化后的信号；
(3)周期拓展：
原始重力异常信号去线性背景后又零均值化后的信号xe(i),i＝1,…,n是周期数为1的信号，将k个该信号首尾相连得到周期数为k的新信号xn(i),i＝1,…,kn。
3.根据权利要求2所述的基于双谱运算的重力异常分离方法，其特征在于：所述第三步是通过下式实现的：
X(f)＝A(f)·ejφ(f)
其中，A(f)和φ(f)分别为频谱X(f)的幅值和相位；其中，频谱X(f)的长度为k·n。
4.根据权利要求3所述的基于双谱运算的重力异常分离方法，其特征在于：所述第四步和第九步都是通过公式(4)计算实现的：
B(f1,f2)＝X(f1)·X(f2)·X*(f1+f2)
＝A(f1)·A(f2)·A(f1+f2)·exp{j[φ(f1)+φ(f2)-φ(f1+f2)]}
＝A(f1,f2)·exp[jφ(f1,f2)] (4)
其中，X(f1),X(f2)分别表示变量f取f1,f2时对应得到的频谱X(f)；X*(f1+f2)是频谱X(f1+f2)的共轭，其中f1,f2是与时间变量t1,t2相对应的频率变量。
5.根据权利要求4所述的基于双谱运算的重力异常分离方法，其特征在于：所述第十一步中是根据式(14a)和式(14b)计算得到谐波项系数ai与相位角
<mrow>
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<msub>
<mi>a</mi>
<mn>1</mn>
</msub>
<mo>=</mo>
<mn>8</mn>
<mi>A</mi>
<mrow>
<mo>(</mo>
<mfrac>
<msub>
<mi>f</mi>
<mn>0</mn>
</msub>
<mn>2</mn>
</mfrac>
<mo>,</mo>
<mfrac>
<msub>
<mi>f</mi>
<mn>0</mn>
</msub>
<mn>2</mn>
</mfrac>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>a</mi>
<mi>i</mi>
</msub>
<mo>=</mo>
<mfrac>
<mrow>
<mn>8</mn>
<mi>A</mi>
<mo>&lsqb;</mo>
<msub>
<mi>f</mi>
<mn>0</mn>
</msub>
<mo>,</mo>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>-</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<msub>
<mi>f</mi>
<mn>0</mn>
</msub>
<mo>&rsqb;</mo>
</mrow>
<mrow>
<msub>
<mi>a</mi>
<mn>1</mn>
</msub>
<msub>
<mi>a</mi>
<mrow>
<mi>i</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
</mrow>
</mfrac>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>,</mo>
<mi>i</mi>
<mo>=</mo>
<mn>2</mn>
<mo>,</mo>
<mn>3</mn>
<mo>,</mo>
<mo>...</mo>
<mo>,</mo>
<mi>M</mi>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>14</mn>
<mi>a</mi>
<mo>)</mo>
</mrow>
</mrow>
所述谐波项数M是由下式求得：
M＝km/nf
其中，km为信号xn(i),i＝1,…,kn的频谱X(f)的长度的1/2，即km＝k·n/2+1；nf是估算得到的基频f0所对应的频率采样序号。
6.根据权利要求5所述的基于双谱运算的重力异常分离方法，其特征在于：所述第十三步是这样实现的：
将第十一步得到的谐波项系数ai、相位角和第十二步得到的重构谐波项数M′代入式(5)，得到重构信号数据：
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201310491933.0A CN104570141B (zh) | 2013-10-18 | 2013-10-18 | 一种基于双谱运算的重力异常分离方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201310491933.0A CN104570141B (zh) | 2013-10-18 | 2013-10-18 | 一种基于双谱运算的重力异常分离方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN104570141A CN104570141A (zh) | 2015-04-29 |
| CN104570141B true CN104570141B (zh) | 2018-01-16 |
Family
ID=53086655
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201310491933.0A Active CN104570141B (zh) | 2013-10-18 | 2013-10-18 | 一种基于双谱运算的重力异常分离方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104570141B (zh) |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101470194B (zh) * | 2007-12-26 | 2012-07-04 | 中国科学院声学研究所 | 一种水雷目标的识别方法 |
| US10575751B2 (en) * | 2008-11-28 | 2020-03-03 | The University Of Queensland | Method and apparatus for determining sleep states |
| CN101997788B (zh) * | 2010-10-15 | 2013-07-31 | 中国石油化工股份有限公司 | 一种信号恢复的优化方法 |
-
2013
- 2013-10-18 CN CN201310491933.0A patent/CN104570141B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN104570141A (zh) | 2015-04-29 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Dai et al. | Multipath mitigation via component analysis methods for GPS dynamic deformation monitoring | |
| Pinnegar et al. | The S-transform with windows of arbitrary and varying shape | |
| Wang | Multichannel matching pursuit for seismic trace decomposition | |
| CN101930598B (zh) | 基于shearlet域非局部均值的自然图像去噪方法 | |
| EP2718746A2 (en) | System and method for seismic data inversion by non-linear model update | |
| CN105549097A (zh) | 一种瞬变电磁信号工频及其谐波干扰消除方法及装置 | |
| CN104360393A (zh) | 一种地震数据重建方法 | |
| Li et al. | Audio magnetotelluric signal-noise identification and separation based on multifractal spectrum and matching pursuit | |
| Maharramov et al. | Robust joint full-waveform inversion of time-lapse seismic data sets with total-variation regularization | |
| Yue et al. | Suppression of periodic interference during tunnel seismic predictions via the Hankel-SVD-ICA method | |
| Wang et al. | Continuous wavelet analysis of matter clustering using the gaussian-derived wavelet | |
| Baziw et al. | Principle phase decomposition: A new concept in blind seismic deconvolution | |
| CN104570141B (zh) | 一种基于双谱运算的重力异常分离方法 | |
| Ghozzi et al. | Peak detection of GPR data with lifting wavelet transform (LWT) | |
| Li et al. | A hybrid filtering method based on a novel empirical mode decomposition for friction signals | |
| CN101997788B (zh) | 一种信号恢复的优化方法 | |
| Baziw | Implementation of the principle phase decomposition algorithm | |
| CN105319594A (zh) | 一种基于最小二乘参数反演的傅里叶域地震数据重构方法 | |
| Hu et al. | A non-iterative approach to inverse elastic scattering by unbounded rigid rough surfaces | |
| Fu et al. | A wavelet multiscale–homotopy method for the inverse problem of two-dimensional acoustic wave equation | |
| Ting et al. | Eeg signal processing based on proper orthogonal decomposition | |
| Zeng et al. | Corresponding relationships between nodes of decomposition tree of wavelet packet and frequency bands of signal subspace | |
| Fichtner et al. | Optimal spherical spline filters for the analysis and comparison of regional-scale tomographic models | |
| Saibaba et al. | A fast Kalman filter for time-lapse electrical resistivity tomography | |
| Kumar et al. | Phasor Estimation and Disturbance Detection in Power Systems by Ramanujan Sum |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |