高光譜遙感數據的處理及其在鈾資源勘查的應用
—以廣西苗兒山地區為例
劉德長 謝紅接 趙英俊 張進業
0引言
高光譜遙感數據是80年代以來發展起來的新一代(第三代)遙感數據。它使遙感技術步入可以同時獲取地球表面物質成分信息和空間分布特征信息的新階段。本文以廣西苗兒山地區機載成像細分紅外高光譜數據為例,探討了高光譜遙感圖像的預處理和深化處理,野外實測光譜曲線庫建立與典型地物的波譜特征,以及高光譜數據在鈾資源勘查中的應用等。
1 高光譜遙感數據的預處理和深化處理
1.1 苗兒山地區高光譜遙感數據特點
上海技術物理所航空紅外細分光譜掃描儀是以地質遙感為主要應用目的開發的,在2~2.5um光譜區間細分為12個工作波段。廣西苗兒山地區高光譜遙感數據,由中南地質局230所和上海技術物理所于1988年11月2日、4日獲取。平臺為伊爾14型飛機,航速為250km/h,航高(絕對航高)達4500m,沿NE30。方向飛行。所選8個通道的中心波長分別為1.600μm、2.143μm、2.200μm、2.250μm、2.300μm、2.330μm、2.370μm和2.450μm。通過分析,苗兒山地區高光譜遙感數據有如下特點:
(1)由于儀器本身的原因,所有數據沿航向具周期性條帶噪聲;
(2)與TM和MSS等圖像不同,高光譜圖像在同一航帶的各波段之間,沿航向和掃描方向有幾個至十幾個象元之位移;
(3)由于航向分辯率(8.68m)和掃描方向分辯率(27m)的不一致,象元呈長方形,其長(航向)寬(掃描向)比為27:8.68=3.11:1;
(4)由于飛行方向、掃描角度、太陽方位角和太陽高度的幾何關系的變化和大氣的影響,掃描圖像存在邊緣輻射差異;
(5)由于飛行速度不等和掃描角度的變化引起圖像幾何差異;
(6)第7航帶數據由于記錄方式的錯亂,引起圖像與正確方式所成的圖像呈鏡向對稱關系;
(7)高光譜圖像的一個顯著特點是能夠從圖像中提取任意象元的光譜率曲線,但由于成像光譜儀獲取的圖像象元亮度值并不能代表地物的反射率值,需要進行反射率轉換。
1.2 高光譜數據的預處理
在研究上述所有特點的基礎上,開發了一系列數據預處理軟件,并進行了條帶去除、波段間配準、航向壓縮、鏡向變換、輻射校正、相對反射率轉換等預處理,最終得到清晰的圖像。
1.2.1去條帶方法研究
(1)原始數據分析 原始數據存在嚴重的條帶干擾,垂直于航線方向分布,且從第一通道到第八通道條帶越來越增強。經分析,這種條帶是由于在成像過程中機械或光學部件不穩定造成的,具有明暗交替的特征,并具有一定的周期性,以9個象素為周期,呈正弦波狀,部分以8個象素為周期。
(2)干擾條帶的消除方法 可以將干擾條帶作為加性噪聲處理以消除條帶,要求盡量減少信息損失,一般有兩種處理方法:
①平滑法,可以采用9點平滑或9×9點,但9×9窗口平滑效果不好,在航線方向上以9點平滑效果好,但也會損失有效信息。
②濾波法,從前面的分析可知干擾條帶的周期為9,通過不斷自動調整波形振幅,可以有效地消除干擾,研究過程中使用了濾波法。
(3)效果分析圖 圖1左圖為原始圖像,右圖為去條帶后的圖像,對比左右兩側的圖像可以明顯地看出干擾條帶被有效地消除了,而原始數據信息的損失比較小。條帶消除后圖像的列數保持不變,但行數減少,如第5航帶行數由5712減為5632。

圖1去條帶影像效果對比圖
(左側為原圖像,右側為處理后圖像)
1.2.2 波段之間配準
一般情況下,TM圖像的波段與波段之間都是完全重合的,不存在象元位移,而苗兒山高光譜遙感圖像存在著明顯的象元位移,一般為幾個→十幾個象元。從5、6、7、9條航帶來看,都是第2波段(2.143μm)和第3波段(2.2μm)之間象元位移量最小,如第5航帶的第2,3波段則完全重合,這可能說明該儀器的第2,3通道工作狀況最好。下面以第9航帶為例來說明波段之間的位移大小及配準方法。以第2波段為基準,其他波段則向左右或上下位移。
以第2波段為基準
第1波段 右移1個 下移 7個
第3波段 左移2個 上移4個
第4波段 左移4個 下移14個
第5波段 左移4個 下移7個
第6波段 左移5個 下移8個
第7波段 右移2個 下移3個
根據以上位移量截取完全重合的圖像部分如表1,得到第9航帶最高配準后圖像大小為列240,行5738。
表1 第9航帶各波段截取部分
|
列數 |
行數 |
第1波段 |
9+2:256-6 |
8:6125-11 |
第2波段 |
9+3:256-5 |
15:6125-4 |
第3波段 |
9+5:256-3 |
19:6125-0 |
第4波段 |
9+4:256-4 |
1:6125-18 |
第5波段 |
9+7:256-1 |
8:6125-11 |
第6波段 |
9+8:256-17 |
7:6125-12 |
第7波段 |
9+1:256-17 |
13:6125-7 |
(說明列數的9是由于每幅圖像都有9象元的亮邊,所以都被裁掉)
1.2.3 航向壓縮和掃描方向拉伸處理
機載成像光譜儀沿飛行方向的分辯率由飛機的飛行速度決定,其速度為250公以/小時(即250000m/3600s),掃描率為每秒8線,2條線之間的距離為25000/3600/8=8.68,即航向分辨率為8.68m。掃描方向分辨率可由瞬時視場(6毫弧度)和飛行高度(4500米)給出,即(6/1000)×4500=27m。由此可見,圖像上每個象元不是呈正方形,而是呈長方形。長、寬分別為航向和掃描方向,長:寬=27:8.68≈3.11:1。從視覺效果來看,顯示出沿航向拉長或沿掃描方向壓縮(扁)的現象。通地沿航向的壓縮或沿掃描方向的拉伸處理(數據重采樣方法用三次樣條法),可得到27×27或8.68×8.68的正方形象元。 1.2.4 邊緣輻射校正處理
由于飛行方向、掃描角度、太陽方位角和太陽高度的幾何關系的變化及大氣的影像,掃描圖像存在邊緣輻射差異。這種差異可以通過幾何光學模型和大氣模型予以校正,但由于模型十分復雜,影響因素很多不易操作,因此,通常采用統計方法進行圖像的邊緣輻射校正。設有圖像矩陣Anm

其中n為行數(沿飛行方向),m為列數(沿掃描方向),對圖像每列求均值得到向量:

反映了整個圖像邊緣輻射差異的平均趨勢,作者以圖像中機下點的輻射量為標準對其余各列像元進行歸一化處理,得差值向量:


校正后得圖像為:

1.2.5 正切校正處理
由于掃描角度的變化,圖像中每個象元對應的地面分辨率隨著變化,這種幾何差異的校正,通常使用正切校正方法。
h NP1 P1 NP2 P2 O 掃描角度變化如圖2:

圖2 掃描角度變化示意圖
其中0為傳感器,NP1,NP2為實際掃描面,P1、P2為地面。傳感器面向正下方等角速地收集地面數據, 為瞬時視場角,h為傳感器與地面之間的距離。
正切校正算法如下:
設圖像象元為等間隔網格上(即正方形象元)的點,則在掃描方向圖像網格,寬度h× 。從機下點開始計算,第i點網格與機下點的距離為i×h× 。實際上,地面網格點與機下點的距離為h×tan(i× )。作者按照網格寬度h× 對圖像進行重采樣,設重采樣后圖像網格點i對應圖像網格點K,則有:
i·h· =h·tan(k· )
求出K:

通常K不會落在網格點上,即K不是整數,可以得出i,滿足:
≤K≤ +1
采用距離倒數插值求出K點的象元值:

校正后的圖像象元值DN(i)=DN(K)
顯然通過上述方法校正后,圖像的列數將增加(圖3)。

圖3 正切校正前后對比影像圖
(左圖為原圖像;右圖為校正后圖像)
1.2.6 鏡像變換處理
從己經處理過的4條航帶數據來看,只有第7航帶記錄錯亂,顯示出來的圖像與實際圖像呈鏡向對應關系。可用下式表達:
f(m,n)= f(M-m,n)
m=0,1,…,M-1; n=0,1,…,N-1
其中m為列數,n為行數,M、N為該圖像總的列數和總的代數,插值方法用的是三次樣條法(圖4)。

圖4 鏡像轉換前后對比影像圖
(左圖為原圖像:右圖為轉換后圖像)
1.2.7 相對反射率轉換處理
高光譜遙感圖像不同于多光譜遙感圖像的一個根本特點在于從高光譜圖像立方體里能夠提取出精確的地物光譜反射率曲線。據此,可以直接識別地物。因此,準確地提取影像中的光譜信息是圖像分析的基礎,也是成像光譜技術研究的一項重要內容。
成像光譜儀在0.4-2.5μm工作波長范圍內主要光譜能量來自地球表面對太陽光的反射。由于太陽光受大氣吸收影響,太陽光的地面輻照度隨波長變化大。因此,成像光譜儀所獲得的圖像DN值不能代表地物的反射率值,需要進行反射率轉換。
采用了2種方法進行相對反射率轉換,下面介紹其原理:
(1)統計法 這是針對沒有進行同步地面反射率測量采用的一種方便、快捷的方法,常用的有對數殘差轉換法和內在平均相對反射率法(IARR)。研究過程編制了對數殘差轉換法軟件。
設有高光譜圖像DN(I,j,b),i=1,2,…M為行數,j=1,2….N為列數,b=1,2,…k為波段數。
在i行j列有光譜量均值:

波段b的灰度均值:

現假設:
①整個圖像區域中各類成份足夠雜亂;
②對征一波段的整個區域取平均,所得的光譜與地面吸收特征無關,則可用對數殘差轉換,得相對反射率;
logRy=logDN(i,j,b)一logSij一logDb
(2)地面定標法 在飛行的同時,用野外光譜儀實地測量典型地物的光譜,與成像光譜圖像中對應地點象元的DN值相對照,求出圖像DN值到地面反射率的轉換公式。理論與實踐都表明圖像DN值與地面反射率有以下近似的關系:

其中DN為波段b給定象元的數值,p(λ)為波段b的波長λ對應的地面反射率,Ab為乘積項,反映了大氣傳輸及儀器設備的放大比;Bb為偏移項,反映了大氣輻射值及儀器的零點偏移。Ab、Bb可以通過地面定標數據及相應位置的圖像DN值來確定。
設有N類地面定標點的反射率 ,則有一組方程:
i=1,2,…,n
求此方程組的最小二乘方解得:


本次使用的數據是1988年11月初的飛行數據,但當時沒有做同步地面測量。作者于1995年11月初進行了補充測量,一方面在于了解工作區的典型地物的光譜特征,另一方面也想做一次似同步野外定標實驗。全部測量數據于1995.11.4-11日獲得,這期間天氣晴朗,采樣時間11-16點,大體與1988.11.2、4日的時相相當。以第五航帶為例,當時選取了資江水,張家礦區7號帶、煙竹河及公路交叉處的路邊白沙及附近正常花崗巖作為定標點。
表2 野外定標點及其圖像坐標表
定 標 點 |
標 品 號 |
圖像坐標點 |
梅溪鎮資江水 |
138 |
1061,1152 |
張家公路邊白沙 |
107 |
981,956 |
張家正常花崗巖 |
110 |
979,957 |
表3 各點在對應波段上的反射率值(%)和圖像DN值(0-255)表
定標點 |
1600
nm |
2143
nm |
2200
nm |
2250
nm |
2300
nm |
2330
nm |
2450
nm |
樣品號138
坐標1061,1153 |
0.187
66 |
0.008
87 |
0.001
81 |
0.001
89 |
0.055
71 |
0.025
73 |
3.859
89 |
樣品號107,
坐標號981,956 |
60
118 |
56
135 |
53
113 |
53
125 |
56
106 |
54
101 |
62
109 |
樣品號110,
坐標979,957 |
44
88 |
41
107 |
39
90 |
40
100 |
40
87 |
39
93 |
51
96 |
這樣依此3點可以得到一條回歸曲線,可將圖像上任一點的DN值轉化為相對反射率值,并編制了相應的軟件。上述3點轉換后的相對反射率值(表4)如下:
表4 轉換后的相對反射率值表
坐標點 |
1600
nm |
2143
nm |
2200
nm |
2250
nm |
2300
nm |
2330
nm |
2450
nm |
1061 1153 |
3.568 |
3.278 |
4.441 |
4.092 |
2.828 |
0.050 |
7.351 |
981 956 |
69 |
65 |
67 |
66 |
54 |
54 |
77 |
979 957 |
31 |
29 |
22 |
30 |
39 |
39 |
32 |
2. 波譜特征
野外地物波譜測量是遙感應用的基礎,而高光譜分辯率數據具有圖譜合一的特征,因此波譜測量對航空高光譜數據的處理和解釋尤其重要。一般情況下,在飛行時應同步進行野外地物波譜測量,為野外定標及數據處理和解釋提供依據。但由于條件限制,如前所述飛行時沒有同步進行野外地物波譜測量。為了彌補這方面的不足,研究過程中擬同步性地補充進行了野外地物波譜測量。在EVNI軟件上建立了野外實測光譜曲線庫,并研究了該區典型地物的波譜特征。
2.1 地物波譜測量及影響巖石反射光譜特征的主要因素
2.1.1 野外地物譜測量
利用美國地球物理環境研究公司(GER)研制的IRIS紅外智能波譜儀對苗兒山地區的典型地物進行了波譜測試,共獲取了150多條波譜曲線,基本上反映了研究區各典型地物的波譜特征,并為野外定標提供了基本數據。
野外光譜測量是在1995年11月4日至11月11日進行的。如前所述,這個時間段與當年獲取數據的時間相當,便于地面所測光譜信息與航測的光譜信息做對應分析。在測錄中用GPS在1:5的地形圖上做了較準確的定位,并且對每一個測點的地物類型做了細致的描述,使以地面波譜測量為基礎的處理比較可靠。為了分析測得的光譜曲線,先分析一下影響巖石反射光譜特征的主要因素。
2.1.2 影響巖石反射光譜特征的主要因素
⑴巖石成分的影響 由于巖石的成因不同,巖石的成分也隨之發生變化,巖石呈現出不同顏色。例如,酸性中花崗巖以長英質礦物為主,巖石呈淺色調;基性巖中的輝長巖以鐵鎂質礦物為主,巖石呈暗色調。前者的光譜反射率高于后者。從基性巖到酸性巖,它們之間的反射率可在70%范圍內變化,深色巖石的反射率要低于淺色巖。可見巖石光譜的反射率的高低與其顏色有關,更確切地說應與巖石的成份有關。
⑵巖石結構、構造的影響 巖石在生成過程中由于物理化學條件不同、生成不同程度、不同結晶度的巖石。巖石的結構、構造反映巖石中不同礦物顆粒大小及它們之間的關系。通常,巖石的結晶度和粒度不同造成同類巖石光譜反射率的不均勻性,如粗粒花崗巖的光譜反射率變化要大于細粒花崗巖。
⑶蝕變和風化作用的影響 在內營力和外營力作用下,巖石發生蝕變或風化,生成蝕變和風化礦物,并引起巖石結構、構造和色調等的變化。
2.2 苗兒山地區典型地物波譜庫的建立及地物特征分析
2.2.1 典型地物波譜庫的建立
從野外獲取的150多條光譜曲線按類別可分成11個大類,并以此建成11個庫。它們分別是正常花崗巖、花崗巖風化而成的沙土、花崗巖弱蝕變破碎物、花崗巖強水云母化、花崗巖與草混合物、礦石、硅化帶、稻田、水體、各類植被、混合植被等。
為了使實測曲線(0.4-2.5um)能與從高光譜圖像立方體中提取的曲線(1.6-2.45um)進個對比分析和匹配研究,須將原曲線按高光譜圖像的七個波段及相應的帶寬(100nm或50nm)進行重采樣。 2.2.2 地物波譜特征分析
根據建成的庫,對其波譜特征進行了分析。取各庫的平均光譜曲線,如S1、S2、S3、S4,分別代表正常花崗巖、花崗巖風化而成的砂土、花崗巖弱蝕變破碎物、花崗巖強水云母化等。
S1代表正常花崗巖 總體反射率較低,在20%左右,曲線較平坦。
S2代表風化后成的沙土 總體反射率在40-500%,在2200處有較明顯的吸收峰。
S3弱蝕變花崗巖 總體反射率比S2稍低,在2200nm處有較明顯的吸收峰。
S4強蝕變花崗巖 總體反射率明顯升高,特別是在低波段區間反射率升高更為明顯,在2200nm處有很明顯的吸收峰,880nm處有較明顯的吸收峰。
S5花崗巖型鈾礦石 總體反射率同正常花崗巖S1,在20%左右,但不同的是在2200nm處有明顯吸收峰。
S6花崗巖與各種草混合后的曲線 在670-760nm之間有明顯的陡坡,反映植被的典型特征,即所謂的“葉綠素紅邊”。植被所占比例越多,這條斜邊越長,反之,越短。
S7硅化帶的特征曲線 總體反射率<40%,在670nm、2200nm有吸收峰(圖5-6)。

圖5 硅化帶的光譜特征曲線
圖6 重采樣后硅化帶的特征曲線
S8稻田反射率特征曲線 其中,S81與植被曲線(S10)有非常近似的特點,在670-760nm之間有明顯的陡坡,在1500nm處有明顯吸收峰,在1650nm和2250nm處有明顯反射峰,反映了田中稻草梗的反射率特點;S82反映的主要是田中泥土的反射率特征,由于它們源于花崗巖風化破碎物,故近似于花崗巖特征曲線。
S9水體 資江水,梅溪水的特征曲線,反射率低。
S10各種植被及混合植被的特征 由于葉綠素在670nm紅波段具強吸收,在760nm近紅外波段具高反射的特點,便得各種植被在670-760nm之間具有明顯的陡坡,即所謂的“葉綠素紅邊”。它是健康植被的一個特征,如果植被或農作物由于受到如重金屬、碳氫化合物、干旱、病蟲害、缺氮等因素引起脅迫(stress),則會使這些植被或農作物的特征曲線產生微小的變化,即所謂的“紅邊藍移”現象。利用“紅邊藍移”現象可以預測病蟲害等的發生,并及時采取補救措施。本區的某些植物,如乍子樹、竹葉、松樹、雜草等都具有此種現象(圖7),說明上述植物受到脅迫。

圖 7 植物受到脅迫后的光譜特征曲線
(上圖為乍子樹、竹葉、松樹;下圖為雜草)
3 高光譜數據的地質應用及效果
3.1 高光譜分辨率圖像上硅化帶的識別及效果
高光譜分辨率數據的一個最重要用途是識別蝕變礦物。盡管由于原始數據和研究區植被覆蓋方面的原因,使應用效果受到了一定的限制,但在硅化帶的識別方面仍取得了成功。
圖8為第5航帶第2波段(中心波長2.143μm)圖像(已經各種處理和精校正),從該圖像上可以清楚地看到許多近東西向的暗色調線性體,這些線性體在對應區域的TM圖像上沒有或僅略有顯示,在原有的大比例尺地質圖上也只有少數幾條被填繪出來。通過對比圖像上直接提取的已知帶(硅化斷裂帶)和未知帶的光譜曲線,發現兩者具有一致的曲線特征。因此,推測這些未知的線性體可能是硅化斷裂帶。
作者對分析結果進行了野外檢驗,證實它們大多數確系硅化斷裂帶,其中馬家南邊的那條硅化帶規模很大,寬10~20m,主要為白色石英脈帶(鈾含量為26 x10-6~45義10—6),帶內有殘留碎裂花崗巖(鈾含量40x10-6,具紅化、高嶺土化和弱硅化),圍巖為碎裂花崗巖(鈾含量為60 x10-6一70 x10-6)。硅化帶內殘留碎裂花崗巖的鈾有丟失,而在圍巖中有富集現象。在豆詐山巖體及其周圍的高光譜圖中,解譯出的線性構造包括了該區地質圖上填出的主要硅化斷裂帶,其中有些硅化斷裂帶直接控制了鈾礦化。

圖8 第5航第2波段(2.143μm)圖像
3.2 高光譜分辨率圖像與其他多源信息的綜合分析及效果
將高光譜數據、TM數據與已有地質圖等多源信息進行融合,可形成一幅集多種信息于一體的新類型圖像(圖9)。

圖9 綜合解譯圖
在該圖像上,地質圖上已有的主要斷裂構造和硅化帶都已被解譯出來,同時,還解譯出很多線狀構造和環狀構造,而這些構造在地質圖上是沒有的。如fl—沙子江-茶坪斷裂構造,f2一斷裂構造,f3一斷裂構造,f4一牛大江硅化帶等;環狀構造有4個,分別以C1、C2、C3、C4表示。 圖9 綜合解譯圖 其中以Cl規模最大。本區是巖漿巖分布區,環狀構造的形成一般與巖漿的侵人活動有關。4個環狀構造都發育在兩期巖體的交界部位,其中C2和C3是豆詐山巖體與香草坪巖體之接觸部位,C4是牛塘坳巖體(γ53-1)、香草坪巖體和苗兒山巖體之接觸部位,可能為更晚期的巖漿侵人活動所致。C1不僅規模最大,且結構復雜,該環形構造西南部處于豆詐山巖體與香草坪巖體之接觸帶;東北部受香草坪巖體、苗兒山巖體和牛塘傲巖體三者的影響;東南部基本上是香草坪巖體和苗兒山巖體的分界線;西北部分還同時與小環C3相交。此環形構造還是本區四條大硅化斷裂帶(F1, F2, F31, F32)的夾持部位。因此,該環形構造區具有很好的構造、鈾源(γ51、γ52-2)和成礦的先決條件,雙滑江礦床和白毛沖礦床就在此環形構造范圍內。前者產于該環形構造的環邊上。本次研究預測的一處鈾礦化遠景區(P1)就在該環的西北部分。由此可以認為,該環形構造的展布范圍應該是很有成礦前景的地區。
經過上述各種處理及綜合解譯,結合前人對本區成礦機理的認識,本次研究預測了三個成礦有利區,通過驗證,其中的兩個地區放射性測量值偏高,并有鈾礦化顯示。
參考文獻
1.鄭蘭芬 王晉年,1992,成像光譜遙感技術及成像光譜信息提取的分析與研究,《環境遙感》,Vol. 7,No. 10。
2.王晉年鄭蘭芬,1996,成像光譜圖像光譜吸收鑒別模型與礦物填圖研究,《環境遙感》》,Vol.11,No. 10。
3.王向軍,1995,成像光譜信息可視化表達系統,中國科學院遙感應用研究所碩士學位論文。 4.王欽敏,1995,衛星遙感技術發展和應用動態,《遙感信息》,Vol.40, No. 4。
5.R. J. P. Lyon, 1996,風化及其它類荒漠漆表面層對高光譜分辨率反射率的影響(一),《環境遙感》,Vol. 11,No. 2。
6.Joseph C.H, Chein-I Chang, 1994, Hyperspectral Image,Classification and Dimenesionality Reduction:An Orthogonal Subspace Projection Approach,《IEEE Transaction on Geoscience and Remote Sensing》,Vo1.32, No.4(779-784)。
7.Geotz, A.,and Herring, M:,1989, The High Resolusion Imaging Spectrometer (HIRIS)for EOS,《IEEE Transaction on Geoscience and Remote Sensing》,Vol.27, No.2。
8.Mazer, A.S.,et al, 1988, Image Processing Software for Imaging Spectrometry Data Analysis.《Remote Sensing of Environment》,Vol.24,No. l(201-210)。
9.William H.F.,et al,1994, Retrieval of Apparent Surface Reflectance from AVRIS data:A Comparison of Empirical line, Radiative Transfer, and Spectral Mixture Methods,《Remote Sensing of Environment》》,Vo1.47(311-321)。
10.ENVI User's Guide, 1996,Research Systems, Inc.,USA
11.A Hyperspectral Remote Sensing Primer, 1996, NASA, USA。
|