付麗華 曾 誠 楊文采 李燕梅

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

0 引言

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

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

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

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

1 方法原理

1.1 表面譜矩

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

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

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

式中：M1、N1分別為滑動窗口在x、y方向包含的點數(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é)方差。類似地，可推導四階譜矩元m40、m31、m22、m13和m04的離散空間域表達式

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

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

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

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

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

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

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

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

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

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

2 理論模型試驗

2.1 球體模型

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

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

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

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

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

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

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

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

2.2 長方體模型

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

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

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

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

3 實際應用

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

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

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

4 結論

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

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

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

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