耿美霞,黃大年,楊慶節(jié)
吉林大學地球探測科學與技術(shù)學院,長春 130026
改進的TRUST方法在航磁數(shù)據(jù)線性特征增強中的應用
耿美霞,黃大年,楊慶節(jié)
吉林大學地球探測科學與技術(shù)學院,長春 130026
基于各向異性擴散濾波的TRUST方法對航磁數(shù)據(jù)線性特征趨勢增強時,未考慮數(shù)據(jù)光滑區(qū)域與異常邊緣區(qū)域的區(qū)別,在光滑區(qū)域擴散持續(xù)進行,導致光滑區(qū)域產(chǎn)生許多虛假小異常。為此,提出一種改進的TRUST方法。該方法對TRUST方法中擴散張量的特征值通過引入相干性的辦法進行了重新構(gòu)造,使其能夠正確識別光滑區(qū)域和異常邊緣區(qū)域,并自動降低光滑區(qū)域的擴散速率,從而避免濾波過程中虛假異常的產(chǎn)生。在模型試驗中,將該方法與TRUST方法處理結(jié)果進行了對比。結(jié)果顯示,采用改進后的TRUST方法,不僅使航磁數(shù)據(jù)圖上“串珠狀”走樣現(xiàn)象被消除,而且光滑區(qū)域的低振幅信息也被很好地保存。最后應用此方法對某地區(qū)的航磁數(shù)據(jù)進行趨勢增強濾波,證實了該方法的實用性。
航磁數(shù)據(jù);串珠狀走樣;趨勢增強;TRUST方法
地球物理數(shù)據(jù)處理往往要求數(shù)據(jù)是間隔均勻的網(wǎng)格數(shù)據(jù),因此需要采用合理的網(wǎng)格化插值方法對原始測量數(shù)據(jù)進行網(wǎng)格化[1]。網(wǎng)格數(shù)據(jù)的質(zhì)量直接影響后期數(shù)據(jù)處理以及解釋圖件的質(zhì)量和可靠性。然而在航空數(shù)據(jù)采集中,測線之間的距離往往比測點之間的距離大很多,使得測量數(shù)據(jù)的不均勻性非常嚴重。這種密度顯著不均勻的數(shù)據(jù)被稱為測線型數(shù)據(jù),對測線型數(shù)據(jù)網(wǎng)格化,常常會帶來假異常[1]。例如,當線性異常的趨勢方向與測線方向的夾角為銳角時,網(wǎng)格化后會出現(xiàn)“串珠狀”假異常,并且夾角越小,這種假異?,F(xiàn)象越明顯。“串珠狀”假異常不僅降低了航磁數(shù)據(jù)圖形光滑度,還使本應明顯的線性異常特征變成了小異常體特征,給后續(xù)處理和解釋工作帶來不便,甚至導致錯誤的解釋。因此,為了滿足當前高精度航空物探數(shù)據(jù)處理和解釋的需要,保持原始測量信息的同時需要消除或減弱“串珠狀”假異?,F(xiàn)象,這已成為一個不能回避的問題。
大量研究[2-5]表明,無論是采用最小曲率、Akima插值還是雙三次樣條插值等方法,都會產(chǎn)生“串珠狀”假異?,F(xiàn)象,即“串珠狀”假異常與常規(guī)網(wǎng)格化方法無關(guān),而是由于測線間采樣不足導致的短波長特征空間混疊效應而產(chǎn)生的。往往當測線間距大于或等于異常寬度的一半時,就會出現(xiàn)這種現(xiàn)象,并且異常走向與測線夾角越小,“串珠狀”假異常越明顯。前人研究[6-9]表明,使用最小曲率網(wǎng)格化過程中,用實際測量梯度數(shù)據(jù)代替計算得到的梯度數(shù)據(jù),在一定程度上能夠增強航磁數(shù)據(jù)中的線性特征,提高異常分辨率。但由于測量設備的限制,梯度數(shù)據(jù)并非總能獲得。
由于“串珠狀”走樣問題的廣泛存在,已經(jīng)有很多學者對此進行了研究。Zhou[10]使用Radon變換的方法進行趨勢增強,但是該方法只能對一個方向進行趨勢增強。Sykes等[11]采用2-D Radon變換,通過在變換域中取閥值,對大于閥值的變換值作一定的增長,這樣反演重建后的圖像中相應的線性特征得到增強。但由于重磁圖像線性特征的灰度相比于“背景”值并無特殊之處,閥值的選取較為困難。另外,該方法不能準確恢復原始圖像的振幅,限制了其在高精度航磁數(shù)據(jù)處理中的應用。Hansen[12]采用各向異性克里金方法對數(shù)據(jù)進行網(wǎng)格化,由于該方法只能對一個固定的方向進行趨勢增強,尚無法在有多個方向的趨勢數(shù)據(jù)中應用。Fitzgerald等[13]沿著異常值方差最小的方向,從測線上抽取數(shù)據(jù)值,并將這些抽取的數(shù)據(jù)作為原始觀測數(shù)據(jù)插入到測線之間,隨后進行樣條或最小曲率網(wǎng)格化;該方法需要確定一個接受度來判定特征是否為異常,網(wǎng)格化結(jié)果受接受度的影響非常大,并且對于走向與測線夾角小于45°的異常,處理效果不是很好。此外,Keating[14]采用最鄰近分析法找出測線上的最大值和最小值,通過連接相鄰測線上最鄰近的最大值(或最小值)來確定插入測線間的數(shù)據(jù),然后再進行網(wǎng)格化;該方法在線性特征走向與測線夾角比較大時能取得不錯的趨勢增強效果,但在構(gòu)造復雜地區(qū)容易產(chǎn)生錯誤的線性趨勢。
此外,基于偏微分方程尤其是各向異性擴散濾波方程的方法也已被應用到航磁數(shù)據(jù)趨勢增強領(lǐng)域中[1, 15-16]。Smith等[16]成功地將各向異性擴散濾波方法應用到航磁數(shù)據(jù)趨勢增強領(lǐng)域。該方法被稱為TRUST(trends using structure)方法,它是在網(wǎng)格化后的數(shù)據(jù)上進行的,能自動對多個走向的趨勢特征進行增強。然而TRUST方法未考慮數(shù)據(jù)光滑區(qū)域與異常邊緣區(qū)域的區(qū)別,擴散在光滑區(qū)域持續(xù)進行,因而在光滑區(qū)域不可避免地產(chǎn)生虛假小異常。這種假異常在解析信號模圖上會更加明顯[16],需要后續(xù)處理來消除。
筆者針對TRUST方法進行趨勢增強時在光滑區(qū)域中產(chǎn)生虛假異常的缺點,基于相干性[17]重新構(gòu)造新的特征值,從而構(gòu)造新的擴散張量。新的特征值能夠正確判斷航磁異常數(shù)據(jù)中的光滑區(qū)域和異常邊緣區(qū)域,以實現(xiàn)對不同的異常特征進行不同程度的擴散濾波,從而在消除“串珠狀”構(gòu)造的同時保護光滑區(qū)域的低振幅信息。
基于偏微分方程的擴散濾波理論最初由Perona等[18]提出,并在數(shù)字圖像處理領(lǐng)域中得到廣泛應用。其構(gòu)建的非線性擴散濾波模型為
其中:t為擴散時間;div為散度算子;g(·)為擴散濾波器的擴散系數(shù),為一非負單調(diào)遞減函數(shù),且滿足g(0)=1,g(∞)=1;u為t時刻的擴散濾波結(jié)果;u0表示t=0時的原始數(shù)據(jù),作為該擴散方程的初始條件。
為了滿足增強圖像中各向異性紋理結(jié)構(gòu)的需要,Weickert[19]通過對圖像數(shù)據(jù)中各向異性結(jié)構(gòu)的分析,將式(1)中的擴散系數(shù)變?yōu)閺埩?,并?gòu)建了如下的張量型擴散模型[15]:
J0(uσ)=uσuσT=
Jρ(uσ)=Gρ*(uσ
第一個特征向量v1=(cosφ,sinφ)T滿足
‖
各向異性擴散張量D與結(jié)構(gòu)張量J的特征向量相同,在Smith等[16]提出的TRUST方法中,其特征值為
式中,λ1和λ2又被稱為擴散速率。D的分量分別為
在實際應用中,式(2)可以用有限差分方法求解,空間導數(shù)通常由中心差分代替[15]。式(2)的有限差分形式為
un+1=un+Δtdiv(Dn
其中:Δt為空間步長,通常情況下,只要Δt<0.25,式(11)在一定迭代次數(shù)內(nèi)就是穩(wěn)定的;un是u(x,y)在位置(x,y)處nΔt時刻的近似。由于擴散張量D在每次迭代過程中變化非常小,因此筆者選擇在經(jīng)過一定迭代次數(shù)后再重新計算一次D,這樣可以大大提高該方法的計算效率。
觀察式(9)可發(fā)現(xiàn),特征值λ2沒有考慮到光滑區(qū)域與邊緣區(qū)域之間的差別,在相對光滑的區(qū)域濾波容易導致虛假異常的產(chǎn)生:在光滑區(qū)域,只要存在梯度各向異性,就會沿著梯度值較小的方向擴散,從而導致虛假異常的出現(xiàn)。
實際上,結(jié)構(gòu)張量的特征值μ1和μ2包含了圖像結(jié)構(gòu)信息[17]。所以特征值μ1、μ2可以用來描述圖像的局部結(jié)構(gòu)特征:在光滑區(qū)域μ1≈μ2≈0;對于構(gòu)造的邊緣區(qū)域μ1?μ2≈0;在拐角的邊界區(qū)域有μ1>μ2?0。由此定義相干性H[17]如下:
根據(jù)局部圖像的相干性,筆者構(gòu)造了新的特征值:
其中,C為常數(shù)。
擴散濾波的過程是將原始的圖像數(shù)據(jù),即航磁數(shù)據(jù)進行濾波。其在消除“串珠狀”假異常的同時,不可避免地造成原始測點上數(shù)據(jù)的改動。這樣的濾波處理違背了高精度航磁異常處理的原則[16]。一個防止原始測線數(shù)據(jù)被改動的簡單方法,就是在迭代濾波過程中將這些測點數(shù)據(jù)固定不變。在實際應用中,另一個問題是迭代停止時間的選擇。迭代次數(shù)過少,圖像上的線性特征得不到足夠的增強;而迭代次數(shù)過多不僅花費的計算時間較長,還會導致圖像被過度平滑。因此,選擇適當?shù)牡螖?shù)對最終趨勢增強效果非常重要。根據(jù)筆者的試驗結(jié)果,一般情況下,迭代次數(shù)取300~500,就能得到比較滿意的趨勢增強效果。
在結(jié)構(gòu)張量構(gòu)造中有2個參數(shù),σ和ρ。其中,參數(shù)σ被稱為噪聲尺度[16]。對原始數(shù)據(jù)處理前,先用噪聲尺度為σ的高斯核對數(shù)據(jù)進行卷積,可以保證在進行求導計算時,尺度小于σ的數(shù)據(jù)不對導數(shù)計算產(chǎn)生影響,從而保證了求導過程的穩(wěn)定性和可靠性。在選取參數(shù)σ時,根據(jù)航磁數(shù)據(jù)的噪聲水平進行選擇,才能達到預期的目的。根據(jù)筆者的試驗結(jié)果,發(fā)現(xiàn)通常σ選擇為1倍測線距左右時,能夠得到較好的濾波效果。
參數(shù)ρ被稱為綜合尺度,它反映的是圖像的紋理特征[16]。它的作用是對結(jié)構(gòu)張量的方向進行濾波,進一步消除尺度規(guī)格小于ρ的噪聲干擾。如果ρ的值過大,就會使濾波后的數(shù)據(jù)過度平滑,丟失原有信息。根據(jù)筆者的試驗結(jié)果,通常ρ選擇為1.5倍的測線距。這樣能夠保證連續(xù)出現(xiàn)在相鄰2條測線上的“串珠”異常得到增強。
a.濾波前;b.TRUST方法濾波后;c.改進的TRUST方法濾波后。圖1 趨勢增強前后航磁異常(ΔT)圖Fig.1 Application of the trend enhancement filter to the aeromagnetic data
a.原始航測數(shù)據(jù)網(wǎng)格化后異常;b.原始異常解析信號模;c.經(jīng)過改進的TRUST方法濾波后的航磁異常;d.經(jīng)過改進的TRUST方法濾波后異常的解析信號模。AS=。圖2 航測數(shù)據(jù)趨勢增強前后異常和相應的解析信號模圖Fig.2 Comparison of the original aeromagnetic data and the data after trend enhancement as well as the corresponding analytic signal maps
為了驗證改進的TRUST方法消除“串珠狀”假異常的有效性,筆者采用了含有高斯白噪聲的組合磁性體模型對趨勢增強效果進行檢驗。組合模型包括1個塊狀和3個不同走向的脈狀體,磁化方向I0=75°,磁偏角A0=25°,磁化強度J=1 A/m。測線方向為南北方向。如圖1a所示,3個脈狀體與測線的夾角分別為30°,60°和90°。測線間距為 200 m,測線上測點間距為 10 m,在原始的測線數(shù)據(jù)中加入了8%的高斯白噪聲。采用最小曲率方法對測線數(shù)據(jù)進行網(wǎng)格化,網(wǎng)格間距為 50 m。圖1a 給出了網(wǎng)格化后的航磁異常。從圖1a 可以看到:當線性磁異常與測線夾角為銳角時,網(wǎng)格化后異常呈現(xiàn)出“串珠狀”假異常特征,并且,夾角越小,這種走樣現(xiàn)象越嚴重;垂直于測線的線性磁異常以及塊狀磁異常,雖然沒有出現(xiàn)走樣現(xiàn)象,但是網(wǎng)格化后異常的邊緣呈“鋸齒狀”,影響了圖像的光滑性。
圖1b顯示的是經(jīng)過TRUST方法擴散濾波后的航磁異常圖。迭代次數(shù)為300。可以看到:經(jīng)過TRUST方法濾波后,與測線夾角為銳角的線性趨勢特征得到增強,“串珠狀”假異常基本被消除;同時,塊狀異常和與測線垂直的線性異常邊緣的“鋸齒”特征也被消除。然而,由于噪聲的存在,濾波后光滑區(qū)域產(chǎn)生了大量小的假異常。此外沿測線方向也有少量的高頻假異常產(chǎn)生。
圖1c顯示的是經(jīng)過改進的TRUST方法擴散濾波后的航磁異常。迭代次數(shù)同樣為300??梢钥吹剑航?jīng)過改進的TRUST方法濾波后,與測線夾角為銳角的線性趨勢特征得到增強,“串珠狀”假異常被完全消除;塊狀異常和與測線垂直的線性異常邊緣的“鋸齒”特征也被很好地消除。此外,由于改進的TRUST方法中擴散張量的特征值能夠在相對平滑區(qū)域自動降低擴散速率,有效地保護了原始低振幅信息,從而避免了平滑區(qū)虛假小異常的產(chǎn)生。
圖2a 為某地區(qū)采集的航空磁異常數(shù)據(jù),該地區(qū)存在大量的西北--東南走向的脈狀磁性體。飛行測線距為 200 m,測點距為 17 m,測線方向為東西方向。采用最小曲率法對原始的航磁數(shù)據(jù)進行網(wǎng)格化,網(wǎng)格間距為 40 m。可以看到,圖2a上存在明顯的“串珠狀”走樣,并且與測線方向夾角越小,“串珠狀”走樣越明顯。圖2b 為該異常的解析信號模。解析信號模類似于航磁化極圖,能夠部分消除斜磁化影響,異常幅度大體反映了場源的相對磁性強度[20],并且對異常體的邊緣非常敏感,因此常被用來估計磁異常體的水平位置。在解析信號模(圖2b)上,“串珠狀”走樣特征更加明顯,給后續(xù)處理和解釋工作帶來不便,甚至可能導致錯誤的解釋。
圖2c 為圖2a 中的數(shù)據(jù)經(jīng)過改進后的TRUST方法濾波后的航磁異常圖。迭代次數(shù)為300??梢钥吹剑?jīng)過濾波之后,“串珠狀”假異常被消除,同時線性特征也得到了增強。由于改進后的TRUST方法能夠在光滑區(qū)域自動降低擴散速率,因此光滑區(qū)域的低振幅信息也被很好地保存下來。圖2d 為圖2c 的解析信號模,可以看到,“串珠狀”走樣特征也已被消除,使后期解釋更加方便。
為了檢驗趨勢增強前后航磁圖上數(shù)據(jù)的變化,從圖2a 上用沿某一線性趨勢方向虛線標出的點提取數(shù)值,并將抽取的數(shù)據(jù)點構(gòu)成一條剖面(圖3)。經(jīng)過300次迭代濾波后,該剖面數(shù)據(jù)再次顯示在圖3中??梢钥吹?,濾波后剖面數(shù)據(jù)上的波峰抖動被消除,得到了一個相對光滑的長波長曲線。圖3中五角星指示的峰值都是航空測線穿越的位置。由于在進行濾波運算的同時固定原始測點數(shù)據(jù)不變,所以濾波后原始測點上的數(shù)據(jù)被保存下來。
圖3 濾波前后抽取的剖面數(shù)據(jù)對比Fig.3 Comparison of the original data on the ridge profile and the data after trend enhance
本文針對TRUST方法擴散張量特征值的選取不適合光滑區(qū)域濾波而導致虛假異常的缺點,引入圖像相干性,重新構(gòu)造了擴散張量的特征值。新的特征值能夠正確區(qū)分航磁異常中的異常邊緣和光滑區(qū)域,并在光滑區(qū)域自動降低擴散速率以便有效地保護低振幅信息,從而避免了假異常的產(chǎn)生。此外,為了保護原始測信息,本文算法在濾波過程中對原始測線上的數(shù)據(jù)點進行了約束保護,使得處理后的數(shù)據(jù)精度得到很大提高,在實際數(shù)據(jù)應用中取得了非常好的效果。
[1] 姚長利,龐旭林. 航磁數(shù)據(jù)線性特征增強濾波方法技術(shù)研究[C]//中國地球物理年會.北京:中國地球物理學會,2009:267. Yao Changli, Pang Xulin. Researeh on the Enhancement Trends of Magnetic Date[C]// Chinese Geophysics Society Annual Meeting. Beijing: Chinese Geophysical Society,2009: 267.
[2] Billings S, Richards D. Quality Control of Gridded Aeromagnetic Data[J]. Exploration Geophysics, 2000, 31 (4): 611-616.
[3] 李麗麗,杜曉娟,馬國慶. 改進的局部波數(shù)法及其在磁場數(shù)據(jù)解釋中的應用[J]. 吉林大學學報:地球科學版,2012,42(4): 1179-1185. Li Lili, Du Xiaojuan, Ma Guoqing. Improved Local Wave Number Methods in Interpretation of Magnetic Fields[J]. Journal of Jilin University: Earth Science Edition, 2012, 42(4): 1179-1185.
[4] 曾昭發(fā),吳燕岡,郝立波,等. 基于泊松定理的重磁異常分析方法及應用[J]. 吉林大學學報:地球科學版,2006,36(2): 279-283. Zeng Zhaofa, Wu Yangang, Hao Libo, et al. The Poisson’s Theorem Based Analysis Method and Application of Magnetic and Gravity Anomalies[J]. Journal of Jilin University: Earth Science Edition, 2006, 36(2): 279-283.
[5] 方東紅,曾昭發(fā),陳家林. 基于小波分析的重磁數(shù)據(jù)求導方法及應用[J]. 吉林大學學報:地球科學版,2008,38(6):1049-1054. Fang Donghong, Zeng Zhaofa, Chen Jialin. The Derivatives Calculation Based on Wavelet Analysis of G/M Data and Its Application[J]. Journal of Jilin University: Earth Science Edition, 2008, 38(6): 1049-1054.
[6] Marcotte D L, Hardwick C D, Lemieux J M, et al. Aeromagnetic Gradiometry Methods: A Study Using Real Data[C]//60th Annual International Meeting. SEG: [s.n.], 1990:584-586.
[7] Cowan D R, Baigent M, Cowan S. Aeromagnetic Gradiometers: A Perspective[J]. Exploration Geophysics, 1995, 26: 241-246.
[8] McMullan S R, McLellan W H, Koosimile D I. Three Dimensional Aeromagnetics[J]. Preview, 1995, 57: 83-91.
[9] Hardwick C. Gradient-Enhanced Total Feld Gridding[C]// 69th Annual International Meeting.SEG: [s.n.], 1999:381-384.
[10] Zhou Y X. Radon Transform Application to the Im-proved Dridding of Airborne Geophysical Survey Data[J]. Geophysical Prospecting, 1993, 41: 459-494.
[11] Sykes M P, Das U C. Directional Filtering for Linear Feature Enhancement in Geophysical Maps[J]. Geophysics, 2000, 65: 1758-1768.
[12] Hansen R O. Interpretive Gridding by Anisotropic Kriging[J]. Geophysics, 1993, 58: 1491-1497.
[13] Fitzgerald D, Yassi N, Dart P. A Case Study on Geophysical Gridding Techniques: Intrepid Perspective[J]. Exploration Geophysics, 1997, 28: 204-208.
[14] Keating P. Automated Trend Reinforcement of Aero-magnetic Data[J]. Geophysical Prospecting, 1997, 45: 521-534.
[15] 孫夕平,杜世通,湯磊. 相干增強各向異性擴散濾波技術(shù)[J].石油地球物理勘探,2004,39(6):651-655. Sun Xiping, Du Shitong, Tang Lei. Coherent-Enhancing Anisotropic Diffusion Filtering Technique[J]. Oil Geophysical Prospecting, 2004, 39(6): 651-655.
[16] Smith R S, O’Connell M D. Interpolation and Grid-ding of Aliased Geophysical Data Using Constrained Anisotropic Diffusion to Enhance Trends[J]. Geophysics, 2005, 70(5): 121-127.
[17] 王大凱,侯榆青,彭進業(yè). 圖像處理的偏微分方程方法[M]. 北京:科學出版社,2008. Wang Dakai, Hou Yuqing, Peng Jinye. Partial Differential Equation in Image Processing[M]. Beijing: Science Press, 2008.
[18] Perona P, Malik J. Scale Space and Edge Detection Using Anisotropic Diffusion[J]. IEEE Trans Pattern Anal Mach Intell, 1990, 12: 629-639.
[19] Weickert J. Anisotropic Diffusion in Image Pro-cessing[M]. Stuttgart: Teubner,1998.
[20] 駱遙,王明,羅鋒,等. 重磁場二維希爾伯特變換:直接解析信號解釋方法[J]. 地球物理學報,2011,54(7):1912-1920. Luo Yao, Wang Ming, Luo Feng,et al. Direct Analytic Signal Interpretation of Potential Field Data Using 2-D Hilbert Transform[J]. Chinese Journal of Geophysics, 2011, 54(7): 1912-1920.
Application of Improved TRUST Method in Enhancing Linear Trends of Aeromagnetic Data
Geng Meixia,Huang Danian,Yang Qingjie
College of GeoExploration Science and Technology, Jilin University, Changchun 130026, China
As for the enhancement of linear Trends, the TRUST method proposed by Smith and O’Connell does not consider the distinctions between the smooth area and other image features. The diffusion is carried on in smooth area and thus inevitably produces artifacts. An improved TRUST method was proposed. The eigenvalues of the structure matrix are reconstructed,so that the method only diffuses the image in specific areas where strong anisotropy is detected. This method is tested on both synthetic and field data sets. The results indicate that the method can be used efficiently to enhance linear structure locally in multiple directions. We have also compared the proposed method with the TRUST method by applying them to the synthetic data set. The results demonstrate that the proposed method can produce better effect with fewer artifacts.
aeromagnetic data;boudinage;enhancement of trends;TRUST method
2014-03-04
國家“863”計劃項目(2014AA06A613);吉林大學研究生創(chuàng)新基金資助項目(2014066)
耿美霞(1986--),女,博士研究生,主要從事航空重磁數(shù)據(jù)處理和反演方面的研究,E-mail:gengmeixia@126.com
黃大年(1958--),男,教授,博士生導師,國家千人計劃特聘專家,主要從事移動平臺探測數(shù)據(jù)處理與解釋及一體化軟件平臺開發(fā),E-mail:dnhuang@jlu.edu.cn。
10.13278/j.cnki.jjuese.201404301.
10.13278/j.cnki.jjuese.201404301
P631.2
A
耿美霞,黃大年,楊慶節(jié).改進的TRUST方法在航磁數(shù)據(jù)線性特征增強中的應用.吉林大學學報:地球科學版,2014,44(4):1333-1339.
Geng Meixia,Huang Danian,Yang Qingjie.Application of Improved TRUST Method in Enhancing Linear Trends of Aeromagnetic Data.Journal of Jilin University:Earth Science Edition,2014,44(4):1333-1339.doi:10.13278/j.cnki.jjuese.201404301.