CN1961221A - 处理电磁数据的Kalman滤波方法 - Google Patents

Abstract

一方法，用于使用Kalman滤波或其他用来跟踪信号幅度随时间的变化和检测所述信号相位随时间的改变的跟踪算法34，来跟踪在噪声地震数据中检测的正弦电磁信号。该方法用于受控的地震源电磁勘查，其中地震源相对于接收器的较远偏移可能引起所述地震源信号显著地衰减，并使其难于从地磁或其他电磁背景31中恢复。

Description

处理电磁数据的Kalman滤波方法

【0001】本申请要求2004年6月1日申请的美国临时专利申请No.60/576,201的权益。

技术领域

【0002】本发明一般涉及地球物理勘探领域，更具体地涉及用来探测碳氢化合物的电磁方法。具体来说，本发明是用于跟踪用于受控的源电磁勘探的电磁源信号，以从噪声恢复该信号的方法。

背景技术

【0003】受控源电磁(“CSEM”)地球物理勘查使用活性(人造)源来产生电磁场激励地球，并配置地球表面，海床上，或钻孔内的接收器仪器，以测量所产生的电场和磁场，即地球对源激励的响应。图1图解说明了海岸CSEM勘查的基本元件。船只在海床13上方拖曳浸入水中的CSEM发射器11。然后，分析由接收器12测量的电场和磁场，以确定表面或海床下面地球结构(地表下地层)的电阻率。该技术已经应用陆上矿物勘探，海洋构造研究，和海岸石油及矿物资源勘探。

【0004】活性电磁源信号可被当作正弦信号的和处理(如，由带有奇次谐波的基频组成的方波信号)。该源的例子是用于许多CSEM工作的水平电偶极子。随着偏移(即该偶极子源11和接收器12之间的距离)的增加，正弦信号可能显著衰减。而且，远偏移对于确定感兴趣的深阻结构来说经常是关键。结果，存在使该正弦信号获得优良信噪比的需求。

【0005】典型的改善该EM数据的信号噪声的处理方法包括将数据分解到多个时间窗口中，在这些时间窗口上使用傅立叶分析或相似的方法，来计算所选频率分量(一个或多个)的幅度和相位。例如，可见Constable和Cox的“Marine controlled-source electromagneticsounding 2.The PEGASUS Experiment，”Journal of Geophysicalresearch 101，5519-5530(1996)。这些窗口不能太长，这是因为信号幅度和相对相位可能在该分析窗口中显著改变。然而，短窗仅允许最小的信噪比改善。当前的方法要求在这两个极端之间平衡。

【0006】现有方法的另一个问题是他们不利用信号和噪声相关性。具体来说，对于活性源海洋EM成像，由于低频大地电磁(“MT”)噪声可伪装成信号，因此低频MT噪声就是一个特别重要的难题。(MT噪声是来自自然非活性源的电磁发射。)不同探测器之间的相关性可用来帮助从这些噪声中分离活性源信号。当前方法中，没有最优地使用其他信号和噪声的相关性(如，两个水平分量信号的相关性)。

【0007】Kalman(卡尔曼)滤波算法起源于导航定位问题并特别适用于跟踪类的问题(Kalman，1960)。最初由Kalman发表于《Trans.Ofthe ASME-J.OF Basic Engr.》35-45(1960)上，还发表了大量文献作为关于基本Kalman滤波的改进和应用的总结，例如在由Brown在纽约John Wiley & Sons(1983)上发表的Introduction to RandomSignal Analysis and Kalman Filtering中。这些改进中的几种对于当前发明的某些实施例是重要的。

【0008】标准Kalman滤波在一个方向上运行，并滤波该方向(或时间)序列中的数据。因此仅前面的数据影响滤波结果。由Rauch等人提出的重要修改给出最优处理，其使用整个时间记录：参看Rauch的“Solutions to the linear smoothing problem”，IEEE Trans.On Auto.Control，AC-8，371(1963)；和Rauch等人的“Maximum likelihoodestimates of linear dynamic systems”，AMA J. 3，1445(1965)。Szelag公开了另一个算法改进，其允许滤波器跟踪已知频率的正弦信号；参见“A short term forecasting algorithm for trunk demand servicing”，TheBell System Technical Journal 61，67-96(1982)。其被发展，用于以电话中继负载值(telephone trunk load values)来跟踪年度周期。

【0009】图2是流程图，其图解说明了Kalman算法。更详细的可参考Brown的论文集第200页。

【0010】La Scala等人公开了公知的用于跟踪时变频率的扩展kalman滤波器的使用。(“Design ofan extended Kalman filter frequencytracker”，IEEE Transaction on Signal Processing 44，No.3，739-742(1996年3月))所述公式假定信号幅度保持恒定。所用的特殊Kalman算法旨在跟踪未知频率的信号，其中频率可经历相当的改变。Lagunas等人公开了一种扩展的Kalman滤波器，用于跟踪有噪声和频率变化的复杂正弦，如Doppler移动。如同La Scala的方法，Lagunas的算法被设计来跟踪未知频率的正弦曲线或正弦波。因此，如果应用来跟踪具有恒定或近似恒定频率的信号，这两种方法都是次优的。如果改变相对小，则Lagunas的方法能够跟踪幅度变化。本发明不是旨在处理使用发射已知频率波形的电磁源的电磁勘探数据。需要使用大窗口或甚至所有电磁数据，来跟踪关于已知正弦曲线的大幅度变化和小相位变化的方法。本发明满足该需要。

发明内容

【0011】在一个实施例中，本发明是一种方法，其用于跟踪在噪声数据中的发射的周期性电磁信号的幅度变化和相位改变，所述噪声数据由至少一个接收器随时间检测，所述信号是以已知频率发射的，所述方法包括的步骤是：a)选择用于跟踪已知频率信号的跟踪算法；b)将所述检测时间分成时间间隔或时间段，在每个时间间隔内检测到的信号和至少一个相关参数被假定不变；c)估计所述检测到的信号和所述至少一个相关参数的初始值，并将这些值分配给第一个时间间隔；d)在时间上在前一个时间间隔估计所述初始信号和每个相关参数的预测；e)使用所述数据和所述跟踪算法来修正步骤d)中的初始估计；以及f)重复步骤d)-e)直到所有数据被处理。

【0012】在本发明的某些实施例中，跟踪算法是Kalman算法，其包括由状态方程和测量方程规定的状态矢量。在这些实施例的部分实施例中，状态矢量具有两个分量，即信号的幅度和正交信号。在其他的实施例中，对于信号经历大幅衰减的情形特别有用，状态矢量具有两个用于更容易地跟踪信号的加性分量：信号包络幅度的改变率和信号的相对相位改变率。

附图说明

【0013】本发明及其优点将通过参考下文详细描述和附图得到更好地理解，在附图中：

图1图解说明了典型的受控源电磁勘查的野外观测系统；

图2是流程图，其显示了Kalman算法中的主要步骤；

图3是流程图，其显示了使用Kalman滤波作为跟踪算法的本发明性方法一个实施例的主要步骤；

图4是本发明更一般的实施例的流程图；

图5-10显示了本发明性方法的结果，在其中Kalman算法被施加到具有0.25Hz信号频率和加性随机噪声的模型数据。

【0014】本发明将结合其优选实施例来描述。然而，下文的具体描述在一定程度上针对本发明的特定实施例或特定应用，其是示意性地而并非限定本发明的范围。相反，本发明趋于覆盖所有的替换，修改和等价表述，其可被包括在所附权利要求的精神和范围之内。

具体实施方式

【0015】本发明是使用诸如Kalman算法的跟踪滤波器来跟踪正弦信号并从电磁噪声中恢复该信号的方法。本文所公开的Kalman滤波方法解决了先前所讨论的现有方法的问题。开始，大部分数据记录可用来获取每个时刻的估计。这很重要，因为相对参考正弦波的低速相位改变，会导致电磁数据在长时间期间的高度相关。换句话说，特定时间处的信息给出了关于其之后很久的信号的信息。另一方面，孤立时间窗口上的傅立叶分析不使用当前窗口以外的任何信息。具体来说，估计的幅度和相位在窗口之间可是不连续的。Kalman方法也可并入信号和噪声特征，例如：远程探测器(或同一探测器上的不同元件)之间的噪声相关性、元件之间的信号相关性、时变信号和噪声幅度改变，以及数据上可预测地质效果。

【0016】为了使用Kalman滤波，该过程必须通过两个线性方程表示：状态方程和测量方程。在适合假设的情形中，Kalman算法给出最小二乘最优信号估计及其关联的误差协方差。线性假设对大多数应用是有效的。线性假设可能失效的例子包括不是加性测量噪声，如某些信号削波或失真或乘性噪声。类似地，例如如果从一个采样到另一个采样信号完全不可预测，或如果信号大小影响转换矩阵或过渡矩阵(transition matrix)而使系统近似不稳定时，则状态方程可能不满足线性假设。轻微非线性的情形仍可使用如下文所做的那样的扩展或其他近似来建模。

【0017】所要求的状态方程含状态矢量x_k，其可以几种电磁数据处理的方式来建立。在最小值处，x_k含两个分量。它们是信号(如在特定位置处的水平电场)及其相应的正交信号。该正交信号是90度相移后的信号。对于正弦信号，正交信号与信号的导数成比例。当正弦曲线是二阶微分方程的解时，需要两个分量。对于每个探测器位置处的每个被估计或估算的信号需要加性分量对。对于每个估计的信号，如果需要，也可建模加性导数。因为对导数的更新给出信号估计的更平滑的校正，因此加性导数是有用的。

【0018】Kalman滤波要求噪声协方差和信号驱动协方差矩阵的指定或规定。例如，噪声协方差项将被用来指示由于MT噪声导致的噪声水平和相关性。信号驱动的协方差项指示所需的滤波调整率和分量之间任意信号相关性。

【0019】在滤波器参数被规定后，在开始滤波算法之前，数据可进行预处理。开始，可缩放数据使得数据预期的信号部分具有相对平坦的幅度。换句话说，使用具有偏移信号衰减速率的粗预测来放大远偏移。滤波算法仅需跟踪来自该预期衰减速率的变化，这比跟踪具有偏移的快速幅度衰减更易管理。由于这些数据被放大来平衡信号，随机噪声也随之被放大。这可在噪声协方差矩阵中规定，以内置在算法中，其对远偏移数据具有更多噪声。

【0020】在本发明的其他实施例中，通过对幅度衰减速率而非幅度自身建模来处理幅度的较大变化。即使信号自身幅度改变几个数量级，衰减(或增加)速率在幅度还是近似的。

【0021】噪声脉冲或丢失数据也可以被识别，因此滤波器将携带正弦信号通过这些区域而无不带所需数据。测量噪声也可被调整从而满足白噪声要求。这些调整包括平滑滤波，去除直流方法，除去谐波噪声的滤波(如果没有被当作信号处理)，和使用分离状态变量(一个或多个)进行有色噪声的建模。

【0022】对于典型方波信号，除了名义基频外还有奇次谐波。这些谐波可被滤出(使用带通滤波器)并作为独立信号处理，或它们可与基频同时建模。如果预期的谐波信号调整要与基频信号调整相关，则同时建模有意义。

【0023】在指定模型和预处理后，可使用Kalman滤波算法作为时间来估计状态矢量分量和相联系的估计误差棒。然后通过比较参数模型或通过将其用作电阻性结构的反演输入(如授予Srnka的美国专利6,603,313中教导的那样)，该最优估计可用在进一步电磁解释(interpretation)中。

【0024】图3显示了本发明一个实施例中的主要步骤。在该实施例中，Kalman滤波器跟踪作为时间函数接收的电磁信号的幅度和相位变化。在下面的讨论中，源和/或接收器被假定是运动的，且因此源相对于接收器的偏移作为时间的函数而变化(如图1中所示，其描述了运动的源)。当应用于静止源和接收器时本发明也能一样良好地发挥作用，即使这不是有效执行CSEM勘查的方式。

【0025】源相对于接收器的偏移的增加导致接收信号的显著衰减。对于跟踪算法，由于预期信号及其修正的幅度可改变几个量级，因此源相对于接收器的偏移的增加是一个潜在的问题。在图3中的步骤31中，执行某些主要步骤从而准备由接收器测量的数据。一个这样的步骤处理了信号在勘查中采用的偏移范围上宽范围变化的问题。

【0026】至少有两种方式处理该问题。在一个方法中，数据可预缩放以补偿典型的偏移幅度降低率。应参考覆盖预期的传导范围的模型结果来确定该幅度衰减。在缩放后，滤波任务被简化，因为其仅跟踪从该基本情况之后的变化，且校正的幅度相对恒定。

【0027】在幅度变化问题的可替换方法中，增加了相应于幅度和相位改变率的状态变量。这些改变率容易建模，因为它们相对于按指数衰减或增加的幅度，在数值上是相对恒定。幅度问题的其他方法也可以设想，包括不对幅度变化做任何处理，且所有方法都在本发明的范围内。

【0028】还是在步骤31中，优选指定或规定噪声协方差矩阵(下面讨论)中的预期噪声，以使得算法可优化使用不同品质的数据。例如，可以是随机噪声脉冲(短时高幅度噪声)。这些可在预处理中作出标记，以使得Kalman滤波器可携带正弦信号通过这些区域而不使用数据。

【0029】另一个必要的预处理步骤包括频率滤波，以平滑噪声谱，这样白加性噪声的假设就是正确的。这通常可能涉及缩小非常低频的分量(即，在基本信号频率之下)，因为环境MT噪声倾向于在这些频率处最大。此外，数据可以被低通滤波以去除较高频的噪声，并允许对较大的样点时间间隔进行在取样或再采样。再采样改进算法的计算需求。

【0030】除了名义基频，典型的方波源信号将包括奇次谐波。这些谐波可被滤出(使用带通滤波器或陷波器)并作为独立信号处理，或它们可与基频同时建模。如果预期谐波信号调整要与基频信号调整相关，则同时建模是优选的。

【0031】对步骤31加以总结，包括上述几种数据准备技术被用在本发明优选的实施例中，但对本发明都不是关键的。

【0032】在图3的步骤32中，建立Kalman滤波器以供应用。在本发明的其他实施例中，可使用不同于Kalman滤波的跟踪算法。图4示出由本发明更一般的实施例执行的基本步骤。在步骤41中，信号和相关联参数的初始估计被输入到算法中。通常，这些初始时间一样点值包括感兴趣的信号及其导数。这些初始值可从具有高信噪比的近偏移数据确定。可替换的，在算法将快速收敛于正确值的假设下，可使用任意初始猜测值。在步骤42中，在使用例如下面结合Kalman滤波讨论的旋转矩阵的例子之前预测(project)这些初始值。在步骤43中，基于测量的数据和特定跟踪算法的指定来调整新样点。步骤44结束循环的一圈，该循环被重复直到所有数据穷尽，即直到对于收集数据的时间期间之前，信号都已经被预测。该解释是简略的，这是因为之后结合Kalman滤波的进一步的细节也将被应用于一般算法。

【0033】Kalman滤波是电磁信号跟踪问题的状态空间公式的优选解决方案。该公式具有两个矩阵方程：“状态”方程和“观测”方程。

状态方程为：

       x_k+1＝Φ_kx_k+w_k          (1)其中x_k是样点k处的状态矢量，Φ_k是状态过渡矩阵，w_k是状态强制函数。时间尺度被分成有限的间隔，且对于每个时间间隔或时间段，每个接收器的测量值被转换为单个数字(下面测量方程中称为z_k)。数据样点k指第k个时间间隔的数字化输出，其中k是表示连续时间段的整数指标或脚标。强制函数是白序列(white sequence)，其表示下一状态矢量样点与应用到当前样点的转换矩阵所预测的样点的差异。在无任何新息情况下，转换矩阵给出下一样点处预测的状态矢量(其中w_k是零)。Szelag方法适用于对振荡信号建模。Szelag使用二元状态矢量，其具有用于振荡信号及其正交信号(与导数成比例)的分量。在本发明的其他实施例中，额外的分量可用来对信号进一步的派生物(derivative)建模。对于两个分量的情形，将在频率f处产生振荡的转

换矩阵如下给出：

$\mathrm{\Φ}=\left[\begin{array}{cc}\cos 2\mathrm{\πfT}& \sin 2\mathrm{\πfT}\\ -\sin 2\mathrm{\πfT}& \cos 2\mathrm{\πfT}\end{array}\right]---\left(2\right)$

其中f是信号频率，且T是样点时间间隔。

【0034】在本发明优选实施例中，该简单的公式被扩展以替代地跟踪幅度和相对相位，这是因为对于典型CSEM问题，幅度和相对相位将随时间逐渐改变(偏移)。该公式使用四元状态矢量：

          x＝[x_s x_q ΔA v]         (3)其中x_s是振荡信号，x_q是正交信号，ΔA是信号的包络幅度改变率，v是信号的相对相位(即频率偏移)改变率。因为幅度和相位相对信号不是线性的，本发明一个实施例中实现了状态方程的小校正线性化。因为v和ΔA的变化预期较小，因此线性假设预期是有效的。

【0035】该线性化处理是通过估计样点(k+1)处的x_s和x_q的值而开始的，给定在样点k处的数值，其是状态矢量的四个元素或要素。如果幅度或相对相位没有变化，方程(2)的简单旋转矩阵Φ给出x_s和x_q在下一样点预测值：

                 x_s(k+1)＝C·x_s(k)+S·x_q(k)，及    (4)

                 x_q(k+1)＝-S·x_s(k)+C·x_q(k)，

其中

                 C＝cos2πfT，及                   (5)

                 S＝sin2πfT。

在本发明的实施例中，接着假设幅度(信号包络)以ΔA每秒的速率增加。换句话说，幅度在进入下一个样点时被乘以(1+T·ΔA)。如果x_s和x_q都被这样的因子所缩放处理，则信号包络将出现该现象：

x_s(k+1)＝(1+T·ΔA)·(C·x_s(k)+S·x_q(k))，及    (6)

x_q(k+1)＝(1+T·ΔA)·(-S·x_s(k)+C·x_q(k))。

接下来，考虑较小的相对相位改变，该相对相位改变率为v弧度/秒。在进入下一个样点时，这将引起相位迁移vT。这可通过修改方程(5)的旋转正弦而并入到状态方程中，如下所示：

$\tilde{C}=\cos (2\mathrm{\πfT}+\mathrm{vT}),$ 及      (7)

$\tilde{S}=\sin (2\mathrm{\πfT}+\mathrm{vT}).$

为了将其线性化，利用vT《1(即vT远小于1)的事实来改写方程(7)如下：

$\tilde{C}=\cos 2\mathrm{\πfT}-\mathrm{vT}\sin 2\mathrm{\πfT},$ 和      (8)

$\tilde{S}=\sin 2\mathrm{\πfT}-\mathrm{vT}\cos 2\mathrm{\πfT}.$

结合方程(6)和(8)并仅保持仅仅一阶校正，从而给出修改的状态方程：

                x_k+1＝Φ_kx_k+w_k        (9)

其中

                x＝[x_s x_q ΔA v]，    (10)

$\mathrm{\Φ}=\left[\begin{array}{cccc}C& S& T({\mathrm{Cx}}_s+{\mathrm{Sx}}_q)& T(-{\mathrm{Sx}}_s+{\mathrm{Cx}}_q)\\ -S& C& T(-S{x}_s+{\mathrm{Cx}}_q)& T(-{\mathrm{Cx}}_s-{\mathrm{Sx}}_q)\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right],$ 及

                w＝[0 0 w_ΔA w_v]应当注意，从恒定幅度正弦曲线的变化仅通过ΔA和v元素发生，因为仅这些元素是w_k非零的。

【0036]也可以指定或规定与w_k关联的协方差矩阵。这是控制滤波器适应速率的方法——w_k中大协方差意味着要求ΔA和v有较大的变化。当几个数据分量被建模时，信号相关性可由状态协方差矩阵中非对角元素指示。

【0037】状态方程的其它修改也是可能的。在最小值处，x_k将包含两个分量。这可能是信号(例如，特殊位置处的水平电场)及其相应的正交信号(与导数成比例)。上面两个加性分量被用于模拟幅度和相位变化。对于每个检测器位置处要估计的每个信号需要加性分量。如果需要，对于每个估计的信号，加性导数可以被建模。加性导数可能是有用的，因为导数的更新可产生对信号估计更平滑的校正。

【0038】上述实施例中Kalman滤波的测量方程由下式给出

           z_k＝H_kx_k+v_k          (11)

其中z_k是在样点k处测量的数据，H＝[1 0 0 0]是测量矩阵，其选择状态矢量中的x_s，且v_k是测量噪声。对于白噪声或近似白噪声，Kalman算法效果很好。如果噪声是窄带的，如正弦的，其可作为独立信号分量建模并被去除。对于v_k，关联的协方差矩阵给出预期的噪声协方差和相关性。方差可以是时变的，如同以缩放的数据工作时那样(即，对于以指数缩放的数据，噪声以指数改变)。特定的噪声区域也可以指定较大的方差以最小化噪声脉冲的效果。对于多分量情形，v_k的协方差矩阵也是包括关于相关的噪声信息的矩阵。例如，当远程探测器含关于MT噪声信息时，这将是有用的。

【0039】在结束图3的步骤32的讨论时，可以看到这里所公开的Kalman滤波方法解决了背景部分指出的问题。作为开始，长有效数据窗口可用来获得每个时刻的估计。由于相对于参考正弦波的低速相位改变，可导致电磁数据在长时间期间的高度相关，因此长有效数据窗口是重要的。换句话说，特定时间处的信息给出关于很长时间以后信号的信息。另一方面，独立时间窗口上的傅立叶分析不使用任何当前窗口外的信息。具体来说，窗口之间的估计的幅度和相位可以是不连续的。Kalman方法也可以并入如下信号和噪声特征：远程探测器(或同一探测器上不同元件)之间的噪声相关性，元件之间的信号识更新，时变信号和噪声幅度变化，以及数据上的地质可预测效应。

【0040】在图3的步骤33中，状态方程，测量方程，和关联的协方差矩阵被用来以图2所示的方式施加Kalman滤波。Kalman滤波通常在一个方向工作。至少有两种方式使用当前样点之前的信息，即时间上在后测量的数据。一个选择是在中途加入数据，并反向运行以获取正向运行的初始估计。另一个选择是采用先前提到的Rauch等人的滤波改进，其包括当前样点之前的所有数据。由于转换矩阵Φ的数据依赖性，Rauch方法中改进的状态方程(9)和(10)存在问题。

【0041】Kalman滤波的输出将是最优的(最小均方误差)状态矢量值和最优的相联系的信号误差协方差矩阵(误差棒)。不同的跟踪算法可使用不同的误差最小化标准。

【0042】再参考图2，以便总结Kalman滤波如何在本发明方法中工作。在步骤21，在某些时间样点处估计状态矢量和其协方差的初始值。在步骤22，Kalman算法通过对所示方程求值，来计算Kalman增益K_k。Kalman增益规定了如何修正测量的数据，以便最好地将其与状态矢量的猜测值相合并或相结合。数据合并是在步骤23中执行的。在步骤24，评估新估计的误差协方差。最后，在步骤25，新估计被用于对下一个样点进行时刻向前预测，然后重复该所述处理。当每个样点被处理时，样点的最小二乘解被在该时刻确定，其通过对求解过程中的方程的简单求解得到。状态方程和测量方程规定了系统状态如何演变及测量如何与系统状态相关联的。图2中协方差矩阵Q_k和R_k相应于状态方程和测量方程中的量w_k和r_k。Q_k的值通常通过反复实验来确定；它们确定滤波器如何快速地对数据变化做出反应。R_k的值是从噪声方差(预期的噪声平方值)和协方差(如果必须，噪声分量之间的相关性)来确定的。除了状态矢量x_k的最小二乘最佳解，Kalman算法也给出与该解相联系的误差协方差。当然，这样的解的最优特性取决于状态方程和测量方程的正确构建，包括用于状态和测量噪声的所需协方差矩阵。

【0043】在解释阶段过程(图3中的步骤34)中，可以几种方式使用状态矢量的最优估计。例如，信号估值x_s可与参数模型比较，以在几种模型化的电阻率结构中选择。可替换地，x_s可用作电阻率结构反演的输入。在参数研究或反演中，ΔA和v也可用来代替x_s使用。

【0044】另一种解释方法是在“快速”1维或2维反演中使用ΔA和v。在这样一个实施例中，这些状态矢量元素被阻挡从而给出分段指数幅度函数，其与简化的地质结构中的单层相对应。

实例

【0045】图5-10中图解说明了模型化的水平电偶极子数据的实例。图5显示了具有加性随机噪声的模型化数据。信号频率是0.25Hz，并从零时间位置进行大致指数衰减。加性噪声从记录的环境地磁噪声中提取。图6示出在约5000秒处，该噪声输入数据(初始信号加加性噪声)61的放大图。在同一图上绘出的是模型数据的无噪声信号部分62和通过使用本发明方法做出的信号的估计63。估计63近似等于无噪声信号62，并难于在图6中与无噪声信号62相区分。由于本发明方法并不“知道”模型信号，因此这种紧密的相似性说明本发明方法取得了成功。其仅被提供噪声数据作为输入，但能恢复信号部分。用在该例子中的本发明实施例采用Kalman滤波作为跟踪算法，否则采用上面所述的一种作为优选方法。相应于图6中所表示的部分的偏移近似为4000米。

【0046】图7示出无噪声信号(实线71)的幅度改变率ΔA和通过使用噪声数据经Kalman滤波获得的ΔA的估计值(锯齿状线72)。应当注意，在信号强度最大的数据中央部分，该估计的质量是最好的。

【0047】图8示出无噪声信号(实线81)相对相位改变率v和通过使用噪声数据经Kalman滤波获得的v的估计值(锯齿状线82)。仍然是在信号强度最大的数据中央部分，该估计的质量是最好的。

【0048】图9示出与Kalman估计(虚线92)相比的无噪声输入数据(实线91)的信号幅度。半对数图显示了信号强度中有3个数量级以上的幅度变化。图10示出输入信号101对Kalman估计102的类似的相对相位比较。

【0049】前面的描述旨在说明本发明的具体实施例，所述具体实施例出于说明本发明的目的。然而，对于本领域技术人员来说，显然这里描述的实施例可能有许多描述和变化。所有这些修正和变化都包含在由所附权利要求所定义的本发明的范围内。

Claims (11)

1.一种方法，其用于跟踪在噪声数据中的发射的周期性电磁信号的幅度变化和相位改变，所述噪声数据由至少一个接收器随时间检测，所述信号是以已知频率发射的，所述方法包括的步骤是：

a)选择用于跟踪已知频率信号的跟踪算法；

b)将所述检测时间分成时间间隔或时间段，在每个时间间隔内检测到的信号和至少一个相关参数被假定不变；

c)估计所述检测到的信号和所述至少一个相关参数的初始值，并将这些值分配给第一个时间间隔；

d)在时间上在前一个时间间隔估计所述初始信号和每个相关参数的预测；

e)使用所述数据和所述跟踪算法来修正步骤d)中的初始估计；以及

f)重复步骤d)-e)直到所有数据被处理。

2.如权利要求1所述的方法，其中所述跟踪算法是Kalman算法，其包括由状态方程和测量方程所规定的状态矢量，所述状态矢量至少具有下列两个分量：所述检测到的信号的幅度和所述至少一个所选择的相关参数。

3.如权利要求2所述的方法，其中所述状态矢量具有两个分量，且所述至少一个相关参数与所述信号的时间导数即正交信号成比例。

4.如权利要求2所述的方法，其中所述状态矢量具有四个分量：所述信号的幅度；所述正交信号幅度；所述信号的包络幅度时变率；以及所述信号的相对相位时变率，用于跟踪信号的所述后两个分量经历显著衰减。

5.如权利要求2所述的方法，其中所述Kalman算法适于在所述修正步骤中使用时间上较晚检测的数据。

6.如权利要求1所述的方法，进一步包括初始数据缩放步骤。

7.如权利要求1所述的方法，其中所述发射的信号是地下地层的电磁勘查中的源信号的傅立叶分量。

8.如权利要求7所述的方法，其包括最终步骤，其用于从信号或至少一个相关参数的估计值来确定所述地下地层的电阻率结构。

9.如权利要求2所述的方法，其中在所述状态方程中实现小校正线性化。

10.如权利要求2所述的方法，其中为每个时间间隔获取所述信号的估计及与其相联系的误差。

11.如权利要求4所述的方法，其中通过使用小校正线性化假设，修正所述Kalman算法以允许将所述信号的包络幅度的时变率A，和所述信号的相对相位的时变率v加到所述状态矢量中，所述线性化假设包括：(i)所述信号包络的所述幅度被假定乘以(1+T·A)，其中A是所述信号包络的幅度时变率，而T是时间间隔，以及(ii)所述信号假设经历相移vT，其中v是所述信号的相对相位改变率，并且假定vT＜＜1。

Images