胡 波，朱谷昌，張遠飛，肖 燃

空間U統(tǒng)計量法在遙感蝕變信息提取中的應(yīng)用研究

胡 波1，2，朱谷昌1，2，張遠飛2，肖 燃3

(1.中南大學(xué)，長沙 410083;2.有色金屬礦產(chǎn)地質(zhì)調(diào)查中心，北京 100012;3.江西應(yīng)用工程職業(yè)學(xué)院，萍鄉(xiāng) 337000)

鑒于傳統(tǒng)蝕變異常信息提取方法的局限性，試圖在二維光譜特征空間(即散點圖)中研究遙感蝕變信息的提取方法。在二維散點圖中，兩個波段灰度值的聯(lián)合分布體現(xiàn)了各向異性的特點，聚類常呈類似橢圓形分布。對散點圖中點的頻數(shù)使用空間U統(tǒng)計量法逐個獲得地物聚類的橢圓參數(shù)，將每個橢圓內(nèi)的點映射到遙感圖像上，經(jīng)過目視解譯，最終獲得蝕變異常信息的空間分布。以青海省巴音山地區(qū)蝕變信息提取為例，對空間U統(tǒng)計量方法進行了應(yīng)用，并將提取的異常區(qū)域與已獲得的各種數(shù)據(jù)進行對比，發(fā)現(xiàn)鐵化和泥化區(qū)域吻合度較高，取得了比主成分分析更好的蝕變信息提取效果。

蝕變異常信息提取;U統(tǒng)計量法;聚類各向異性;二維散點圖;異常聚類

0 引言

當(dāng)前，用于遙感蝕變異常信息提取的方法主要有基于主成分分析法和波段比值法兩種。由于這兩種方法采用的是相對固定的處理流程和分析手段，往往不能夠因地制宜地提取遙感蝕變異常信息。張玉君［1］等人對幾種傳統(tǒng)蝕變信息提取方法做了總結(jié)，認(rèn)為應(yīng)用比值法得到的結(jié)果不可避免地會出現(xiàn)難以解釋的假異常，故提出采用以主成分分析為主，輔以比值法和光譜角法的“去干擾異常主分量門限化技術(shù)流程”;張遠飛［2］等人的研究發(fā)現(xiàn)，基于掩模+主成分分析+比值法的方法適用環(huán)境存在局限性。從光譜特征空間來說，傳統(tǒng)的方法對于偏離出背景的蝕變信息捕捉效果較好。

筆者嘗試在二維光譜特征空間中研究遙感蝕變異常信息的提取。二維散點圖是二維光譜特征空間的一種平面表現(xiàn)形式，它所反映的是2個波段灰度值的聯(lián)合分布。同類地物無論在空間上分布如何，其在散點圖當(dāng)中總是形成聚類，并且呈現(xiàn)出近似橢圓的形態(tài);而不同地物組合在一起則形成橢圓的套疊結(jié)構(gòu)。找到蝕變異常信息所構(gòu)成的聚類橢圓，然后將其映射到遙感圖像上，就可以獲得蝕變異常信息的空間分布。

空間U統(tǒng)計量法是一種用于多重分形的局部奇異性分析方法，它能夠獲取局部位置的各向異性參數(shù)［3，4］，并且可以有效減少地物聚類橢圓之間交疊區(qū)域被誤判的概率。在散點圖中使用該方法可以獲得各個聚類橢圓的參數(shù)，依次找到并剔除背景和各種干擾，最終提取到蝕變異常信息。

本文以青海省都蘭縣巴音山地區(qū)為例，采用空間U統(tǒng)計量法提取了鐵化和泥化蝕變異常信息，并將其與該地區(qū)的已知數(shù)據(jù)進行對比，驗證方法的應(yīng)用效果。

1 光譜特征空間遙感數(shù)據(jù)的結(jié)構(gòu)特點

1.1 光譜特征空間的相關(guān)定義

光譜特征空間是指在一幅多光譜圖像中假設(shè)有n個波段，則每一個像元在各個波段上的灰度值將構(gòu)成一個向量，用 X=(x1，x2，…，xn)T表示。包含所有X的n維空間稱為光譜特征空間［5］。本文重點研究的是二維光譜特征空間，如TM1與TM3的組合，那么每個像元矢量即為X=(xTM1，yTM3)T。以該像元在TM1波段上的灰度值作為x坐標(biāo)，在TM3波段上的灰度值作為y坐標(biāo)，包含所有X的空間就是本文所說的二維光譜特征空間(即二維散點圖，文中簡稱散點圖)。散點圖中(x，y)處的值代表了遙感圖像中TM1=x且TM3=y這種聯(lián)合分布的頻數(shù)。把頻數(shù)作為z軸，可以在三維空間上觀察地物聚類的形態(tài)和分布特征，文中稱立體散點圖。

1.2 背景、干擾與蝕變異常的定義

文獻［6］基于光譜特征空間給出了背景、干擾與蝕變異常的定義。背景是指在多波段遙感圖像數(shù)據(jù)點集空間中的超橢球體(實際上為近似超平面體)，而游離在超橢球體之外的“小類”點集則為干擾或蝕變異常目標(biāo)。從統(tǒng)計角度來看，遙感圖像是由散點圖中高頻數(shù)的背景和相對低頻數(shù)的異常組成。背景是圖像上頻數(shù)最高的區(qū)域，在立體散點圖當(dāng)中呈“主峰”形態(tài);異常是相對于背景定義的，在“主峰”周邊有時會出現(xiàn)一些相對獨立的“次峰”形態(tài)，這就是異常。本次研究中蝕變異常信息提取的目標(biāo)就是首先在光譜特征空間中定位蝕變異常，然后使用一定的技術(shù)方法獲得其空間分布。

1.3 光譜特征空間遙感數(shù)據(jù)的幾何特征

在二維散點圖中，地物灰度的聯(lián)合分布會在不同方向上產(chǎn)生變形，形成的聚類體現(xiàn)出各向異性的特點，往往近似一個橢圓。依照多元高斯分布等值線描述［7］:多元高斯密度的等值線是中心在均值向量m處的超橢球平面。該超橢球面的主軸平行于特征向量e1，而特征向量a是相應(yīng)的方差，如圖1所示。

圖1 多元高斯密度等值線Fig.1 Isoline of multi－ Gaussian density

假設(shè)地物灰度的概率密度服從一定的正態(tài)分布，在2個波段組合中，這種地物的灰度概率密度服從二元高斯分布。散點圖上，每類地物會形成橢圓形的聚類。

眾多地區(qū)的散點圖上都出現(xiàn)了橢圓的套疊結(jié)構(gòu)，頻數(shù)最高區(qū)域?qū)?yīng)的是背景，而周邊的頻數(shù)相對較低，變化較緩的橢圓對應(yīng)的是異常。假設(shè)每一類地物的灰度服從正態(tài)分布，那么每個波段的灰度直方圖是疊加之后的多種地物灰度的分布?？梢宰鞒鲆韵峦茢?

(1)在散點圖當(dāng)中，各類地物分別服從二元高斯分布，但由于受到地物混合影響，其聯(lián)合分布形態(tài)發(fā)生變化，但仍近似于橢圓形;

(2)遙感圖像是多種地物的整體映射，故在散點圖中應(yīng)表現(xiàn)為橢圓的套疊。

基于二維光譜特征空間中遙感信息呈現(xiàn)的近似橢圓形的空間分布特征，使用空間U統(tǒng)計量法獲得各向異性參數(shù)，并且定位蝕變異常的聚類橢圓，最終便可提取到蝕變異常信息。

2 空間U統(tǒng)計量法

U統(tǒng)計量的構(gòu)造為

式中，z0表示在分類范圍內(nèi)異常總體和背景總體的分類閾值，通常采用局部的均值。一般情況下，當(dāng) U［B(x，y)(r，β，θ)］＞ 0 時，認(rèn)為樣本屬于異?？傮w;反之可以判定樣本屬于背景總體。|U|越大，則表明分類效果越好;而|U|取值接近0，則意味著樣本內(nèi)存有混合現(xiàn)象。

最佳的窗口參數(shù)應(yīng)當(dāng)是U值取得最大時所對應(yīng)的 U*［B(x，y)(r，β，θ)］，即

此時的橢圓參數(shù)(r，β，θ)就是所求。隨著窗口半徑的加大，U*［B(x，y)(r，β，θ)］的值應(yīng)該不斷增大才有意義。如果取值減小后又繼續(xù)增大，就意味著窗口已經(jīng)跨越了邊界，應(yīng)該在值出現(xiàn)減小時就終止運算，所以，有必要用生成的U*－r曲線圖作為輔助(圖2)。

圖2 U*－r曲線示例Fig.2 Instance of U*－ r curve

圖2上曲線1表示異?？傮w，在|U*|值減小前也就是rA處取得極值;曲線2表示背景總體，在|U*|值減小前也就是rB處取得極值，當(dāng)|U*|接近0時的r'，表明背景和異常有混雜現(xiàn)象。

根據(jù)空間U統(tǒng)計量的構(gòu)造方法可以看出，對于一個樣本，它的期望μ沒有改變，而方差σ2［U］=σ2/n。故，當(dāng)判斷一個樣本是屬于背景總體還是異?？傮w時，顯著縮小了出現(xiàn)判定錯誤的區(qū)間(圖3)。

圖3 用樣本值和U作為統(tǒng)計量的分布函數(shù)Fig.3 Distribution function by statistics of sample value and U value

圖3 中α表示A類樣本被判定為B類的概率;β表示B類樣本被判定為A類的概率。不難看出，采用U統(tǒng)計量之后，出現(xiàn)兩類錯誤的概率顯著減小了。

3 應(yīng)用實例

以青海省都蘭縣巴音山銅、鉛、鋅、銀多金屬礦集區(qū)為實驗區(qū)，該區(qū)位于柴北緣火山弧裂陷構(gòu)造帶、東昆侖北構(gòu)造帶及鄂拉山構(gòu)造帶三者交匯處，那里巖漿活動強烈且?guī)r類復(fù)雜，成礦物質(zhì)來源廣泛，成礦作用類型多種，礦產(chǎn)種類較多。

該區(qū)出露的地層主要有古元古界沙柳河群(Pt1dk)片巖、麻巖，古生界上奧陶統(tǒng)灘澗山群(O3tn)以基性、中基性、酸性火山巖及其碎屑巖類為特點的細碧角斑巖類，間夾長英質(zhì)碎屑巖、硅質(zhì)巖、大理巖透鏡體的巖性組合，變質(zhì)后一般為綠片巖和絹云母石英片巖、千枚巖及大理巖等。其上不整合覆蓋有泥盆系上統(tǒng)砂巖、礫巖和千枚巖。

本區(qū)構(gòu)造以柴北緣火山弧裂陷—裂谷構(gòu)造帶及其后期繼承性近東西向構(gòu)造線為主導(dǎo)格局。區(qū)內(nèi)構(gòu)造較為復(fù)雜，為礦床定位造就了良好條件。除下古生界廣泛發(fā)育雙峰式火山巖系外，加里東晚期和海西晚期均有花崗巖體沿東西向斷裂或?qū)娱g斷裂帶侵入，為成礦物質(zhì)流體的活動提供了熱動力條件。

綜合分析，阿爾茨托山東段的巴音山一帶是重要的多金屬成礦帶。該地區(qū)有著詳實的野外考察資料和豐富的工作總結(jié)，擁有充分的地質(zhì)、物探、化探及遙感等數(shù)據(jù)儲備，因而選擇該地區(qū)作為已知地區(qū)進行新蝕變信息提取方法的實驗和評價。

傳統(tǒng)方法使用 TM1、TM3、TM4、TM5等波段組合進行主成分提取，并通過判定確定PC4是代表鐵化蝕變的主分量，利用最優(yōu)彩色密度分割的方法，最終取得鐵化蝕變的區(qū)域。通過與各種數(shù)據(jù)比較，發(fā)現(xiàn)提取到的鐵化蝕變區(qū)域與巴音山頂、東山及烏龍灘南部等幾個大型礦區(qū)吻合良好，但是在沙柳河?xùn)|岸并處于山脊線西側(cè)的幾個小型礦點，與之對應(yīng)的蝕變卻全部出現(xiàn)在山脊線東側(cè)，偏離較大。同樣，使用TM1、TM4、TM5、TM7等波段提取泥化蝕變信息，其聚類橢圓特點與鐵化信息的類似。

［8］，筆者對于巴音山地區(qū)采用TM4與TM3波段組合提取鐵化異常信息，用TM7與TM5波段組合提取泥化異常信息。

從TM4和TM3的立體散點圖(圖4)不難看出，此波段組合下的聯(lián)合灰度分布有明顯的聚類特征。圖像中的“主峰”對應(yīng)背景信息，其周邊相對低矮的“次峰”對應(yīng)異常信息。

圖4 TM3與TM4波段組合立體散點圖Fig.4 3D scatter plot by TM3 and TM4

圖5 是在TM4和TM3散點圖中提取各類地物的流程，其中藍—黃—綠—紅—品紅的過渡表示頻數(shù)由低至高。

通過使用空間U統(tǒng)計量法，首先可以得到聚類橢圓I。這是圖像上頻數(shù)最高的部分(即背景)，為避免其對之后提取聚類的影響，故先將其剔除。通過剔除局部最高的峰值，之前被高峰所掩蓋的較低峰的區(qū)域變成最高峰，變?yōu)榧t色顯示出來。當(dāng)剔除背景橢圓I之后，圖像中另外兩個橢圓逐漸顯露出來。如此通過不斷地“捕捉橢圓—去除橢圓—捕捉橢圓”，層層剝離背景和干擾，最終使得遙感蝕變異常這類弱信息“水落石出”。

圖5 巴音山TM3與TM4組合散點圖中的聚類提取流程Fig.5 Extraction process of clusters in scatter plot by TM3 and TM4 in Bayinshan

經(jīng)過上述提取流程，最終確定了6個聚類橢圓。圖6是散點圖上各個橢圓的位置。表1列出了各個橢圓的參數(shù)。將各橢圓映射到圖像上進行目視解譯，對聚類橢圓 I、III、IV、VI分別定性判斷為背景、淺層植被、裸露河床、深厚植被。橢圓II所對應(yīng)的主要是沿沙柳河、烏龍灘東側(cè)地區(qū)以及巴音山東側(cè)河流沿岸地區(qū)，該地區(qū)徑流密布(基本干涸);橢圓V對應(yīng)沿巴音山的山脊線兩側(cè)分布的區(qū)域、東山頂部以及哈莉哈德山頂部區(qū)域。

圖6 TM3與TM4提取的6個橢圓Fig.6 6 ovals extracted from TM3 and TM4

表1 TM3與TM4散點圖中的聚類橢圓參數(shù)Tab.1 Cluster ovals’parameters of scatter plot by TM3 and TM4

圖7中黃色曲線框中的異常區(qū)域?qū)?yīng)橢圓V，其余對應(yīng)橢圓II。這兩個橢圓所對應(yīng)地區(qū)共同特點是礦產(chǎn)豐富且有多處已探明的礦床。從已知礦點分布來看，橢圓II和橢圓V同屬蝕變異常的信息聚類，但卻屬于不同的聚類橢圓。據(jù)考察資料顯示，這些地區(qū)蝕變均以硅化、絹云母化、綠泥石化及褐鐵礦化等為主，蝕變礦物一致，所以蝕變異常出現(xiàn)兩個聚類不是由于蝕變礦物不同造成的，而是受河流泥沙類物質(zhì)影響。研究區(qū)域河流基本干涸，河床裸露，河道中的泥沙類物質(zhì)反射率較高，使得裸露河床出現(xiàn)在散點圖右上方(橢圓IV)。野外記錄顯示，河道周邊確實存在一定量的蝕變物質(zhì)并且與泥沙類物質(zhì)混合，這就造成沿沙柳河及烏龍灘局部鐵化蝕變異常區(qū)域反射率較高。而沒有受到泥沙類物質(zhì)影響的巴音山頂部、東山頂部和哈莉哈德山頂部反射率相對較低。從整體形態(tài)來看，橢圓II和橢圓V同在散點圖對角線上側(cè)，符合Fe3+在TM3波段反射率高，在TM4波段反射率低的物理特征。

另外，從圖6中還可以看出，有兩個植被的聚類，淺層植被對應(yīng)植被覆蓋稀疏地區(qū)，深層植被代表植被茂盛地區(qū)。淺層植被地區(qū)在散點圖中的位置連接著深層植被與背景，充分反映了其植物、巖石及土壤混雜的特點。而對于深層植被地區(qū)，因為植被十分茂盛幾乎遮蔽了其余地物的反射光，故其光譜特征較為獨立，距離背景橢圓I較遠。

將橢圓II、橢圓V映射到遙感圖像上便獲得了鐵化蝕變異常區(qū)域。對比該地區(qū)物探和化探數(shù)據(jù)，發(fā)現(xiàn)提取結(jié)果與之高度對應(yīng)，并且與野外記錄相符。圖7是提取到的鐵化蝕變異常與已知礦點分布疊加圖。圖上藍色—黃色—綠色—紅色—品紅過渡色表示蝕變異常由弱至強。經(jīng)統(tǒng)計，在21個已知礦點當(dāng)中，有19個毗鄰或者存在于鐵化異常區(qū)，命中率達90.4%。距離鐵化蝕變區(qū)域較遠的有2個，分別是吉給申溝鉛礦床(野外勘察顯示該處礦體規(guī)模較小，為共生礦產(chǎn))和藏碑溝鉛、鋅、銀多金屬礦點(未見詳細記錄)。

圖7 鐵化蝕變與已知礦點分布疊加圖Fig.7 Stacking chart of ironed alteration and known deposits

利用TM7與TM5影像組合的散點圖提取泥化蝕變異常過程與上述提取鐵化信息的類似。泥化信息也出現(xiàn)了兩個聚類橢圓，與鐵化蝕變信息的聚類橢圓呈同樣的特點，此處不再贅述。

將蝕變異常信息映射在遙感圖像上，經(jīng)統(tǒng)計，在21個礦點中，有20個在泥化異常區(qū)內(nèi)，命中率達95.2%，僅吉給申溝鉛礦床未出現(xiàn)在蝕變區(qū)域內(nèi)。

疊加兩類異常，發(fā)現(xiàn)在礦點周圍基本都有異常分布，從而肯定了提取結(jié)果的正確性。在沙柳河?xùn)|岸、巴音山北部、哈莉哈德山及東山南部等地區(qū)(圖7中紅色框范圍)，泥化和鐵化均有廣泛分布，但是礦體分布有限，需要進一步工作查明原因。而烏龍灘北部等其余藍框范圍內(nèi)雖然有大片蝕變分布，但這是由于河流淤積造成的蝕變礦物聚集，考察意義不大。

對比主成分分析方法的提取效果，本文方法提取結(jié)果更加準(zhǔn)確，尤其是在沙柳河?xùn)|岸的幾個礦點，異常與礦點吻合度較高。

4 結(jié)論

(1)從光譜特征空間分析入手，地物分布特征表現(xiàn)的更加直觀。在散點圖中，原本空間上分散的地物能夠形成聚類，并且在立體散點圖中形成局部的“高峰”。通過觀察二維散點圖和立體散點圖的圖形特征，蝕變異常信息提取可以做到有的放矢。

(2)蝕變異常信息在散點圖上可能出現(xiàn)不止一個聚類，這取決于蝕變礦物的光譜以及諸如礦物粒度、地物混合情況、含水量及陰影等綜合因素的影響。有必要加強對研究地區(qū)地理現(xiàn)象的整理和分析，以避免蝕變異常提取過程中錯提、漏提現(xiàn)象。

(3)應(yīng)用空間U統(tǒng)計量法提取蝕變異常信息是一種新嘗試，還要進一步加強光譜混合機制和遙感數(shù)據(jù)空間結(jié)構(gòu)的研究，加深對地物聚類規(guī)律的把握，以提升蝕變異常信息提取的可信度。

參考文獻:

［1］ 張玉君，曾朝銘，陳 薇.ETM+(TM)蝕變遙感異常提取方法研究與應(yīng)用——方法選擇和技術(shù)流程［J］.國土資源遙感，2003(2):44－50.

［2］ 張遠飛，吳健生.基于遙感圖像提取礦化蝕變信息［J］.有色金屬礦產(chǎn)與勘查，1999，8(6):604 －606.

［3］ Cheng Q M，Agterberg F P，Bonham － Carter G F.A Spatial Analysis Method for Geochemical Anomaly Separation［J］.Journal of Geochemical Exploration，1996，56(3):183 －195.

［4］ 陳志軍.多重分形局部奇異性分析方法及其在礦產(chǎn)資源信息提取中的應(yīng)用［D］.武漢:中國地質(zhì)大學(xué)，2007.

［5］ 中國知網(wǎng).CNKI概念知識元庫［DB/OL］.http//define.cnki.net/，2006.

［6］ 張遠飛，吳德文，朱谷昌，等.遙感蝕變信息檢測中背景與干擾問題的研究［J］.國土資源遙感，2008(2):22－26.

［7］ Therrien C W.Discrete Random Signals and Statistical Signal Processing［M］//周宗潭，董國華，徐馨，等.獨立成分分析.北京:電子工業(yè)出版社，2007:24.

［8］ 張玉君，楊建民，陳 薇.ETM+(TM)蝕變遙感異常提取方法研究與應(yīng)用——地質(zhì)依據(jù)和波譜前提［J］.國土資源遙感，2002(4):30－36.

The Application of Spatial U－static Method to the Extraction of Alteration Anomalies

HU Bo1，2，ZHU Gu － chang1，2，ZHANG Yuan － fei2，XIAO Ran3
(1.Central South University，Changsha 410083，China;2.China Non － ferrous Metals Resource Geological Survey，Beijing 100012，China;3.Jiangxi Application Engineering Vocational College，Pingxiang 337000，China)

The authors tried to extract alteration anomalies in spectral characteristic space(scatter plot)in view of the limitation of the traditional methods.The scatter plot takes on an anisotropic feature in associated distribution of RS data’s grey scale.The distribution is usually combined by oval clusters.Parameters of oval clusters are acquired sequentially by applying U －Statistic method in the frequency of the scatter plot.Through mapping the points inside the oval into RS image and interpreting visually，the spatial distribution of alteration anomalies is obtained eventually.In this paper，this new method was described with the instance of Bayinshan area in Qinghai province.Other data acquired were also comparatively studied，and it is found that the anomalies of ferruginization and argillation are consistent well with each other.This new method has a better performance than PCA in the study area.

Alteration information extraction;U－static method;Anisotropy of cluster;2D scatter plot;Anomaly cluster

TP 753

A

1001－070X(2011)03－0071－06

2010－11－30;

2011－02－05

胡 波(1986－)，男，碩士研究生，地圖學(xué)與地理信息系統(tǒng)專業(yè)，主要從事GIS和RS的應(yīng)用研究。

(責(zé)任編輯:刁淑娟)