當前位置:編程學習大全網 - 編程語言 - 重、磁數據常規處理方法

重、磁數據常規處理方法

所謂常規數據處理是指在重、磁數據處理中經常使用到的位場轉換或濾波處理,如上延、求導數、化極等。下面對本次研究中使用的常規數據處理方法做壹簡述。

空間域異常處理、轉換的基本公式均可以寫成如下的褶積形式

東北地球物理場與地殼演化

式中:fa、fb分別表示原始異常處理和轉換前、後的異常;?為權函數。不同的處理、轉換間只是它們的權函數不同而已。

由於上述兩式與電子電工學中描述濾波器濾波特性的卷積的形式完全相同,因此異常的處理和轉換又稱為異常的濾波。由卷積定理有

東北地球物理場與地殼演化

式中: 分別為原始異常和處理、轉換後異常的波譜; 為權函數φ的波譜(也稱為處理與轉換因子、波數響應、濾波算子)。

將(4-5)、(4-6)、(4-7)式比較可知,空間域的處理和轉換是褶積運算,而波數域是乘積運算。而且,波數譜的連乘可以完成連續的多種變換。因此數域的轉換方法要簡單得多。隨著電子計算機的廣泛應用,特別是快速傅裏葉變換算法的問世,使區域重磁資料數據處理中的波數域方法成為主要的方法。

在本次研究工作中,數據處理的計算工作是在波數域中進行的。十分明顯,為實現波數域的異常處理、轉換,必須已知或設計出處理、轉換因子 。

根據計算,重、磁異常波譜公式是壹些獨立因子的乘積,其通式為

東北地球物理場與地殼演化

式中:A為參數因子,只與地質體的剩余密度或剩余質量有關;B為參數因子,只與地質體的磁化強度或磁矩有關;H(?,h)為深度因子,僅與地質體的深度有關;S(?,a,b)為水平尺寸因子,僅與地質體的水平尺寸有關;L(ls,ms,ns,mt,nt)為方向因子,只與磁化強度和地磁場的方向有關;D(?,ξ,η)為位移因子,它是由於坐標原點任選而增加的因子。

(壹)解析延拓

根據異常波譜特征的計算,可知無限延深直角棱柱體異常波譜深度因子H=e(-h?)。若原始異常體的深度為h1,解析延拓後異常體的深度為h2,則有

東北地球物理場與地殼演化

式中:Δh=h2-h1。

於是 因子延拓因子應為

東北地球物理場與地殼演化

式中:?x、?y分別為與x軸、y軸的對應的圓波數, 為徑向圓波數。

上延時Δh>0,下延時Δh<0。上延可壓制高波數成分(即突出低波數成分),屬低通濾波;下延可放大高波數成分(即突出高波數成分),屬高通濾波,但對低波數成分無壓制作用。

(二)導數計算

重、磁異常的導數計算廣泛應用於異常的處理和解釋。原因在於:①異常的導數在不同形狀的地質體上有不同的特征,有助於對異常的解釋和分類;②異常導數可以突出反映淺部地質因素,而壓制區域深部地質因素的影響,在壹定程度上可以劃分不同深度和大小異常源產生的疊加異常;③在利用歐拉反褶積對重、磁異常進行構造反演計算時,要用到異常水平壹階導數和垂向壹階導數。

1.垂向導數

由於重、磁異常函數f(x,y,z)的n階垂向導數可以用下列公式表示

nf(x,y,z)/?n=-?nf(x,y,h)/?hn,因此由深度因子H=e(-h?)可以求出異常的n階垂向導數因子。

東北地球物理場與地殼演化

顯然壹、二階垂向導數因子分別為:

東北地球物理場與地殼演化

重、磁異常垂直導數可放大高波數成分(即突出高波數成分),但對低波數成分有壓制作用。

2.水平導數

由傅裏葉變換(FT)的微分性質可知,沿x和y方向的異常水平導數因子分別為

東北地球物理場與地殼演化

式中: 當n取1,即為壹階水平導數的轉換因子。如果異常f(x,y,h)對任意水平方向l的導數為: (其中,α為l與x軸的夾角)。依FT的微分性質可得到異常的方向導數因子。

由上式可知,異常的水平導數可突出某壹方向上異常特征(或構造線),如α=45°時,能突出135°方向的構造線。

(三)化地磁極

球體總強度磁異常T的譜Δ 為:

東北地球物理場與地殼演化

式中:qs=j(lscosθ+mssinθ)+ns,qt=j(ltcosθ+mtsinθ)+nt;而ls,ms,ns為磁化強度M的方向余弦;lt,mt,nt為地磁場T0的方向余弦;θ為徑向圓波數的極角。 為引力位的譜;G為萬有引力常數;ρ為密度;?為圓頻率;μ0為真空導有磁率。

令qt=nt=1時,即垂直磁異常的譜為:

東北地球物理場與地殼演化

再令qs=ns=1時,即化極後垂直磁異常的譜為:

東北地球物理場與地殼演化

比較以上各式,便可得到化極轉換因子為:

東北地球物理場與地殼演化

從(4-18)式看出,化極因子與?無關,因此磁異常化極無濾波作用。另外,從(4-18)式可知,磁異常化極需已知磁化強度方向。然而剩磁與感磁方向不壹致時,磁化強度的方向是難以確定的,尤其在大面積磁測資料處理時,區內磁體很多,更無法了解它們的磁化強度方向,因此往往假設磁化強度的方向與地磁場方向壹致。另外,實踐中還認為,在研究區範圍內的地磁場方向是相同的。這種假設在測區不大的情況下對結果的影響較小。圖4-2顯示出化磁時磁化傾角對結果的影響。根據這壹計算結果,若測區的南北緯度差在10°之內,用統壹的地磁場方向余弦來作化磁極運算,其結果受影響不大。因此,對壹般成礦預測為目的的區域性資料研究問題並不突出。但當作深部地質研究或大面積區域地質研究時,就要註意這種情況,有些學者已開始了測點地磁場方向余弦各異時的化磁研究工作。

圖4-2 不同磁化方向化磁極後的曲線對比

(四)譜分析方法

譜分析方法作為重、磁異常數據處理、轉換的重要方法,有著廣泛的應用。利用徑向平均對數能譜分析可以估算重、磁場源的平均深度,為進壹步的處理和解釋提供基礎信息。

下面對徑向平均對數能譜分析和平均深度估算的原理進行簡介:

經計算可知,球體重力異常的波譜為:

東北地球物理場與地殼演化

則球體重力異常功率譜為:

東北地球物理場與地殼演化

對數徑向功率頻譜為:

東北地球物理場與地殼演化

式中A=2πGρv=2πGm。上式表明,球體重力異常對數徑向功率譜與徑向圓波數中呈線性關系,見圖4-3a,故可利用lnE(?)的擬合直線斜率求解出球體中心深度:

東北地球物理場與地殼演化

式中:f為波數,而?=2πf。

另外,由球體磁異常功率譜也可計算深度。把引力位譜( =2πGme-h?/?)代入(4-15)式,經整理可求得球體垂直磁化(qs=1),垂直磁異常(qt=1)的功率譜為:

東北地球物理場與地殼演化

式中B=2πμ0VM/4π=2πμ0m/4π;m為磁矩;V為球體體積;M為磁化強度。

東北地球物理場與地殼演化

上式說明,球體垂直磁化垂直磁異常的對數徑向功率譜與圓波數呈非線性關系(圖4-3b)。但是,在高波數段近於線性關系,可用(4-21)式計算深度。

圖4-3 對數徑向功率譜

(五)重、磁對應分析

基於泊松定理發展起來的重磁異常對應分析方法,是重磁數據綜合解釋的重要方法,能對重磁異常的相關性進行定量研究,有效地將重磁信息進行綜合,對重磁資料定性地賦予地質意義,並突出地質目標的反映,為重磁資料的地質解釋提供有用的信息,特別是在強磁性火山巖解釋中具有重要的作用。

重磁異常對應分析方法的基本理論如下:

由同壹場源引起的重力異常和磁異常間的關系可以簡單地用泊松方程描述。當垂直磁化時,泊松方程可表示為:

東北地球物理場與地殼演化

式中:Δz⊥為垂直磁化的垂直磁異常;M為場源磁化強度;G為萬有引力常數;Δρ為場源剩余密度;Δg為重力異常; 為重力異常的垂向壹階導數。

上式表明垂直磁化的垂直磁異常與重力場的垂向壹階導數滿足線性關系,而且擬合直線的截距為零。

由於原始資料不可避免地存在某些幹擾因素,通常進行重磁異常的線性回歸分析時,選用如下稍加推廣的泊松方程:

東北地球物理場與地殼演化

式中:b為斜率;A為截距。

將Δz⊥與 作線性回歸分析則可得到斜率b與截距A的估計值。兩個離散序列的相關導數可以由下式求得:

東北地球物理場與地殼演化

式中:Cxy(k) 為兩個離散序列x(t)={x1,x2,…,xn}和y(t)={y1,y2,…,yn}的相關函數,k為延遲時間。當x(t)=y(t)時,稱為自相關函數Cxx(k)或Cyy(k)。

計算處理時,給定適當大小的分析窗口,將窗口內各點垂直磁化磁異常和重力異常的垂向壹階導數進行最小二乘線性回歸,求得中心點的相關系數R、斜率b和截距A。

相關系數R反映了在給定窗口內重磁異常的線性相關程度,即宏觀地反映了重磁異常的“同源性程度”。相關系數絕對值接近於1的窗口區間重磁異常的“同源性”較好,它們或者同源、或者都離場源較遠、或者同處異常的拐點等。其中R接近+1時,重磁異常正相關;R接近-1時,重磁異常負相關。當R絕對值較小時,重磁異常相關性差,重磁異常可能不同源,或存在鄰近異常幹擾,或是存在方向不同於地磁場的強剩磁磁性體等。

斜率b反映了所有場源泊松比的加權平均值,稱為廣義泊松比。只有在重磁異常同源的前提下,回歸所得的斜率b才有意義。僅由b不能直接確定M和Δσ,但若在解釋中結合其他地質、地球物理信息,就能從中獲得關於物性分布的有用信息,從而為進壹步的定量解釋提供依據。

截距A反映了實測資料中的長波長成分,它主要反映重磁異常數據的背景變化。在重磁異常完全同源的理想情況下,A=0。

由於重磁異常對應分析是對場源之間的相關系數進行定性和半定量研究的方法,它能分離和鑒別不同類型的異常,從而勾畫出與異常場源相對壹致的地質單元和構造分區,不相關區說明重磁異常不同源或存在鄰近異常幹擾。

(六)歐拉反褶積與構造反演

歐拉反褶積方法使用歐拉(Euler)齊次關系,對經方向譜分析過的數據快速估計重、磁場源的位置和深度,是壹種既能夠利用重磁網格數據,又對剖面數據有效地確定地質體位置(邊界)和深度的定量反演方法(Reid等,1990)。這種方法並不需要已知地質信息(密度、磁化率等)的控制。使用該方法可以將位場及其梯度以及場源位置之間的關系用歐拉齊次方程表示,而場源的不同形狀即地質構造的差異則表現為方程的齊次程度,就是所謂的地質構造指數,地質構造指數或齊次程度實質上表現了場隨離開場源距離的衰減率。模型研究和應用實例表明,這個方法對確定斷層、磁性接觸帶、巖脈、噴出巖體等構造位置或勾繪它們的輪廓有較高的精度。

位場的歐拉方程是由Thompson推導的。首先建立壹個直角坐標系,取觀測平面為z=0,z軸向下為正,x軸指北,y軸指東。考慮在此坐標系中的任壹函數f(x,y,z),如果

東北地球物理場與地殼演化

則稱函數是n階齊次的。此外可證明,如果f(x,y,z)是n階齊次的,則滿足下列方程

東北地球物理場與地殼演化

此偏微分方程稱為歐拉齊次方程,或稱歐拉方程。

對於位於(x0,y0,z0)的點磁源,在觀測平面上任壹點(x,y,z)處的總磁場強度具有如下形式:

東北地球物理場與地殼演化

式中 N=1,2,3,……。G不依賴於(x,y,z)。對於(3-23)這樣的函數,其歐拉方程可寫成

東北地球物理場與地殼演化

方程(4-29)是n=-N階齊次的。三個坐標方向的梯度值可以利用空間域或波數域的壹般位場變換計算出來。如果梯度值通過觀測獲得,直接用於方程(4-29)則更可取。

方程(4-29)雖然是根據磁源異常推導的,但對重力異常也同樣適用。該方程用於平面網格重、磁異常數據的反演計算。如果假定方程中橫向梯度?ΔT/?y為零,則可得到適用於剖面數據計算的方程。這對於眾多走向方向不變的二度情況很顯然就是這樣。

齊次度N被定義為“構造指數”,它是重、磁異常場源深度變化“陡緩”的量度。特定的地質構造具有特定的衰減率(即:構造指數)。例如:傾斜斷層的磁場、水平薄巖脈的磁場按線性的規律變化,構造指數就為1。表4-1列出了構造指數對應的地質構造。

利用不同坐標點(x,y,z)上的場值ΔT及其三個方向上的梯度值以及方程(4-29)組成的線性方程組,最後可以解出未知變量x0、y0、z0,進而確定構造形跡及位置。

表4-1 歐拉構造指數表

但是,直接用方程(4-29)及其變換的二度形式解決構造問題,會使解的精度極不可靠和不穩定。主要原因有如下幾個方面:

(1)很難知道磁場ΔT的絕對水平,區域場或鄰近磁異常的影響幾乎總是存在的。

(2)根據線性方程組與系數的關系,較低的構造指數才會有較好的深度估計值。但大多數磁異常是偶極性的,有較高的構造指數。同時又有許多線性構造的指數接近於零而使反褶積發散。

(3)實測異常是多種構造指數特征的復雜疊加,很難用壹些簡單模型來模擬,亦很難將具有線性特征的構造識別與分離出來。

為克服上述三個方面的問題,釆用下面的壹些辦法:

從觀測數據中消除偏差是通過網格數據進行窗口計算解決的。對網格數據假定異常在方程(4-29)求值的窗口範圍內有壹常量偏差,觀測值為

東北地球物理場與地殼演化

這裏B是常數。從方程(4-30)中解出ΔT,代入(4-29),整理得

東北地球物理場與地殼演化

如果構造指數小於0.5,即構造指數接近零時,這樣可能造成對深度值的過低估計。為此需要提供壹個補償值A,使得歐拉齊次方程在構造指數較低時寫為

東北地球物理場與地殼演化

式中:A是與場幅值有關的壹個參數。不同的構造形體有不同的A,A可以通過將已知的某壹構造的解析式代入歐拉方程(4-29)而求出。

圖4-4 重力異常歐拉反褶積計算結果示意圖

圖4-4中的曲線為重力異常等值線,圓圈為反演解的構造位置。圓圈直徑的大小代表了不同的構造深度。

(七)重、磁人機交互剖面正反演

該項技術的優點是便於將重、磁異常的處理、轉換方法得出的結果和其他地質、地球物理方法獲取的先驗信息輸入到模型裏,形成初始模型。並且根據計算結果和實際重、磁異常的差異,隨時方便地修改模型,直觀地監督和指導正反演過程。重、磁人機交互剖面正反演流程見圖4-5。

圖4-5 人機交互正反演流程圖

1.重力人機交互正反演技術

重力人機交互正反演技術(Gamble,1979)主要是依據A截面為多邊形的二度體重力異常計算方法來實現的。通過對初始模型計算出的重力效應與測線上的布格重力異常進行對比,不斷修正模型,直至達到計算出的重力效應與測線上的布格重力異常之差滿足預定精度。重力人機交互正反演流程見圖4-5。

圖4-6 二度體A截面

A截面為多邊形的二度體重力異常計算方法:

假設二度體的剩余密度為σ,以計算點作為坐標原點,x軸與二度體走向垂直,z軸鉛垂向下(圖4-6)。若n邊形第k個頂點的坐標為(ξk,ζk),其中k=1,2,...,n。則(ξk,ζk)與(ξk+1,ζk+1)兩個頂點連線上ξ與ζ有如下關系:

東北地球物理場與地殼演化

引用解Δg正問題的基本公式,首先對ξ求積分,得

東北地球物理場與地殼演化

式中:s為多邊形的A截面積;l為A截面的周長。

將式(4-33)代入(4-34)得

東北地球物理場與地殼演化

對上式積分可得到如下結果

東北地球物理場與地殼演化

或寫成下面的形式

東北地球物理場與地殼演化

(4-36)式或(4-37)可以編成計算機程序,用以計算A截面為任意多邊形的二度體的重力異常,進而可以進行重力人機交互正反演。在具體編程計算時應註意以下幾個問題:

(1)因多邊形的邊數為n,故ξn+1=ξ1,ζn+1=ζ1;

(2)(4-36)和(4-37)兩式是假設計算點位於原點時導出的,因此,當任意計算點P(x,y)的重力異常時,式中的ξk和ζk應以ξk-x和ζk-z來代替;

(3)在(4-36)和(4-37)式中,反正切函數的取值範圍應在-π到π之間,即當ξk+1>ξk時,反正切函數在0到π之間取值;反之,則在-π到0之間取值。

2.磁法人機交互正反演技術

磁法人機交互正反演技術主要是依據A截面為多邊形的二度體磁力異常計算方法來實現的。其基本思想同重力人機交互正反演技術相壹致。

由於V2=Δg,所以根據公式(4-37)可以求出引力位的二階導數

東北地球物理場與地殼演化

將(4-38)和(4-39)式代入下面二度體的磁異常公式,就可以利用該式進行磁力人機交互正反演計算。

東北地球物理場與地殼演化

式中:μ0為真空的磁導率;MS為有效磁化強度;is為有效磁化傾角;I0為地磁場傾角;A'為x軸與磁北的夾角。在具體編程序上機計算時應註意的問題等方面與重力人機交互方法相同(圖4-5)。

  • 上一篇:C C++數學的實際用途
  • 下一篇:如何用單片機壹個按鈕開關控制3個燈3種效果。流水。閃爍。還有壹直亮。用keil編程。
  • copyright 2024編程學習大全網