專(zhuān)利名稱(chēng)：：處理電磁數(shù)據(jù)的Kalman濾波方法
技術(shù)領(lǐng)域：
：0002本發(fā)明一般涉及地球物理勘探領(lǐng)域，更具體地涉及用來(lái)探測(cè)碳?xì)浠衔锏碾姶欧椒?。具體來(lái)說(shuō)，本發(fā)明是用于跟蹤用于受控的源電磁勘探的電磁源信號(hào)，以從噪聲恢復(fù)該信號(hào)的方法。
背景技術(shù)：
：0003受控源電磁(“CSEM”)地球物理勘查使用活性(人造)源來(lái)產(chǎn)生電磁場(chǎng)激勵(lì)地球，并配置地球表面，海床上，或鉆孔內(nèi)的接收器儀器，以測(cè)量所產(chǎn)生的電場(chǎng)和磁場(chǎng)，即地球?qū)υ醇?lì)的響應(yīng)。圖1圖解說(shuō)明了海岸CSEM勘查的基本元件。船只在海床13上方拖曳浸入水中的CSEM發(fā)射器11。然后，分析由接收器12測(cè)量的電場(chǎng)和磁場(chǎng)，以確定表面或海床下面地球結(jié)構(gòu)(地表下地層)的電阻率。該技術(shù)已經(jīng)應(yīng)用陸上礦物勘探，海洋構(gòu)造研究，和海岸石油及礦物資源勘探。0004活性電磁源信號(hào)可被當(dāng)作正弦信號(hào)的和處理(如，由帶有奇次諧波的基頻組成的方波信號(hào))。該源的例子是用于許多CSEM工作的水平電偶極子。隨著偏移(即該偶極子源11和接收器12之間的距離)的增加，正弦信號(hào)可能顯著衰減。而且，遠(yuǎn)偏移對(duì)于確定感興趣的深阻結(jié)構(gòu)來(lái)說(shuō)經(jīng)常是關(guān)鍵。結(jié)果，存在使該正弦信號(hào)獲得優(yōu)良信噪比的需求。0005典型的改善該EM數(shù)據(jù)的信號(hào)噪聲的處理方法包括將數(shù)據(jù)分解到多個(gè)時(shí)間窗口中，在這些時(shí)間窗口上使用傅立葉分析或相似的方法，來(lái)計(jì)算所選頻率分量(一個(gè)或多個(gè))的幅度和相位。例如，可見(jiàn)Constable和Cox的“Marinecontrolled-sourceelectromagneticsounding2.ThePEGASUSExperiment，”JournalofGeophysicalresearch101，5519-5530(1996)。這些窗口不能太長(zhǎng)，這是因?yàn)樾盘?hào)幅度和相對(duì)相位可能在該分析窗口中顯著改變。然而，短窗僅允許最小的信噪比改善。當(dāng)前的方法要求在這兩個(gè)極端之間平衡。0006現(xiàn)有方法的另一個(gè)問(wèn)題是他們不利用信號(hào)和噪聲相關(guān)性。具體來(lái)說(shuō)，對(duì)于活性源海洋EM成像，由于低頻大地電磁(“MT”)噪聲可偽裝成信號(hào)，因此低頻MT噪聲就是一個(gè)特別重要的難題。(MT噪聲是來(lái)自自然非活性源的電磁發(fā)射。)不同探測(cè)器之間的相關(guān)性可用來(lái)幫助從這些噪聲中分離活性源信號(hào)。當(dāng)前方法中，沒(méi)有最優(yōu)地使用其他信號(hào)和噪聲的相關(guān)性(如，兩個(gè)水平分量信號(hào)的相關(guān)性)。0007Kalman(卡爾曼)濾波算法起源于導(dǎo)航定位問(wèn)題并特別適用于跟蹤類(lèi)的問(wèn)題(Kalman，1960)。最初由Kalman發(fā)表于《Trans.OftheASME-J.OFBasicEngr.》35-45(1960)上，還發(fā)表了大量文獻(xiàn)作為關(guān)于基本Kalman濾波的改進(jìn)和應(yīng)用的總結(jié)，例如在由Brown在紐約JohnWiley&Sons(1983)上發(fā)表的IntroductiontoRandomSignalAnalysisandKalmanFiltering中。這些改進(jìn)中的幾種對(duì)于當(dāng)前發(fā)明的某些實(shí)施例是重要的。0008標(biāo)準(zhǔn)Kalman濾波在一個(gè)方向上運(yùn)行，并濾波該方向(或時(shí)間)序列中的數(shù)據(jù)。因此僅前面的數(shù)據(jù)影響濾波結(jié)果。由Rauch等人提出的重要修改給出最優(yōu)處理，其使用整個(gè)時(shí)間記錄參看Rauch的“Solutionstothelinearsmoothingproblem”，IEEETrans.OnAuto.Control，AC-8，371(1963)；和Rauch等人的“Maximumlikelihoodestimatesoflineardynamicsystems”，AMAJ.3，1445(1965)。Szelag公開(kāi)了另一個(gè)算法改進(jìn)，其允許濾波器跟蹤已知頻率的正弦信號(hào)；參見(jiàn)“Ashorttermforecastingalgorithmfortrunkdemandservicing”，TheBellSystemTechnicalJournal61，67-96(1982)。其被發(fā)展，用于以電話(huà)中繼負(fù)載值(telephonetrunkloadvalues)來(lái)跟蹤年度周期。0009圖2是流程圖，其圖解說(shuō)明了Kalman算法。更詳細(xì)的可參考Brown的論文集第200頁(yè)。0010LaScala等人公開(kāi)了公知的用于跟蹤時(shí)變頻率的擴(kuò)展kalman濾波器的使用。(“DesignofanextendedKalmanfilterfrequencytracker”，IEEETransactiononSignalProcessing44，No.3，739-742(1996年3月))所述公式假定信號(hào)幅度保持恒定。所用的特殊Kalman算法旨在跟蹤未知頻率的信號(hào)，其中頻率可經(jīng)歷相當(dāng)?shù)母淖儭agunas等人公開(kāi)了一種擴(kuò)展的Kalman濾波器，用于跟蹤有噪聲和頻率變化的復(fù)雜正弦，如Doppler移動(dòng)。如同LaScala的方法，Lagunas的算法被設(shè)計(jì)來(lái)跟蹤未知頻率的正弦曲線(xiàn)或正弦波。因此，如果應(yīng)用來(lái)跟蹤具有恒定或近似恒定頻率的信號(hào)，這兩種方法都是次優(yōu)的。如果改變相對(duì)小，則Lagunas的方法能夠跟蹤幅度變化。本發(fā)明不是旨在處理使用發(fā)射已知頻率波形的電磁源的電磁勘探數(shù)據(jù)。需要使用大窗口或甚至所有電磁數(shù)據(jù)，來(lái)跟蹤關(guān)于已知正弦曲線(xiàn)的大幅度變化和小相位變化的方法。本發(fā)明滿(mǎn)足該需要。
發(fā)明內(nèi)容0011在一個(gè)實(shí)施例中，本發(fā)明是一種方法，其用于跟蹤在噪聲數(shù)據(jù)中的發(fā)射的周期性電磁信號(hào)的幅度變化和相位改變，所述噪聲數(shù)據(jù)由至少一個(gè)接收器隨時(shí)間檢測(cè)，所述信號(hào)是以已知頻率發(fā)射的，所述方法包括的步驟是a)選擇用于跟蹤已知頻率信號(hào)的跟蹤算法；b)將所述檢測(cè)時(shí)間分成時(shí)間間隔或時(shí)間段，在每個(gè)時(shí)間間隔內(nèi)檢測(cè)到的信號(hào)和至少一個(gè)相關(guān)參數(shù)被假定不變；c)估計(jì)所述檢測(cè)到的信號(hào)和所述至少一個(gè)相關(guān)參數(shù)的初始值，并將這些值分配給第一個(gè)時(shí)間間隔；d)在時(shí)間上在前一個(gè)時(shí)間間隔估計(jì)所述初始信號(hào)和每個(gè)相關(guān)參數(shù)的預(yù)測(cè)；e)使用所述數(shù)據(jù)和所述跟蹤算法來(lái)修正步驟d)中的初始估計(jì)；以及f)重復(fù)步驟d)-e)直到所有數(shù)據(jù)被處理。0012在本發(fā)明的某些實(shí)施例中，跟蹤算法是Kalman算法，其包括由狀態(tài)方程和測(cè)量方程規(guī)定的狀態(tài)矢量。在這些實(shí)施例的部分實(shí)施例中，狀態(tài)矢量具有兩個(gè)分量，即信號(hào)的幅度和正交信號(hào)。在其他的實(shí)施例中，對(duì)于信號(hào)經(jīng)歷大幅衰減的情形特別有用，狀態(tài)矢量具有兩個(gè)用于更容易地跟蹤信號(hào)的加性分量信號(hào)包絡(luò)幅度的改變率和信號(hào)的相對(duì)相位改變率。0013本發(fā)明及其優(yōu)點(diǎn)將通過(guò)參考下文詳細(xì)描述和附圖得到更好地理解，在附圖中圖1圖解說(shuō)明了典型的受控源電磁勘查的野外觀(guān)測(cè)系統(tǒng)；圖2是流程圖，其顯示了Kalman算法中的主要步驟；圖3是流程圖，其顯示了使用Kalman濾波作為跟蹤算法的本發(fā)明性方法一個(gè)實(shí)施例的主要步驟；圖4是本發(fā)明更一般的實(shí)施例的流程圖；圖5-10顯示了本發(fā)明性方法的結(jié)果，在其中Kalman算法被施加到具有0.25Hz信號(hào)頻率和加性隨機(jī)噪聲的模型數(shù)據(jù)。0014本發(fā)明將結(jié)合其優(yōu)選實(shí)施例來(lái)描述。然而，下文的具體描述在一定程度上針對(duì)本發(fā)明的特定實(shí)施例或特定應(yīng)用，其是示意性地而并非限定本發(fā)明的范圍。相反，本發(fā)明趨于覆蓋所有的替換，修改和等價(jià)表述，其可被包括在所附權(quán)利要求的精神和范圍之內(nèi)。具體實(shí)施例方式0015本發(fā)明是使用諸如Kalman算法的跟蹤濾波器來(lái)跟蹤正弦信號(hào)并從電磁噪聲中恢復(fù)該信號(hào)的方法。本文所公開(kāi)的Kalman濾波方法解決了先前所討論的現(xiàn)有方法的問(wèn)題。開(kāi)始，大部分?jǐn)?shù)據(jù)記錄可用來(lái)獲取每個(gè)時(shí)刻的估計(jì)。這很重要，因?yàn)橄鄬?duì)參考正弦波的低速相位改變，會(huì)導(dǎo)致電磁數(shù)據(jù)在長(zhǎng)時(shí)間期間的高度相關(guān)。換句話(huà)說(shuō)，特定時(shí)間處的信息給出了關(guān)于其之后很久的信號(hào)的信息。另一方面，孤立時(shí)間窗口上的傅立葉分析不使用當(dāng)前窗口以外的任何信息。具體來(lái)說(shuō)，估計(jì)的幅度和相位在窗口之間可是不連續(xù)的。Kalman方法也可并入信號(hào)和噪聲特征，例如遠(yuǎn)程探測(cè)器(或同一探測(cè)器上的不同元件)之間的噪聲相關(guān)性、元件之間的信號(hào)相關(guān)性、時(shí)變信號(hào)和噪聲幅度改變，以及數(shù)據(jù)上可預(yù)測(cè)地質(zhì)效果。0016為了使用Kalman濾波，該過(guò)程必須通過(guò)兩個(gè)線(xiàn)性方程表示狀態(tài)方程和測(cè)量方程。在適合假設(shè)的情形中，Kalman算法給出最小二乘最優(yōu)信號(hào)估計(jì)及其關(guān)聯(lián)的誤差協(xié)方差。線(xiàn)性假設(shè)對(duì)大多數(shù)應(yīng)用是有效的。線(xiàn)性假設(shè)可能失效的例子包括不是加性測(cè)量噪聲，如某些信號(hào)削波或失真或乘性噪聲。類(lèi)似地，例如如果從一個(gè)采樣到另一個(gè)采樣信號(hào)完全不可預(yù)測(cè)，或如果信號(hào)大小影響轉(zhuǎn)換矩陣或過(guò)渡矩陣(transitionmatrix)而使系統(tǒng)近似不穩(wěn)定時(shí)，則狀態(tài)方程可能不滿(mǎn)足線(xiàn)性假設(shè)。輕微非線(xiàn)性的情形仍可使用如下文所做的那樣的擴(kuò)展或其他近似來(lái)建模。0017所要求的狀態(tài)方程含狀態(tài)矢量xk，其可以幾種電磁數(shù)據(jù)處理的方式來(lái)建立。在最小值處，xk含兩個(gè)分量。它們是信號(hào)(如在特定位置處的水平電場(chǎng))及其相應(yīng)的正交信號(hào)。該正交信號(hào)是90度相移后的信號(hào)。對(duì)于正弦信號(hào)，正交信號(hào)與信號(hào)的導(dǎo)數(shù)成比例。當(dāng)正弦曲線(xiàn)是二階微分方程的解時(shí)，需要兩個(gè)分量。對(duì)于每個(gè)探測(cè)器位置處的每個(gè)被估計(jì)或估算的信號(hào)需要加性分量對(duì)。對(duì)于每個(gè)估計(jì)的信號(hào)，如果需要，也可建模加性導(dǎo)數(shù)。因?yàn)閷?duì)導(dǎo)數(shù)的更新給出信號(hào)估計(jì)的更平滑的校正，因此加性導(dǎo)數(shù)是有用的。0018Kalman濾波要求噪聲協(xié)方差和信號(hào)驅(qū)動(dòng)協(xié)方差矩陣的指定或規(guī)定。例如，噪聲協(xié)方差項(xiàng)將被用來(lái)指示由于MT噪聲導(dǎo)致的噪聲水平和相關(guān)性。信號(hào)驅(qū)動(dòng)的協(xié)方差項(xiàng)指示所需的濾波調(diào)整率和分量之間任意信號(hào)相關(guān)性。0019在濾波器參數(shù)被規(guī)定后，在開(kāi)始濾波算法之前，數(shù)據(jù)可進(jìn)行預(yù)處理。開(kāi)始，可縮放數(shù)據(jù)使得數(shù)據(jù)預(yù)期的信號(hào)部分具有相對(duì)平坦的幅度。換句話(huà)說(shuō)，使用具有偏移信號(hào)衰減速率的粗預(yù)測(cè)來(lái)放大遠(yuǎn)偏移。濾波算法僅需跟蹤來(lái)自該預(yù)期衰減速率的變化，這比跟蹤具有偏移的快速幅度衰減更易管理。由于這些數(shù)據(jù)被放大來(lái)平衡信號(hào)，隨機(jī)噪聲也隨之被放大。這可在噪聲協(xié)方差矩陣中規(guī)定，以?xún)?nèi)置在算法中，其對(duì)遠(yuǎn)偏移數(shù)據(jù)具有更多噪聲。0020在本發(fā)明的其他實(shí)施例中，通過(guò)對(duì)幅度衰減速率而非幅度自身建模來(lái)處理幅度的較大變化。即使信號(hào)自身幅度改變幾個(gè)數(shù)量級(jí)，衰減(或增加)速率在幅度還是近似的。0021噪聲脈沖或丟失數(shù)據(jù)也可以被識(shí)別，因此濾波器將攜帶正弦信號(hào)通過(guò)這些區(qū)域而無(wú)不帶所需數(shù)據(jù)。測(cè)量噪聲也可被調(diào)整從而滿(mǎn)足白噪聲要求。這些調(diào)整包括平滑濾波，去除直流方法，除去諧波噪聲的濾波(如果沒(méi)有被當(dāng)作信號(hào)處理)，和使用分離狀態(tài)變量(一個(gè)或多個(gè))進(jìn)行有色噪聲的建模。0022對(duì)于典型方波信號(hào)，除了名義基頻外還有奇次諧波。這些諧波可被濾出(使用帶通濾波器)并作為獨(dú)立信號(hào)處理，或它們可與基頻同時(shí)建模。如果預(yù)期的諧波信號(hào)調(diào)整要與基頻信號(hào)調(diào)整相關(guān)，則同時(shí)建模有意義。0023在指定模型和預(yù)處理后，可使用Kalman濾波算法作為時(shí)間來(lái)估計(jì)狀態(tài)矢量分量和相聯(lián)系的估計(jì)誤差棒。然后通過(guò)比較參數(shù)模型或通過(guò)將其用作電阻性結(jié)構(gòu)的反演輸入(如授予Srnka的美國(guó)專(zhuān)利6,603,313中教導(dǎo)的那樣)，該最優(yōu)估計(jì)可用在進(jìn)一步電磁解釋(interpretation)中。0024圖3顯示了本發(fā)明一個(gè)實(shí)施例中的主要步驟。在該實(shí)施例中，Kalman濾波器跟蹤作為時(shí)間函數(shù)接收的電磁信號(hào)的幅度和相位變化。在下面的討論中，源和/或接收器被假定是運(yùn)動(dòng)的，且因此源相對(duì)于接收器的偏移作為時(shí)間的函數(shù)而變化(如圖1中所示，其描述了運(yùn)動(dòng)的源)。當(dāng)應(yīng)用于靜止源和接收器時(shí)本發(fā)明也能一樣良好地發(fā)揮作用，即使這不是有效執(zhí)行CSEM勘查的方式。0025源相對(duì)于接收器的偏移的增加導(dǎo)致接收信號(hào)的顯著衰減。對(duì)于跟蹤算法，由于預(yù)期信號(hào)及其修正的幅度可改變幾個(gè)量級(jí)，因此源相對(duì)于接收器的偏移的增加是一個(gè)潛在的問(wèn)題。在圖3中的步驟31中，執(zhí)行某些主要步驟從而準(zhǔn)備由接收器測(cè)量的數(shù)據(jù)。一個(gè)這樣的步驟處理了信號(hào)在勘查中采用的偏移范圍上寬范圍變化的問(wèn)題。0026至少有兩種方式處理該問(wèn)題。在一個(gè)方法中，數(shù)據(jù)可預(yù)縮放以補(bǔ)償?shù)湫偷钠品冉档吐?。?yīng)參考覆蓋預(yù)期的傳導(dǎo)范圍的模型結(jié)果來(lái)確定該幅度衰減。在縮放后，濾波任務(wù)被簡(jiǎn)化，因?yàn)槠鋬H跟蹤從該基本情況之后的變化，且校正的幅度相對(duì)恒定。0027在幅度變化問(wèn)題的可替換方法中，增加了相應(yīng)于幅度和相位改變率的狀態(tài)變量。這些改變率容易建模，因?yàn)樗鼈兿鄬?duì)于按指數(shù)衰減或增加的幅度，在數(shù)值上是相對(duì)恒定。幅度問(wèn)題的其他方法也可以設(shè)想，包括不對(duì)幅度變化做任何處理，且所有方法都在本發(fā)明的范圍內(nèi)。0028還是在步驟31中，優(yōu)選指定或規(guī)定噪聲協(xié)方差矩陣(下面討論)中的預(yù)期噪聲，以使得算法可優(yōu)化使用不同品質(zhì)的數(shù)據(jù)。例如，可以是隨機(jī)噪聲脈沖(短時(shí)高幅度噪聲)。這些可在預(yù)處理中作出標(biāo)記，以使得Kalman濾波器可攜帶正弦信號(hào)通過(guò)這些區(qū)域而不使用數(shù)據(jù)。0029另一個(gè)必要的預(yù)處理步驟包括頻率濾波，以平滑噪聲譜，這樣白加性噪聲的假設(shè)就是正確的。這通?？赡苌婕翱s小非常低頻的分量(即，在基本信號(hào)頻率之下)，因?yàn)榄h(huán)境MT噪聲傾向于在這些頻率處最大。此外，數(shù)據(jù)可以被低通濾波以去除較高頻的噪聲，并允許對(duì)較大的樣點(diǎn)時(shí)間間隔進(jìn)行在取樣或再采樣。再采樣改進(jìn)算法的計(jì)算需求。0030除了名義基頻，典型的方波源信號(hào)將包括奇次諧波。這些諧波可被濾出(使用帶通濾波器或陷波器)并作為獨(dú)立信號(hào)處理，或它們可與基頻同時(shí)建模。如果預(yù)期諧波信號(hào)調(diào)整要與基頻信號(hào)調(diào)整相關(guān)，則同時(shí)建模是優(yōu)選的。0031對(duì)步驟31加以總結(jié)，包括上述幾種數(shù)據(jù)準(zhǔn)備技術(shù)被用在本發(fā)明優(yōu)選的實(shí)施例中，但對(duì)本發(fā)明都不是關(guān)鍵的。0032在圖3的步驟32中，建立Kalman濾波器以供應(yīng)用。在本發(fā)明的其他實(shí)施例中，可使用不同于Kalman濾波的跟蹤算法。圖4示出由本發(fā)明更一般的實(shí)施例執(zhí)行的基本步驟。在步驟41中，信號(hào)和相關(guān)聯(lián)參數(shù)的初始估計(jì)被輸入到算法中。通常，這些初始時(shí)間一樣點(diǎn)值包括感興趣的信號(hào)及其導(dǎo)數(shù)。這些初始值可從具有高信噪比的近偏移數(shù)據(jù)確定?？商鎿Q的，在算法將快速收斂于正確值的假設(shè)下，可使用任意初始猜測(cè)值。在步驟42中，在使用例如下面結(jié)合Kalman濾波討論的旋轉(zhuǎn)矩陣的例子之前預(yù)測(cè)(project)這些初始值。在步驟43中，基于測(cè)量的數(shù)據(jù)和特定跟蹤算法的指定來(lái)調(diào)整新樣點(diǎn)。步驟44結(jié)束循環(huán)的一圈，該循環(huán)被重復(fù)直到所有數(shù)據(jù)窮盡，即直到對(duì)于收集數(shù)據(jù)的時(shí)間期間之前，信號(hào)都已經(jīng)被預(yù)測(cè)。該解釋是簡(jiǎn)略的，這是因?yàn)橹蠼Y(jié)合Kalman濾波的進(jìn)一步的細(xì)節(jié)也將被應(yīng)用于一般算法。0033Kalman濾波是電磁信號(hào)跟蹤問(wèn)題的狀態(tài)空間公式的優(yōu)選解決方案。該公式具有兩個(gè)矩陣方程“狀態(tài)”方程和“觀(guān)測(cè)”方程。狀態(tài)方程為xk+1＝Φkxk+wk(1)其中xk是樣點(diǎn)k處的狀態(tài)矢量，Φk是狀態(tài)過(guò)渡矩陣，wk是狀態(tài)強(qiáng)制函數(shù)。時(shí)間尺度被分成有限的間隔，且對(duì)于每個(gè)時(shí)間間隔或時(shí)間段，每個(gè)接收器的測(cè)量值被轉(zhuǎn)換為單個(gè)數(shù)字(下面測(cè)量方程中稱(chēng)為zk)。數(shù)據(jù)樣點(diǎn)k指第k個(gè)時(shí)間間隔的數(shù)字化輸出，其中k是表示連續(xù)時(shí)間段的整數(shù)指標(biāo)或腳標(biāo)。強(qiáng)制函數(shù)是白序列(whitesequence)，其表示下一狀態(tài)矢量樣點(diǎn)與應(yīng)用到當(dāng)前樣點(diǎn)的轉(zhuǎn)換矩陣所預(yù)測(cè)的樣點(diǎn)的差異。在無(wú)任何新息情況下，轉(zhuǎn)換矩陣給出下一樣點(diǎn)處預(yù)測(cè)的狀態(tài)矢量(其中wk是零)。Szelag方法適用于對(duì)振蕩信號(hào)建模。Szelag使用二元狀態(tài)矢量，其具有用于振蕩信號(hào)及其正交信號(hào)(與導(dǎo)數(shù)成比例)的分量。在本發(fā)明的其他實(shí)施例中，額外的分量可用來(lái)對(duì)信號(hào)進(jìn)一步的派生物(derivative)建模。對(duì)于兩個(gè)分量的情形，將在頻率f處產(chǎn)生振蕩的轉(zhuǎn)換矩陣如下給出Φ=cos2πfTsin2πfT-sin2πfTcos2πfT---(2)]]>其中f是信號(hào)頻率，且T是樣點(diǎn)時(shí)間間隔。0034在本發(fā)明優(yōu)選實(shí)施例中，該簡(jiǎn)單的公式被擴(kuò)展以替代地跟蹤幅度和相對(duì)相位，這是因?yàn)閷?duì)于典型CSEM問(wèn)題，幅度和相對(duì)相位將隨時(shí)間逐漸改變(偏移)。該公式使用四元狀態(tài)矢量x＝[xsxqΔAv](3)其中xs是振蕩信號(hào)，xq是正交信號(hào)，ΔA是信號(hào)的包絡(luò)幅度改變率，v是信號(hào)的相對(duì)相位(即頻率偏移)改變率。因?yàn)榉群拖辔幌鄬?duì)信號(hào)不是線(xiàn)性的，本發(fā)明一個(gè)實(shí)施例中實(shí)現(xiàn)了狀態(tài)方程的小校正線(xiàn)性化。因?yàn)関和ΔA的變化預(yù)期較小，因此線(xiàn)性假設(shè)預(yù)期是有效的。0035該線(xiàn)性化處理是通過(guò)估計(jì)樣點(diǎn)(k+1)處的xs和xq的值而開(kāi)始的，給定在樣點(diǎn)k處的數(shù)值，其是狀態(tài)矢量的四個(gè)元素或要素。如果幅度或相對(duì)相位沒(méi)有變化，方程(2)的簡(jiǎn)單旋轉(zhuǎn)矩陣Φ給出xs和xq在下一樣點(diǎn)預(yù)測(cè)值xs(k+1)＝C·xs(k)+S·xq(k)，及(4)xq(k+1)＝-S·xs(k)+C·xq(k)，其中C＝cos2πfT，及(5)S＝sin2πfT。在本發(fā)明的實(shí)施例中，接著假設(shè)幅度(信號(hào)包絡(luò))以ΔA每秒的速率增加。換句話(huà)說(shuō)，幅度在進(jìn)入下一個(gè)樣點(diǎn)時(shí)被乘以(1+T·ΔA)。如果xs和xq都被這樣的因子所縮放處理，則信號(hào)包絡(luò)將出現(xiàn)該現(xiàn)象xs(k+1)＝(1+T·ΔA)·(C·xs(k)+S·xq(k))，及(6)xq(k+1)＝(1+T·ΔA)·(-S·xs(k)+C·xq(k))。接下來(lái)，考慮較小的相對(duì)相位改變，該相對(duì)相位改變率為v弧度/秒。在進(jìn)入下一個(gè)樣點(diǎn)時(shí)，這將引起相位遷移vT。這可通過(guò)修改方程(5)的旋轉(zhuǎn)正弦而并入到狀態(tài)方程中，如下所示C~=cos(2πfT+vT),]]>及(7)S~=sin(2πfT+vT).]]>為了將其線(xiàn)性化，利用vT《1(即vT遠(yuǎn)小于1)的事實(shí)來(lái)改寫(xiě)方程(7)如下C~=cos2πfT-vTsin2πfT,]]>和(8)S~=sin2πfT-vTcos2πfT.]]>結(jié)合方程(6)和(8)并僅保持僅僅一階校正，從而給出修改的狀態(tài)方程xk+1＝Φkxk+wk(9)其中x＝[xsxqΔAv]，(10)Φ=CST(Cxs+Sxq)T(-Sxs+Cxq)-SCT(-Sxs+Cxq)T(-Cxs-Sxq)00100001,]]>及w＝應(yīng)當(dāng)注意，從恒定幅度正弦曲線(xiàn)的變化僅通過(guò)ΔA和v元素發(fā)生，因?yàn)閮H這些元素是wk非零的。0036]也可以指定或規(guī)定與wk關(guān)聯(lián)的協(xié)方差矩陣。這是控制濾波器適應(yīng)速率的方法——wk中大協(xié)方差意味著要求ΔA和v有較大的變化。當(dāng)幾個(gè)數(shù)據(jù)分量被建模時(shí)，信號(hào)相關(guān)性可由狀態(tài)協(xié)方差矩陣中非對(duì)角元素指示。0037狀態(tài)方程的其它修改也是可能的。在最小值處，xk將包含兩個(gè)分量。這可能是信號(hào)(例如，特殊位置處的水平電場(chǎng))及其相應(yīng)的正交信號(hào)(與導(dǎo)數(shù)成比例)。上面兩個(gè)加性分量被用于模擬幅度和相位變化。對(duì)于每個(gè)檢測(cè)器位置處要估計(jì)的每個(gè)信號(hào)需要加性分量。如果需要，對(duì)于每個(gè)估計(jì)的信號(hào)，加性導(dǎo)數(shù)可以被建模。加性導(dǎo)數(shù)可能是有用的，因?yàn)閷?dǎo)數(shù)的更新可產(chǎn)生對(duì)信號(hào)估計(jì)更平滑的校正。0038上述實(shí)施例中Kalman濾波的測(cè)量方程由下式給出zk＝Hkxk+vk(11)其中zk是在樣點(diǎn)k處測(cè)量的數(shù)據(jù)，H＝[1000]是測(cè)量矩陣，其選擇狀態(tài)矢量中的xs，且vk是測(cè)量噪聲。對(duì)于白噪聲或近似白噪聲，Kalman算法效果很好。如果噪聲是窄帶的，如正弦的，其可作為獨(dú)立信號(hào)分量建模并被去除。對(duì)于vk，關(guān)聯(lián)的協(xié)方差矩陣給出預(yù)期的噪聲協(xié)方差和相關(guān)性。方差可以是時(shí)變的，如同以縮放的數(shù)據(jù)工作時(shí)那樣(即，對(duì)于以指數(shù)縮放的數(shù)據(jù)，噪聲以指數(shù)改變)。特定的噪聲區(qū)域也可以指定較大的方差以最小化噪聲脈沖的效果。對(duì)于多分量情形，vk的協(xié)方差矩陣也是包括關(guān)于相關(guān)的噪聲信息的矩陣。例如，當(dāng)遠(yuǎn)程探測(cè)器含關(guān)于MT噪聲信息時(shí)，這將是有用的。0039在結(jié)束圖3的步驟32的討論時(shí)，可以看到這里所公開(kāi)的Kalman濾波方法解決了背景部分指出的問(wèn)題。作為開(kāi)始，長(zhǎng)有效數(shù)據(jù)窗口可用來(lái)獲得每個(gè)時(shí)刻的估計(jì)。由于相對(duì)于參考正弦波的低速相位改變，可導(dǎo)致電磁數(shù)據(jù)在長(zhǎng)時(shí)間期間的高度相關(guān)，因此長(zhǎng)有效數(shù)據(jù)窗口是重要的。換句話(huà)說(shuō)，特定時(shí)間處的信息給出關(guān)于很長(zhǎng)時(shí)間以后信號(hào)的信息。另一方面，獨(dú)立時(shí)間窗口上的傅立葉分析不使用任何當(dāng)前窗口外的信息。具體來(lái)說(shuō)，窗口之間的估計(jì)的幅度和相位可以是不連續(xù)的。Kalman方法也可以并入如下信號(hào)和噪聲特征遠(yuǎn)程探測(cè)器(或同一探測(cè)器上不同元件)之間的噪聲相關(guān)性，元件之間的信號(hào)識(shí)更新，時(shí)變信號(hào)和噪聲幅度變化，以及數(shù)據(jù)上的地質(zhì)可預(yù)測(cè)效應(yīng)。0040在圖3的步驟33中，狀態(tài)方程，測(cè)量方程，和關(guān)聯(lián)的協(xié)方差矩陣被用來(lái)以圖2所示的方式施加Kalman濾波。Kalman濾波通常在一個(gè)方向工作。至少有兩種方式使用當(dāng)前樣點(diǎn)之前的信息，即時(shí)間上在后測(cè)量的數(shù)據(jù)。一個(gè)選擇是在中途加入數(shù)據(jù)，并反向運(yùn)行以獲取正向運(yùn)行的初始估計(jì)。另一個(gè)選擇是采用先前提到的Rauch等人的濾波改進(jìn)，其包括當(dāng)前樣點(diǎn)之前的所有數(shù)據(jù)。由于轉(zhuǎn)換矩陣Φ的數(shù)據(jù)依賴(lài)性，Rauch方法中改進(jìn)的狀態(tài)方程(9)和(10)存在問(wèn)題。0041Kalman濾波的輸出將是最優(yōu)的(最小均方誤差)狀態(tài)矢量值和最優(yōu)的相聯(lián)系的信號(hào)誤差協(xié)方差矩陣(誤差棒)。不同的跟蹤算法可使用不同的誤差最小化標(biāo)準(zhǔn)。0042再參考圖2，以便總結(jié)Kalman濾波如何在本發(fā)明方法中工作。在步驟21，在某些時(shí)間樣點(diǎn)處估計(jì)狀態(tài)矢量和其協(xié)方差的初始值。在步驟22，Kalman算法通過(guò)對(duì)所示方程求值，來(lái)計(jì)算Kalman增益Kk。Kalman增益規(guī)定了如何修正測(cè)量的數(shù)據(jù)，以便最好地將其與狀態(tài)矢量的猜測(cè)值相合并或相結(jié)合。數(shù)據(jù)合并是在步驟23中執(zhí)行的。在步驟24，評(píng)估新估計(jì)的誤差協(xié)方差。最后，在步驟25，新估計(jì)被用于對(duì)下一個(gè)樣點(diǎn)進(jìn)行時(shí)刻向前預(yù)測(cè)，然后重復(fù)該所述處理。當(dāng)每個(gè)樣點(diǎn)被處理時(shí)，樣點(diǎn)的最小二乘解被在該時(shí)刻確定，其通過(guò)對(duì)求解過(guò)程中的方程的簡(jiǎn)單求解得到。狀態(tài)方程和測(cè)量方程規(guī)定了系統(tǒng)狀態(tài)如何演變及測(cè)量如何與系統(tǒng)狀態(tài)相關(guān)聯(lián)的。圖2中協(xié)方差矩陣Qk和Rk相應(yīng)于狀態(tài)方程和測(cè)量方程中的量wk和rk。Qk的值通常通過(guò)反復(fù)實(shí)驗(yàn)來(lái)確定；它們確定濾波器如何快速地對(duì)數(shù)據(jù)變化做出反應(yīng)。Rk的值是從噪聲方差(預(yù)期的噪聲平方值)和協(xié)方差(如果必須，噪聲分量之間的相關(guān)性)來(lái)確定的。除了狀態(tài)矢量xk的最小二乘最佳解，Kalman算法也給出與該解相聯(lián)系的誤差協(xié)方差。當(dāng)然，這樣的解的最優(yōu)特性取決于狀態(tài)方程和測(cè)量方程的正確構(gòu)建，包括用于狀態(tài)和測(cè)量噪聲的所需協(xié)方差矩陣。0043在解釋階段過(guò)程(圖3中的步驟34)中，可以幾種方式使用狀態(tài)矢量的最優(yōu)估計(jì)。例如，信號(hào)估值xs可與參數(shù)模型比較，以在幾種模型化的電阻率結(jié)構(gòu)中選擇?？商鎿Q地，xs可用作電阻率結(jié)構(gòu)反演的輸入。在參數(shù)研究或反演中，ΔA和v也可用來(lái)代替xs使用。0044另一種解釋方法是在“快速”1維或2維反演中使用ΔA和v。在這樣一個(gè)實(shí)施例中，這些狀態(tài)矢量元素被阻擋從而給出分段指數(shù)幅度函數(shù)，其與簡(jiǎn)化的地質(zhì)結(jié)構(gòu)中的單層相對(duì)應(yīng)。實(shí)例0045圖5-10中圖解說(shuō)明了模型化的水平電偶極子數(shù)據(jù)的實(shí)例。圖5顯示了具有加性隨機(jī)噪聲的模型化數(shù)據(jù)。信號(hào)頻率是0.25Hz，并從零時(shí)間位置進(jìn)行大致指數(shù)衰減。加性噪聲從記錄的環(huán)境地磁噪聲中提取。圖6示出在約5000秒處，該噪聲輸入數(shù)據(jù)(初始信號(hào)加加性噪聲)61的放大圖。在同一圖上繪出的是模型數(shù)據(jù)的無(wú)噪聲信號(hào)部分62和通過(guò)使用本發(fā)明方法做出的信號(hào)的估計(jì)63。估計(jì)63近似等于無(wú)噪聲信號(hào)62，并難于在圖6中與無(wú)噪聲信號(hào)62相區(qū)分。由于本發(fā)明方法并不“知道”模型信號(hào)，因此這種緊密的相似性說(shuō)明本發(fā)明方法取得了成功。其僅被提供噪聲數(shù)據(jù)作為輸入，但能恢復(fù)信號(hào)部分。用在該例子中的本發(fā)明實(shí)施例采用Kalman濾波作為跟蹤算法，否則采用上面所述的一種作為優(yōu)選方法。相應(yīng)于圖6中所表示的部分的偏移近似為4000米。0046圖7示出無(wú)噪聲信號(hào)(實(shí)線(xiàn)71)的幅度改變率ΔA和通過(guò)使用噪聲數(shù)據(jù)經(jīng)Kalman濾波獲得的ΔA的估計(jì)值(鋸齒狀線(xiàn)72)。應(yīng)當(dāng)注意，在信號(hào)強(qiáng)度最大的數(shù)據(jù)中央部分，該估計(jì)的質(zhì)量是最好的。0047圖8示出無(wú)噪聲信號(hào)(實(shí)線(xiàn)81)相對(duì)相位改變率v和通過(guò)使用噪聲數(shù)據(jù)經(jīng)Kalman濾波獲得的v的估計(jì)值(鋸齒狀線(xiàn)82)。仍然是在信號(hào)強(qiáng)度最大的數(shù)據(jù)中央部分，該估計(jì)的質(zhì)量是最好的。0048圖9示出與Kalman估計(jì)(虛線(xiàn)92)相比的無(wú)噪聲輸入數(shù)據(jù)(實(shí)線(xiàn)91)的信號(hào)幅度。半對(duì)數(shù)圖顯示了信號(hào)強(qiáng)度中有3個(gè)數(shù)量級(jí)以上的幅度變化。圖10示出輸入信號(hào)101對(duì)Kalman估計(jì)102的類(lèi)似的相對(duì)相位比較。0049前面的描述旨在說(shuō)明本發(fā)明的具體實(shí)施例，所述具體實(shí)施例出于說(shuō)明本發(fā)明的目的。然而，對(duì)于本領(lǐng)域技術(shù)人員來(lái)說(shuō)，顯然這里描述的實(shí)施例可能有許多描述和變化。所有這些修正和變化都包含在由所附權(quán)利要求所定義的本發(fā)明的范圍內(nèi)。權(quán)利要求1.一種方法，其用于跟蹤在噪聲數(shù)據(jù)中的發(fā)射的周期性電磁信號(hào)的幅度變化和相位改變，所述噪聲數(shù)據(jù)由至少一個(gè)接收器隨時(shí)間檢測(cè)，所述信號(hào)是以已知頻率發(fā)射的，所述方法包括的步驟是a)選擇用于跟蹤已知頻率信號(hào)的跟蹤算法；b)將所述檢測(cè)時(shí)間分成時(shí)間間隔或時(shí)間段，在每個(gè)時(shí)間間隔內(nèi)檢測(cè)到的信號(hào)和至少一個(gè)相關(guān)參數(shù)被假定不變；c)估計(jì)所述檢測(cè)到的信號(hào)和所述至少一個(gè)相關(guān)參數(shù)的初始值，并將這些值分配給第一個(gè)時(shí)間間隔；d)在時(shí)間上在前一個(gè)時(shí)間間隔估計(jì)所述初始信號(hào)和每個(gè)相關(guān)參數(shù)的預(yù)測(cè)；e)使用所述數(shù)據(jù)和所述跟蹤算法來(lái)修正步驟d)中的初始估計(jì)；以及f)重復(fù)步驟d)-e)直到所有數(shù)據(jù)被處理。2.如權(quán)利要求1所述的方法，其中所述跟蹤算法是Kalman算法，其包括由狀態(tài)方程和測(cè)量方程所規(guī)定的狀態(tài)矢量，所述狀態(tài)矢量至少具有下列兩個(gè)分量所述檢測(cè)到的信號(hào)的幅度和所述至少一個(gè)所選擇的相關(guān)參數(shù)。3.如權(quán)利要求2所述的方法，其中所述狀態(tài)矢量具有兩個(gè)分量，且所述至少一個(gè)相關(guān)參數(shù)與所述信號(hào)的時(shí)間導(dǎo)數(shù)即正交信號(hào)成比例。4.如權(quán)利要求2所述的方法，其中所述狀態(tài)矢量具有四個(gè)分量所述信號(hào)的幅度；所述正交信號(hào)幅度；所述信號(hào)的包絡(luò)幅度時(shí)變率；以及所述信號(hào)的相對(duì)相位時(shí)變率，用于跟蹤信號(hào)的所述后兩個(gè)分量經(jīng)歷顯著衰減。5.如權(quán)利要求2所述的方法，其中所述Kalman算法適于在所述修正步驟中使用時(shí)間上較晚檢測(cè)的數(shù)據(jù)。6.如權(quán)利要求1所述的方法，進(jìn)一步包括初始數(shù)據(jù)縮放步驟。7.如權(quán)利要求1所述的方法，其中所述發(fā)射的信號(hào)是地下地層的電磁勘查中的源信號(hào)的傅立葉分量。8.如權(quán)利要求7所述的方法，其包括最終步驟，其用于從信號(hào)或至少一個(gè)相關(guān)參數(shù)的估計(jì)值來(lái)確定所述地下地層的電阻率結(jié)構(gòu)。9.如權(quán)利要求2所述的方法，其中在所述狀態(tài)方程中實(shí)現(xiàn)小校正線(xiàn)性化。10.如權(quán)利要求2所述的方法，其中為每個(gè)時(shí)間間隔獲取所述信號(hào)的估計(jì)及與其相聯(lián)系的誤差。11.如權(quán)利要求4所述的方法，其中通過(guò)使用小校正線(xiàn)性化假設(shè)，修正所述Kalman算法以允許將所述信號(hào)的包絡(luò)幅度的時(shí)變率A，和所述信號(hào)的相對(duì)相位的時(shí)變率v加到所述狀態(tài)矢量中，所述線(xiàn)性化假設(shè)包括(i)所述信號(hào)包絡(luò)的所述幅度被假定乘以(1+T·A)，其中A是所述信號(hào)包絡(luò)的幅度時(shí)變率，而T是時(shí)間間隔，以及(ii)所述信號(hào)假設(shè)經(jīng)歷相移vT，其中v是所述信號(hào)的相對(duì)相位改變率，并且假定vT＜＜1。全文摘要一方法，用于使用Kalman濾波或其他用來(lái)跟蹤信號(hào)幅度隨時(shí)間的變化和檢測(cè)所述信號(hào)相位隨時(shí)間的改變的跟蹤算法34，來(lái)跟蹤在噪聲地震數(shù)據(jù)中檢測(cè)的正弦電磁信號(hào)。該方法用于受控的地震源電磁勘查，其中地震源相對(duì)于接收器的較遠(yuǎn)偏移可能引起所述地震源信號(hào)顯著地衰減，并使其難于從地磁或其他電磁背景31中恢復(fù)。文檔編號(hào)G01V3/12GK1961221SQ200580018019公開(kāi)日2007年5月9日申請(qǐng)日期2005年4月26日優(yōu)先權(quán)日2004年6月1日發(fā)明者S·奧恩博斯特爾,呂新友申請(qǐng)人:?？松梨谏嫌窝芯抗?