用于地殼弧形構(gòu)造信息提取的四階譜矩分析

付麗華 曾 誠 楊文采 李燕梅

（①中國地質(zhì)大學(xué)（武漢）數(shù)學(xué)與物理學(xué)院，湖北武漢430074；②西安電子科技大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，陜西西安710071；③浙江大學(xué)地球科學(xué)學(xué)院，浙江杭州310027；④哈爾濱工業(yè)大學(xué)數(shù)學(xué)學(xué)院，黑龍江哈爾濱150001）

0 引言

在位場數(shù)據(jù)解釋中，場源邊界識別是不可或缺的內(nèi)容［1］。事實(shí)上，構(gòu)造單元劃分、斷裂構(gòu)造帶位置的確定、巖性識別與地層分布解釋、物性填圖等，均需確定邊界位置［2］。

邊界識別方法的發(fā)展可分為三個階段。第一階段的研究大多基于位場異常表面不同階次、不同方向上的導(dǎo)數(shù)及其組合，包括垂直導(dǎo)數(shù)［3］、水平總導(dǎo)數(shù)和二維解析信號［4］等。由于位場表面的導(dǎo)數(shù)隨場源深度的增加衰減較快，且處理實(shí)際數(shù)據(jù)時求導(dǎo)計(jì)算會放大噪聲，所以該階段方法的主要缺點(diǎn)是對異常幅值較小的場源識別能力較差且抗噪能力差。第二階段的方法主要針對邊界均衡識別，包括傾斜角法（Tilt angle）［5］、傾斜角的水平總導(dǎo)數(shù)法［6］、自動增益控制法［7］等。這類方法能同時識別淺部與深部的地質(zhì)體，但是提取的邊界比較發(fā)散，定位不夠準(zhǔn)確。在第三階段，學(xué)者們致力于在均衡識別的基礎(chǔ)上增強(qiáng)所提取的邊界，同時提高方法的抗噪能力，如Theta map法［8］、歸一化標(biāo)準(zhǔn)差（NSTD）法［9］、增強(qiáng)型均衡濾波器法［10］、水平總導(dǎo)數(shù)的傾斜角（TAHG）法［11］等。由此可見，邊界識別方法向更均衡、更收斂和更穩(wěn)定的方向發(fā)展。

各階譜矩及其統(tǒng)計(jì)不變量能夠詳細(xì)地描述表面幾何特征，早先被廣泛應(yīng)用于工程領(lǐng)域的表面形貌識別，目前已被引入地球物理學(xué)領(lǐng)域并取得了良好的應(yīng)用效果。孫艷云等［12］和楊文采等［13-14］基于二階譜矩定義了脊形化系數(shù)和邊界脊形化系數(shù)，從重力場成功地提取了地殼變形帶信息，其中邊界脊形化系數(shù)在構(gòu)造邊界定位中效果顯著。付麗華等［15］基于二階譜矩提出了一種描述表面數(shù)據(jù)粗糙度的特征因子，可有效反映數(shù)據(jù)波動與變異特征。付麗華等［16］基于四階譜矩，應(yīng)用算術(shù)平均頂點(diǎn)曲率提取磁異常數(shù)據(jù)特征，并將所提取的信息用于場源深度反演，取得了很好的效果。因此，利用譜矩能夠從位場數(shù)據(jù)中提取出與場源邊界、埋深等有關(guān)的幾何特征信息。

相較于其他方法，基于譜矩的位場幾何特征分析方法能夠更均衡、收斂地識別地質(zhì)體的幾何特征。二階譜矩可用于提取線形刻痕［12］，但不能很好地識別其他刻痕。應(yīng)用四階譜矩的幾何特征分析方法則可以進(jìn)一步完善對弧形等刻痕信息的識別。為此，提出四階弧形刻痕系數(shù)，理論和真實(shí)數(shù)據(jù)試驗(yàn)證明了利用此方法提取弧形和交會線形刻痕的有效性。最后結(jié)合二階譜矩的脊形化系數(shù)對塔里木盆地航磁數(shù)據(jù)進(jìn)行特征提取，更為全面地解釋了該區(qū)域的地殼構(gòu)造信息。

1 方法原理

1.1 表面譜矩

在笛卡爾直角坐標(biāo)系中，可用二維隨機(jī)過程函數(shù)z（x，y）表示位場表面，其中z為異常強(qiáng)度，（x，y）是平面直角坐標(biāo)位置。連續(xù)位場表面的p＋q階譜矩定義為［17-19］

式中：G（fx，fy）是位場表面z（x，y）的二維功率譜密度；fx、fy分別表示x、y方向上的空間頻率。

實(shí)際勘探中得到的位場數(shù)據(jù)是離散的、有限的。假定研究區(qū)域共有M×N個采樣點(diǎn)，在x、y方向上分別以Δx、Δy等間距采樣，此時位場表面函數(shù)可表示為z（xi，yj）（i＝1，2，…，M；j＝1，2，…，N），對這樣的離散位場表面，需要用滑動窗口提取局部特征信息。孫艷云等［12］推導(dǎo)出離散位場表面的p＋q階譜矩的表達(dá)式為

式中：M1、N1分別為滑動窗口在x、y方向包含的點(diǎn)數(shù)；F（fu，fv）為z（x，y）的二維離散傅里葉變換，F(xiàn)*（fu，fv）是其共軛。F（fu，fv）的傅里葉逆變換有如下性質(zhì)

式中IDFT表示離散傅里葉逆變換。因此窗口面元的二階譜矩元可用褶積形式表示為

式中：m20為表面x方向斜率的方差；m02為y方向斜率的方差；m11是x方向斜率與y方向斜率的協(xié)方差。類似地，可推導(dǎo)四階譜矩元m40、m31、m22、m13和m04的離散空間域表達(dá)式

四階譜矩元中，m40為x方向二階導(dǎo)數(shù)的方差，m04為y方向二階導(dǎo)數(shù)的方差，m31為與的協(xié)方差，m13為的協(xié)方差，m22為的方差或者的協(xié)方差。將計(jì)算結(jié)果定位于窗口中心點(diǎn)，逐點(diǎn)移動滑動窗口，即可得到每一個數(shù)據(jù)點(diǎn)的二階和四階譜矩元，從而提取到離散位場表面的特征信息。

1.2 基于譜矩的邊界識別方法

對于二階譜矩，位場表面可用統(tǒng)計(jì)不變量M2和Δ2表征［17，20-21］

式中M2表示位場表面斜率的方差，能夠反映位場異常的刻痕強(qiáng)弱，因此定義為刻痕的強(qiáng)度系數(shù)；Δ2主要與場在局部區(qū)域的各向異性有關(guān)，反映表面刻痕脊形化的程度，且通常Δ2≥0［12，22-23］。

為從位場表面準(zhǔn)確識別地殼變形帶，孫艷云等［12］提出了刻痕的脊形化系數(shù)

Λ2變化區(qū)間為［0，1］，與位場表面各向異性程度呈正相關(guān)，與刻痕強(qiáng)度系數(shù)M2成反比，可增強(qiáng)對異常幅值較小的場源的識別能力［12，22］。

類似地，四階譜矩的統(tǒng)計(jì)不變量M4和Δ4的計(jì)算公式［20］為

式中：M4表示表面曲率的方差，可反映表面曲率變化的劇烈程度，數(shù)值大說明表面明顯地尖銳，峰頂呈尖峰狀，數(shù)值小說明表面有較小的粗糙度變化，峰頂較圓滑；Δ4度量表面曲率的方向性效應(yīng)，對應(yīng)位場表面呈現(xiàn)圓弧形或者交會線形刻痕的幾何特征，值小表明位場表面這種特征不明顯。

與二階譜矩的脊形化系數(shù)不同，四階譜矩反映的是表面曲率變化，可提取圓弧形或者交會線形刻痕的幾何特征，這種特征不受位場幅度大小的影響。因此，可以用四階譜矩提取位場中的弧形刻痕信息，尤其是弧形的弱信號。為此，弧形刻痕系數(shù)必須與統(tǒng)計(jì)不變量Δ4正相關(guān)，與M4負(fù)相關(guān)。根據(jù)歸一化的要求，用四階譜矩的統(tǒng)計(jì)不變量M4和Δ4定義弧形刻痕系數(shù)

Λ4的變化區(qū)間亦為［0，1］。Λ4與成正比，表示與位場表面圓弧形刻痕呈正相關(guān)；與M4成反比，可增強(qiáng)對弱刻痕的識別能力。

2 理論模型試驗(yàn)

2.1 球體模型

為檢驗(yàn)譜矩方法提取弧形刻痕信息的有效性，設(shè)計(jì)圖1a所示模型。在大小為120km×120km、網(wǎng)格間距為1km×1km的區(qū)域內(nèi)有兩個半徑均為15km、中心埋深均為20km的球狀磁源體。球體A和B的球心坐標(biāo)分別為（－20km，－20km）和（20km，20km），磁化強(qiáng)度分別為0.8A/m和0.3A/m。該模型的理論磁異常見圖1b。

球體在地面的磁異常具有完整的弧形刻痕特征，是應(yīng)用四階譜矩提取位場中弧形刻痕信息的理想模型。采用3×3大小的滑動窗口對理論磁異常數(shù)據(jù)（圖1b）進(jìn)行特征提取，結(jié)果如圖1c～圖1f所示。圖1c為刻痕的強(qiáng)度系數(shù)M2，由式（13）可知水平導(dǎo)數(shù)越大其值越大，主要反映較強(qiáng)磁異常的邊界刻痕，對弱異常的識別能力較差，且提取的刻痕信息比較發(fā)散，不利于邊界的準(zhǔn)確定位。圖1d為脊形化系數(shù)Λ2，由圖可知，Λ2能均衡識別強(qiáng)弱磁異常，可以確定兩球體的中心位置，但是無法準(zhǔn)確定位其弧形邊界。圖1e為表面曲率的方差M4，主要反映較強(qiáng)磁異常的形狀與位置信息，對弱異常的識別能力較差。由圖1f可知，弧形刻痕系數(shù)Λ4能均衡識別強(qiáng)弱磁異常，其高值清晰收斂地刻畫出磁源體的邊界。圖中球體B因磁化強(qiáng)度較小而產(chǎn)生的異常較弱，但在Λ4的提取結(jié)果中，它的弧形刻痕信息和球體A的一樣清晰，證明了應(yīng)用四階譜矩能夠提取位場中的弧形刻痕的弱信號。

實(shí)際情況中地質(zhì)體的埋深是未知的，因此邊界識別方法對深度越不敏感越好，說明地質(zhì)體埋深對提取效果的影響不大。為進(jìn)一步檢驗(yàn)弧形刻痕系數(shù)Λ4識別弧形邊界的有效性，改變圖1a所示模型中球體的埋深，測試Λ4對深度的敏感程度。圖2a～圖2d是球體中心埋深分別為20、22、24、26km時的理論磁異常和Λ4提取結(jié)果。由于計(jì)算結(jié)果的范圍變化較大，圖中根據(jù)場值的大小采用不同的色標(biāo)。

圖1 球體模型及譜矩方法不同異常參數(shù)提取結(jié)果

圖2 不同深度球體模型的理論磁異常（左）和弧形刻痕系數(shù)Λ4的提取結(jié)果（右）

由圖2可知，隨著深度增加，原本獨(dú)立的兩磁異常范圍漸漸擴(kuò)大并產(chǎn)生交疊，場源強(qiáng)度逐漸變小，增加了邊界識別的難度。分析Λ4的提取結(jié)果可知，弧形分布信息受場源強(qiáng)度影響較小，能穩(wěn)定有效地反映不同深度地質(zhì)體的弧形邊界信息，對深度的敏感度較低。

為對比弧形刻痕系數(shù)Λ4與常用邊界識別方法的效果，本文選取水平總導(dǎo)數(shù)（TDX）［9］、水平總導(dǎo)數(shù)的傾斜角TAHG［11］、Theta map［8］和歸一化標(biāo)準(zhǔn)差（NSTD）方法［9］對模型（圖1a）進(jìn)行邊界信息提取，結(jié)果如圖3a～圖3d所示。

由圖3可知，四種傳統(tǒng)方法都能提取到兩個球狀磁源體的邊界，但是均存在邊界發(fā)散的問題，除TDX外的另外三種方法的提取結(jié)果中還存在干擾刻痕，難以準(zhǔn)確分辨地質(zhì)體的邊界位置。對比圖1f可知，Λ4相較于傳統(tǒng)方法能更準(zhǔn)確、收斂地刻畫球狀磁源體的邊界，且沒有干擾刻痕。試驗(yàn)結(jié)果表明，Λ4的邊界識別效果明顯優(yōu)于傳統(tǒng)邊界識別方法。

2.2 長方體模型

為檢驗(yàn)譜矩方法提取線形刻痕信息的有效性，設(shè)計(jì)圖4a所示包含兩個長方體狀磁源體（C和D）的組合模型進(jìn)行理論試驗(yàn)。研究范圍為300 km×300km、網(wǎng)格間距為1km×1km。磁源體C的中心坐標(biāo)為（0，5.32km，9.00km），頂面埋深為4km，長、寬、高分別為150.00、20.22、10.00km，磁化強(qiáng)度為0.5A/m；磁源體D的中心坐標(biāo)為（0，－79.79km，13.00km），頂面埋深為8km，長、寬、高分別為179.78、20.21、10.00km，磁化強(qiáng)度為0.2A/m。圖4b為該模型理論磁異常。

據(jù)圖4可見，應(yīng)用Λ2能夠提取位場中的弱板形刻痕信號，應(yīng)用Λ4也能夠提取位場中的弱板形刻痕信號，不過邊界不如Λ2清晰，這是因?yàn)棣?主要反映的是磁異常表面斜率的變化，而Λ4主要反映的是曲率的變化，在長方體模型的平直邊界處，磁異常斜率變化劇烈，曲率變化相對較小。僅僅利用Λ4，則可清晰地反映板形體的端角位置，而不是邊界。板形體的端角呈現(xiàn)交會線形刻痕的幾何特征，能夠作為位場中的弧形刻痕信號被提取出來。

圖3 常用邊界識別方法對球狀地殼模型的提取結(jié)果（圖中黑色虛線是球體邊界）

圖4 長方體模型及譜矩方法不同異常參數(shù)提取結(jié)果

3 實(shí)際應(yīng)用

楊文采等［22］曾經(jīng)介紹過塔里木盆地重力場的高階譜矩分析，現(xiàn)在利用四階譜矩的弧形刻痕系數(shù)分析該地區(qū)的航磁場。航空物探數(shù)據(jù)已經(jīng)過標(biāo)定和同化處理［24］，經(jīng)度和緯度方向點(diǎn)距分別為0.00886°和0.03380°，采樣點(diǎn)數(shù)為678×356。磁力化極異常如圖5a所示。采用3×3的滑動窗口對其進(jìn)行特征提取，Λ2和Λ4的提取結(jié)果分別見圖5b和圖5c。對比Λ2與Λ4的提取結(jié)果可見，二者有明顯的區(qū)別，Λ2主要反映盆地基底中的線形構(gòu)造，而Λ4主要反映出盆地內(nèi)部主體為等軸狀構(gòu)造，如二疊紀(jì)古火山群等。

圖5 塔里木盆地磁異常（a）、Λ2（b）和Λ4（c）分布圖

由于異常疊加及反演固有的多解性，準(zhǔn)確解釋區(qū)域磁異常非常困難，計(jì)算磁異常源的埋藏深度可以對地質(zhì)成果的解釋提供更準(zhǔn)確的依據(jù)［25-27］。楊文采等［24］應(yīng)用三維歐拉反褶積計(jì)算了區(qū)域磁異常源的埋藏深度。由于克拉通盆地地層通常具有上新下老的規(guī)律，將磁異常源分解為三個深度層次，分別為淺層2～5km、中層5～10km和深層10～20km。淺層磁異常源分布在盆地邊緣，由中生代構(gòu)造運(yùn)動和海西期玄武巖侵位共同作用而成，在圖5a和圖5b中盆地邊緣的線性構(gòu)造中有明顯反映。中層磁異常源分布在盆地內(nèi)部，主要是海西期玄武巖侵位形成的，盆地內(nèi)部5～10km的海相地層中常有發(fā)現(xiàn)。深層磁異常源主要是古老基底變質(zhì)作用形成的，反映為圖5a和圖5b中盆地南部的線性構(gòu)造。在磁異常圖（圖5a）中，盆地內(nèi)部深度為5～10km的海相地層中玄武巖侵位圍繞古火山群分布，表現(xiàn)為等軸狀構(gòu)造的弱異常，淹沒在古老變質(zhì)基底的深層磁異常中，難以定位。在圖5c中Λ4強(qiáng)異常表現(xiàn)為紅色圓圈狀刻痕和圈點(diǎn)（如圖中箭頭所示），可能反映古火山群和其他等軸形狀的地質(zhì)體，說明應(yīng)用四階譜矩分析和信息提取技術(shù)成功地揭示了新的地殼組構(gòu)。

4 結(jié)論

（1）基于譜矩的位場幾何特征分析方法是近些年提出的一種新的自動邊界識別方法，相較于傳統(tǒng)方法能夠更均衡、收斂地識別地質(zhì)體的幾何特征。二階譜矩可用于提取線形刻痕，但不能很好地識別其他刻痕。應(yīng)用四階譜矩的幾何特征分析方法可以進(jìn)一步完善弧形等刻痕信息的識別。

（2）與刻痕的脊形化系數(shù)不同，弧形刻痕系數(shù)反映的是圓弧形或者交會線形刻痕的幾何特征，這種特征獨(dú)立于位場的幅度。因此，可以應(yīng)用四階譜矩提取位場中的弧形刻痕信息，尤其是弧形的弱信號。

（3）球體在地面的磁異常具有完整的弧形刻痕特征，是應(yīng)用四階譜矩來提取位場中弧形刻痕信息的理想模型。理論模型試驗(yàn)證明了應(yīng)用四階譜矩能夠提取位場中的弧形刻痕的弱信號。

（4）在塔里木盆地磁異常圖中，盆地內(nèi)部深度為5～10km的海相地層中玄武巖侵位圍繞古火山群分布，表現(xiàn)為等軸狀構(gòu)造的弱異常，淹沒在古老變質(zhì)基底的深層磁異常中，難以定位?；⌒慰毯巯禂?shù)強(qiáng)異常可反映古火山群和其他等軸形狀的地質(zhì)體，證明應(yīng)用四階譜矩分析和信息提取技術(shù)可以揭示特殊的地殼組構(gòu)。