一種大地電磁阻抗估計的時頻分析方法與流程
2023-07-15 11:44:51 3
本發明涉及大地電磁測深（Magnetotelluric，MT）資料處理領域，具體說是一種能提高大地電磁阻抗估計精度和可解釋性的時頻分析方法。
背景技術：
大地電磁測深（MT）是一種以天然交變電磁場為場源的電磁勘探方法，其電磁場信號弱、頻帶寬、極易受各種噪聲的幹擾，是典型的非線性、非平穩信號。處理MT信號最傳統的方法是Fourier變換，而Fourier變換是以平穩信號為理論基礎的，與MT信號的非平穩特性不相適應，因此用Fourier變換分析MT信號有著明顯的缺陷。在大地電磁測深法中，數據分析、定量反演和地質解釋都必須建立在正確可靠的大地電磁響應函數的基礎上才有意義，在沒有確定MT響應函數是否可靠之前，去追求定量反演的精度意義不大。怎樣獲得無偏的阻抗估算，是學者們研究MT測深的重要課題。所幸MT信號的處理技術與其它學科的發展息息相關，現代科學技術每發展一步，都可能會帶動MT處理技術前進一大步。
Hilbert-Huang變換（Hilbert-Huang transform，HHT）是最新發展起來的處理非線性、非平穩信號的時頻分析技術，主要是利用經驗模式分解( Empirical Mode Decomposition，簡稱 EMD)對Hilbert變換中的大多數非平穩信號進行平穩化的時頻處理。HHT在許多領域得到了不同程度的研究與應用，已被證明是分析非平穩信號的強有力工具，也為分析非平穩的大地電磁信號提供了一種選擇。提高大地電磁阻抗估計精度和可解釋性提供條件。
技術實現要素：
針對上述現有技術存在的不足，本發明的目的是提供一種基Hilbert-Huang變換的能提高大地電磁阻抗估計精度的的時頻分析方法。
為實現上述目的，本發明採用的技術方案如下：
一種大地電磁阻抗估計的時頻分析方法，包括如下步驟：
（1）將大地電磁信號的4個分量Hx、Ey、Hy、Ex分別都進行如下操作：採用Hilbert-Huang變換中的經驗模態分解技術，將大地電磁信號的分量進行分解，自適應地得到各階模態函數（intrinsic mode function，IMF），然後對各階IMF進行Hilbert變換，計算得到各階模態函數的瞬時幅度ai(t)和瞬時頻率ωi(t)，綜合各階IMF的瞬時譜，得到大地電磁的時頻譜H(ω, t )，這樣就可以分別得到大地電磁信號4個分量的時頻譜。
（2）根據時頻譜H(ω, t )，求得任意時頻點上信號的瞬時幅值，將得到的瞬時幅值代替傳統傅立葉方法中的傅氏譜，帶入到最小二次方程，求解時頻點上的阻抗張量，數學模型如下：
（1），
（2），
上述， ,,,是具有相同瞬時頻率的4個電磁場分量在時頻點上的瞬時幅值；和是上述4個電磁場分量中任意兩道數據的共軛複數；
（3）瞬時阻抗函數的分量根據公式（1）-（2）可求得， （3），
類同公式（3），其他幾個分量,,也同理得到；然後計算出各個時頻點的阻抗張量元素，把大地電磁信號的時頻功率譜圖轉變為時頻阻抗譜圖；時頻阻抗譜圖中，將具有相同瞬時頻率(ω0)的各時頻點的響應函數值累加起來，除以時頻點的個數，其均值則為頻率值(ω0)的響應函數估計值：
（4），
然後利用、分別求得電阻率和相位，式中ij表示下標xx, xy, yx, yy；f表示頻率。
相比現有技術而言，本發明的有益效果是：1、利用時頻分析從瞬時譜上統計估算大地電磁的響應參量，比Fourier分析更有利於實現穩健估計；2、本發明沒有要求信號頻率成分在整個時間段都平穩，而其時頻阻抗值是由真實存在的各時頻點譜值計算得到的，相比傅立葉譜，減小了信號非平穩性對譜估計帶來的誤差，更接近大地電磁信號的真實譜，由此方法計算得到的阻抗元素更符合地層結構的實際；3、不用對大地電磁信號進行窄帶濾波，瞬時頻率的頻率解析度更高；4、在時頻阻抗譜圖中，將具有相同瞬時頻率(ω0)的各時頻點的響應函數值累加起來，求其均值作為頻率值(ω0)的響應函數估計值，其過程其實質是個整體平均的過程，算法上就有抑制噪聲的能力，所以基於Hilbert-Huang的時頻阻抗估計方法比基於傅立葉譜的傳統方法有更強的噪聲抑制能力，估算的結果更具穩健性。
附圖說明
圖1為本發明工藝流程圖
圖2為本發明所述經驗模態分解法流程圖
圖3為本發明所述大地電磁4個分量Hy、Ex、Hx、Ey的EMD分解圖，其中（a）為Hy分量上的EMD分解圖，(b)為Ex分量上的EMD分解圖，(c)為Hx分量上的EMD分解圖，(d)為Ey分量上的EMD分解圖。
圖4為大地電磁4個分量的Hilbert時頻譜：其中（a）為Hy分量上的時頻譜，(b)為Ex分量上的時頻譜，(c)為Hx分量上的時頻譜，(d)為EY分量上的時頻譜。
圖5為估算的電阻率ρxy：(a)為電阻率的時頻分布圖，表徵了其時變特徵；(b)為(a)的局部放大圖(time: 0~0.0041s; frequency: 0~337.5Hz ); (c)是由(a)獲得的電阻率曲線（由具有相同瞬時頻率的各點視電阻率疊加再求平均得到）；(d)為從圖(b)中選出來的47.5Hz瞬時頻率的時變電阻率。
圖6為估算的相位Φxy： (a) 為相位的時頻分布，表徵現了其時變特徵； (b) 為(a)的局部放大圖 (time: 0~0.0041s; frequency: 0~337.5Hz ); (c) 由圖(a)獲得的相位曲線(由具有相同瞬時頻率的各點相位疊加再求平均得到)； (d) 從圖 (b)中選出來的47.5Hz瞬時頻率的時變相位。
圖7為傳統Fourier 變換方法和本發明的處理結果對比：(a) 視電阻率曲線 ρxy的對比； (b) 視電阻率曲線 ρyx 的對比；(c) 相位曲線 Φxy 的對比； (d) 相位曲線Φyx 的對比。
具體實施方式
現結合具體實施例，來對本發明作進一步描述。以下實施例僅為本發明的優選方式，並非只用於限制本發明的保護範圍，任意與本發明相同或相似的技術方案均視為本發明的保護範圍。
實施例一
如圖1所示，大地電磁阻抗估計的時頻分析方法，包括如下步驟：
（1）採用經驗模態分解法，將大地電磁信號的4個分量Hx、Ey、Hy、Ex分別進行分解。分解過程如圖2所示，分別設定大地電磁信號每個分量上的時間序列信號為，首先確定出上的所有極值點，然後將所有極值點用曲線擬合，得到信號的上、下包絡曲線，分別設為和，則上、下包絡的平均曲線為：； 用減去後，得到剩餘部分設為，即：。
理想情況下，應該是一個固有模態函數，因為的構造過程就是使它滿足固有模態函數的條件。但由於包絡曲線逼近的過衝和俯衝作用，會產生新的極值，影響原來極值的位置與大小，因此，分解得到的並不完全滿足IMF條件。雖然這種影響是間接的，因為參與篩選的是上、下包絡線的平均曲線。即使擬合過程做得非常好，信號斜坡上的一個微小凸點或凹點，也會把局部零均值從直線變為曲線時被放大為新的局部極值點，這些新的極值點正是前一次篩選過程中遺漏的。
因此，需要反覆篩選以恢復所有低幅值的疊加波。篩選過程的目的有兩個：一方面消除信號上的騎波，另一方面使波形輪廓更加對稱。用代替，用相應的上、下包絡線和重複上述過程，即：
。
上述篩分過程重複k次，當hk(t)滿足以下兩個條件：(1)整個信號中，零點數與極點數相等或至多相差1；(2)信號上任意一點，由局部極大值點確定的包絡線和由局部極小值點確定的包絡線的均值均為零，即信號關於時間軸局部對稱，則hk(t)就可定義為原始信號的第一個IMF，設為C1（t），信號的剩餘部分設為，即：
;
是從原始數據中處理得到的第1個固有模態函數，為原始信號中最短的周期分量，即頻率最高的分量。為把從原始信號中分離出來後得到的剩餘分量。 把剩餘部分作為新的原信號重複以上過程，繼續進行EMD分解，依次提取第2個固有模態函數，如此重複n次，直到所得的剩餘部分為一單調信號或前後兩個剩餘項和的標準偏差(standard deviation，SD) 滿足0.05<SD<0.3時分解完畢，SD定義如下：
此時，得到第n個固有模態函數和剩餘分量：
…
；
這樣，原始的信號可表示為所有的IMF及剩餘量之和：。大地電磁信號的4個分量Hy、Ex、Hx、Ey的EMD分解結果如圖3所示。
（2）然後對大地電磁信號的4個分量Hy、Ex、Hx、Ey分別分解得到的模態函數進行Hilbert變換：對於某一個模態函數，設為，首先求出其Hilbert變換及相應的解析信號：
；
和共同組合為一解析信號 ，採用極坐標形式表示：
，
這裡瞬時幅度； 。
解析信號的極坐標形式反映了Hilbert變換的物理含義：它是通過一正弦曲線的頻率和幅值調製獲得信號局部最佳逼近。其相應的瞬時頻率為：。對每一分量上的模態函數作上述Hilbert變換，並求解出相應解析信號的瞬時幅度ai(t)和瞬時頻率ωi(t)，從而原始信號可以表示為：，其中Re表示取實部。這裡不考慮剩餘分量，因為它是一個平均趨勢或者是一個常量。相同的數據若用Fourier形式展開是 。
（3）綜合各分量IMF的瞬時譜，得到大地電磁的時頻譜H(ω, t )，，如圖4所示。
（4）根據時頻譜H(ω, t )，求得任意時頻點上信號的瞬時幅值，將得到的瞬時幅值代替傳統傅立葉方法中的傅氏譜，帶入到最小二次方程，求解時頻點上的阻抗張量，公式如下：（1） ；
（2）；
上述， ,,,是具有相同瞬時頻率的4個電磁場分量在時頻點上的瞬時幅值；和是上述4個電磁場分量中任意兩道數據的共軛複數；是和他的共軛複數在時頻點的乘積, 定義為分量在時頻點的瞬時自功率譜。相似的， 是和在時頻點的瞬時互功率譜，對於其他分量，自功率譜和互功率譜的計算方法也是一致的。
（5）瞬時阻抗函數的分量根據公式（1）-（2）所得，（3）；
其他幾個分量,,也同理得到。根據公式（3），[微軟用戶5] 計算出各個時頻點的阻抗張量元素，把大地電磁信號的時頻功率譜圖轉變為時頻阻抗譜圖。
（6）時頻阻抗譜圖中，將具有相同瞬時頻率(ω0)的各時頻點的響應函數值累加起來，除以時頻點的個數，其均值則為頻率值(ω0)的響應函數估計值：
（4）；
然後利用、分別求得電阻率和相位，式中ij表示下標xx, xy, yx, yy；f表示頻率。
從圖5-圖6上所知，用Hilbert 瞬時譜代替傳統的Fourier譜，進行MT 響應參數的估算，將時頻功率譜轉化為時頻阻抗譜、時頻相位譜，可在時頻域上來統計估算MT響應函數。這種方法並沒有要求信號頻率成分在整個時間段都平穩，時頻響應參數是由真實存在的各時頻點譜值計算得到的。相比傅立葉譜，減小了信號非平穩性對譜估計帶來的誤差，更接近大地電磁信號的真實譜，所以由此方法計算得到的響應函數更符合實際。從HHT瞬時譜上統計估算參量，比Fourier分析，更有利於實現穩健估計, 最小化了大地電磁信號非穩性帶來的估算偏差。
從圖7可知，本發明中的阻抗時頻分析與傳統Fourier方法相比：兩種算法得到的曲線離散值的基本形態是一致的，而本發明處理結果比傳統算法的數據更集中，曲線更平穩一些。就整個曲線形態來看，規律是明顯的，估算的參量魯棒性更好，曲線異常的情況得到緩解，緩解了個別頻點魯棒性差，異常突出的現象，估算的結果更具穩健性，資料的可解釋性得到明顯提高。