熊文成,申文明,王 橋,史園莉,肖如林,付 卓
(環境保護部衛星環境應用中心,北京 100094)
我國經濟社會的快速發展帶來了巨大的環境問題,如水體污染、城市垃圾、大氣污染、生態破壞等。現在環境問題已經嚴重制約了我國經濟社會的發展。要解決環境問題,首先必須了解我國的環境狀況,這就需要對環境質量進行有效監測。環境一號衛星是我國自主研制的小衛星星座,其目標之一就是對我國環境質量狀況進行嚴密監測[1]。綜合有效地利用環境一號衛星各載荷數據,對實現我國生態破壞、環境污染進行大范圍、全天候、全天時動態監測,對生態環境質量變化的過程和趨勢的預測都有著重要意義。要實現對各載荷數據的綜合應用,業務化地數據幾何精校正是一個重要的基礎[2~4]。
一般來說,圖像的幾何精校正過程可以分為地面控制點選取、模型糾正、重采樣3個步驟。傳統的控制點選取方式主要是人工在圖像上有明顯的、清晰的定位識別標志處選點,如道路交叉點、河流叉口、建筑邊界、農田界線[4]。自動選點則主要是基于基準影像和配準影像相似性關系,自動匹配選點的過程。自動匹配在減少人工選點的主觀性以及提高工作效率方面有著重要意義。自動匹配的方法在國內外研究較多,基本可以分為基于灰度的匹配和基于特征的匹配[5,6]。基于灰度的圖像匹配算法主要有相關系數法、互信息法、傅里葉變換法、最小二乘法[7~11];基于特征的圖像匹配算法主要有基于角點檢測的圖像配準方法和基于圖像特征描述符的配準方法。角點檢測方法有Harris角點檢測算法和Susan角點檢測算法;常用的特征描述符圖像配準算法有尺度不變特征變換匹配(SIFT)算法和加速穩健特征(SURF)算法[12~15]。轉換模型(校正方法模型)要根據遙感圖像幾何畸變的性質來確定。現有的糾正模型主要有物理模型、有理函數模型、多項式擬合或三角網擬合等[2]。重采樣主要有3種方法,即最鄰近點法、雙向線性插值法、三次褶積法。
為實現環境一號衛星數據業務化運行,本文從影像自動配準的角度,研究環境一號A/B星數據的自動幾何精校正技術的設計與實現。
2008年9月6日,環境一號衛星“2+1”星座中的兩顆光學小衛星(A星和B星)“一箭雙星”發射成功。環境一號衛星A星搭載了兩個charge-coupled device(CCD)相機(CCD1、CCD2)及一個高光譜成像儀(HSI);環境一號衛星B星搭載了兩個CCD相機及一個紅外相機(IRS)。具體參數可見文獻[1]。
環境一號衛星A/B星有以下特點:a.A星與B星上的4臺CCD相機參數基本一致,這意味著4臺CCD相機的圖像具有良好的相關性;b.A星上CCD與高光譜成像儀是同平臺不同載荷,B星上CCD與紅外多光譜相機是同平臺不同載荷,同平臺數據具有較好的幾何相關性,具體關系可見圖1、圖2。

圖1 環境一號衛星A星2級數據產品幾何關系圖Fig.1 Geometric relationship of HJ-1A satellite

圖2 環境一號衛星B星2級數據產品幾何關系圖Fig.2 Geometric relationship of HJ-1B satellite
環境一號衛星地面處理分系統向外分發的環境一號衛星數據主要為2級產品,2級產品為經過系統幾何校正和相對輻射校正的產品。由于環境一號衛星的小衛星平臺、大幅寬等特點,其數據幾何畸變現象較為突出。本文從環境一號衛星2級數據的外部定位精度及內部幾何畸變分析環境一號衛星數據幾何特點及變形規律。
定位精度是指衛星產品的幾何定位精度,即2級產品圖像上地理位置和真實位置之間的差異。外部定位精度評價以LandsatETM幾何精糾正數據產品為參考圖像,在環境一號衛星的2級圖像以及參考影像上,選擇多個均勻分布的匹配同名點,分別獲取環境一號衛星數據和參考數據的圖像坐標,然后計算出匹配點(環境一號衛星圖像坐標和參考圖像坐標)的差值。通過計算多景不同地區(東北、西南、西北、中部)2級圖像的位置誤差值,并把總的均方根誤差作為環境一號衛星2級數據的定位誤差。經多景統計,環境一號衛星A/B星的幾何定位精度評價結果見表1、表2。

表1 環境一號衛星A星幾何定位精度評價結果Table1 Geometricalpositioning accuracy evaluation resultsof HJ-1A satellite m

表2 環境一號衛星B星幾何定位精度評價結果Table2 Geometricalpositioning accuracy evaluation resultsof HJ-1B satellite m
綜合來看,環境一號衛星A/B星的幾何定位精度大約為1 100m。
遙感圖像的幾何形變是指圖像上像元在圖像坐標中的坐標與其在地圖坐標系統等參考系統中的坐標之間的差異[2]。形變可由傳感器結構等因素引起,如攝影機的焦距變動、像主點偏移、鏡頭畸變等;也可由遙感器以外的各因素所造成的誤差引起,如傳感器的外方位變化、傳感器介質的不均勻、地球曲率、地形起伏、地球旋轉等因素所引起的誤差。
圖像內部幾何變形可以歸納為長度變形、角度變形和放射變形等指標,評價的內容包括變形的絕對量和整幅圖像變形的一致性。具體的方法如下所述。
1)從圖像上選擇同名控制點20~30對,并從4個方向(垂直軌道、沿著軌道和對角方向)大致進行配對與分組。每組類別中有20組左右的配對點。
2)計算圖像上點對距離與地圖上各同名點對距離的相對差距。點對距離是指圖像上兩點的距離,即

式(1)中,(L1,P1)和 (L2,P2)分別為圖像上兩個點的坐標;Δp為圖像分辨率。圖像上點對距離diimg與地圖上同名點對距離diref的相對差距為

統計每組配對點間的距離差的均值作為變形誤差,即

3)計算8個方位的控制點到景中心的距離變化,得到相對變化程度和矢量方向。計算方法與上一步類似。
將4個方向控制點間相對幾何誤差以及各方位控制點相對于景中心的幾何定位誤差取平均得到,環境一號衛星A/B星的內部幾何相對誤差是0.12%~0.17%,絕對誤差是100~400m。
從分析結果來看,環境一號衛星數據內部幾何誤差較大,從各個方向的統計值來看,也沒有展示出明顯的規律。這對于進行下一步的幾何精校是一個挑戰。
基于環境一號衛星A/B星的參數特點,幾何定位以及內部形變等特點,本文設計了環境一號衛星A/B星自動幾何精校正的技術流程。
環境一號衛星多載荷數據總體幾何校正流程設計如圖3所示。

圖3 環境一號衛星多載荷數據總體幾何校正流程Fig.3 Geometric correction of HJ-1 satellite multi-sensor data
總體流程的設計基于以下3個方面。
1)利用同平臺的高分辨CCD對高光譜數據實現幾何自動匹配。高光譜成像儀在光譜復原后任一波段的切片數據相當于一幅多光譜CCD相機單波段數據,兩者數據的輻射機理是相同的。同星平臺獲取的CCD和高光譜數據具有較好的相似性。
2)利用同平臺的高分辨CCD對紅外相機數據實現幾何自動匹配。CCD數據第四波段和紅外相機第一波段,波長較為接近,輻射機理相同。且同星平臺的CCD與紅外相機數據具有較好的相似性。
3)利用環境一號衛星CCD精校正影像作為控制庫。相對于基于控制點影像庫的方法,控制影像易于維護以及匹配更多的點,而控制點庫可能由于地物的變化而需要進行更新,并且影像局部受云雨污染時難以進行匹配。另外,同源的環境星CCD較異源影像(如ETM/TM控制庫)具有更好的匹配性,自動匹配的精度和準確性更高。
因此本流程的關鍵是對環境一號衛星CCD影像實現自動幾何精校正。
實現CCD影像的自動幾何精校正的關鍵也就是實現待糾正的CCD影像與CCD控制庫影像的自動匹配。
基于灰度和特征的圖像配準方法各有優缺點[5]。基于灰度的圖像配準算法簡單,計算量大,耗時長,對光照條件敏感。當待匹配點位于低反差區內,即在該窗口內信息貧乏,信噪比很小,基于灰度的匹配算法成功率不高。而基于特征的圖像配準算法可以在一定程度上克服此缺點,但在精度上遜于基于灰度的圖像匹配算法。因此,環境一號衛星數據的校正將結合兩種方法的優勢,以達到快速、準確的配準效果。
環境一號衛星CCD數據匹配策略利用了基于灰度的匹配和基于特征的匹配相結合的方式。利用金字塔影像結構,將上一層影像的特征匹配結果傳到下一層作為初始值,并考慮對粗差的剔除或改正。最后以特征匹配結果為“控制”,對其他點進行匹配或內插。由于基于特征匹配是以“整像素”精度定位,因而對需要高精度的情況,將其結果作為近似值,再利用最小二乘影像匹配進行精確匹配,取得“亞像元”精度[15,16]。環境一號衛星配準策略流程見圖4,可以分為以下幾個步驟。

圖4 環境一號衛星配準策略流程Fig.4 CCD registration flow of HJ-1 satellite
1)建立金字塔分層影像。環境一號衛星CCD數據特征匹配方案應用了金字塔分層影像數據。金字塔影像的建立按照2×2像元變換成一個像元,此時,上一層的結果與下一層2×2個像元的公共角點相對應。利用三級金字塔逐級精化,利用頂級金字塔數據量小的特點進行快速初始定位,利用頂級金字塔初始定位引導中級金字塔進行初始定位精化,利用中級金字塔精化的初始定位引導底層金字塔或原始影像精匹配。
2)特征提取。在特征點檢測階段優先采用Fast、Harris、Susan等速度較快的角點檢測算法。本研究采用的是Fast特征點檢測算法。
3)粗差剔除。在一個小范圍內利用二次曲面為模型進行擬合,將殘差大于某一閾值的點作為粗差剔除。當所有錯誤的匹配點作為粗差剔除后,即得到與目標模型一致的匹配點對。有效的粗差剔除決定了影像配準是否成功,粗差剔除是匹配關鍵關節。必須采用穩健的粗差剔除方法保證配準成功。本研究考慮小概率事件原理粗差剔除,按照5%小概率事件發生原理,可認為超過中誤差2倍的匹配點為粗差,進行迭代剔除可達到粗差剔除目的。
4)校正模型。針對環境一號衛星2級數據,幾何校正模型可采用多項式校正模型。遙感圖像的總體變形是圖像的平移、縮放、旋轉、仿射、偏扭、彎曲以及更高次的基本變形的綜合作用結果,這樣的變形用一個嚴格的數學表達式來描述是困難的,可以用一個適當的多項式來描述校正前后圖像相應點之間的坐標關系。
5)重采樣。經空間變換后輸出的新圖像像元,在多數情況下會落在原始圖像陣列中的幾個像元之間(即共軛位置),因此輸出圖像的像元灰度值,必須通過適當的方法把該點四周鄰近的若干個整數點上的像元灰度值對該點的灰度值貢獻累積起來進行計算,這個過程稱為數字圖像的重采樣。校正空間網格上的亮度值等于原始圖像空間共軛點的亮度值,即

式(4)中,(ε,η)為校正空間某點坐標;(X,Y)為該點在原始圖像空間上的共軛點坐標。
通常有3種方法確定共軛點的亮度值f(X,Y),即最鄰近點法、雙線性插值法,三次褶積法。
環境一號衛星紅外相機數據和高光譜數據以同步成像的CCD數據進行匹配校正。由于是同步成像,匹配的成功率一般較高。
4.3.1 高光譜與CCD圖像配準
在特征點檢測階段優先采用Fast、Harris、Susan等速度較快的角點檢測算法。以CCD圖像為基準地圖,特征匹配階段建立2層金字塔。并且,由于CCD與高光譜同時搭載在A星上,受外方位元素引起的形變以及地形引起的形變具有很大的相似性,兩個載荷之間的幾何差異主要是平移和旋轉等差異。因此,校正模型采用一次多項式即可。參與匹配的波段盡量采用CCD第四波段和高光譜數據波長類似的波段。通過配準,其配準精度可以小于1個像素。
4.3.2 紅外與CCD圖像配準
在特征點檢測階段優先采用Fast、Harris、Susan等速度較快的角點檢測算法。以CCD圖像為基準地圖,特征匹配階段建立2層金字塔。并且,由于CCD與紅外相機同時搭載在B星上,因此受外方位元素引起的形變以及地形引起的形變具有很大的相似性,兩個載荷之間的幾何差異主要是平移和旋轉等差異。因此,校正模型采用一次多項式即可。參與匹配的波段盡量采用CCD第四波段和紅外相機的第一波段,因為這兩個波段之間的相關性最高。通過配準,其配準精度可以小于1個像素。
根據以上方法、流程設計,環境一號衛星數據處理系統已經實現了環境一號衛星數據的自動幾何精校正。本研究在全國范圍內各個方向上選取10景影像。利用自動幾何精校正方法進行自動匹配選點、校正,并進行精度結果抽樣檢驗,結果見表3,典型效果見圖5~圖7。

表3 多項式幾何精校正Table3 Polynom ialgeometric precision correction

圖5 CCD糾正細節Fig.5 Correction detailsof CCD image

圖6 CCD與高光譜配準影像細節Fig.6 CCD and HSI im age registration

圖7 CCD與紅外配準影像細節Fig.7 CCD and infrared image registration
本文分析了環境一號衛星2級數據(系統幾何校正后的數據)的幾何特點。系統幾何校正后的環境一號衛星數據外部定位精度大概是1 km;另外,由于環境一號衛星平臺為小衛星平臺,并且圖幅的幅寬很大,幾何變形規律性不強。其內部畸變也可以達到10多個像素,這給準確的幾何精校正帶來很高的難度。
基于環境一號衛星數據的特點,本文設計了自動幾何校正方法流程。對于CCD影像,筆者利用自動匹配的多項式方法進行校正,其優勢是快速且成功率較高,缺點是對部分地區(尤其是山地)的校正精度有限。對于紅外相機和高光譜數據的幾何校正,由于其與同星CCD的高度相關性,可以利用同星CCD進行自動匹配校正。通過實際數據評測,結果表明,環境一號衛星CCD經自動幾何精校正后,圖像精度可從1 km改善到2個像素(60m)以內。同星不同載荷圖像之間自動匹配方法可以獲取小于1個像素的圖像匹配精度。
總之,基于自動匹配的選點和精校正對環境一號衛星數據大批量處理具有重要意義,一方面可以減少人工選點的誤判而帶來的誤差,另一方面可以大大節省人力的投入。
[1] 王 橋,魏 斌,王昌佐,等.基于環境一號衛星的生態環境遙感監測[M].北京:科學出版社,2010.
[2] 朱述龍,張占睦.遙感圖像獲取與分析[M].北京:科學出版社,2000.
[3] 張永生.遙感圖像信息系統[M].北京:科學出版社,2000.
[4] 趙英時.遙感應用分析原理與方法[M].北京:科學出版社,2003.
[5] 余 濤,王 橋,魏 斌.環境一號衛星遙感數據處理[M].北京:科學出版社,2013.
[6] Barbara Zitova,Jan Flusser.Image registration methods:A survey[J].Imageand Vision Computing,2003,21:977-1000.
[7] Li Xiaoming,Zhao Xunpo,Zheng Lian,et al.An image registration technique based on Fourier-Mellin transform and its extended applications[J].Chinese Journal of Computers,2006,29(3):116-122.
[8] Castro ED,MorandiC.Registration of translated and rotated imagesusing finite Fourier transform[J].IEEE Transactions on Pattern Analysisand Machine Intelligence,1987(9):700-703.
[9] Viola P,WellsW.Alignment by maximization ofmutual information[J].International Journal of Computer Vision,1997,24(2):137-154.
[10] Thevenaz P,Unser M.A pyram id approach to sub-pixel image fusion based on mutual information[C]//Proceedings of the IEEE International Conference on Image Processing.Lausanne,Sw itzerland,1996,1:265-268.
[11] Thevenaz P,Unser M.An efficientmutual information optim izer for multiresolution image registration[C]//Proceedings of the IEEE International Conference on Image Processing.Chicago,Illinois,USA,1998:833-837.
[12] Goshtasby A.Description and discrimination of planar shapes using shapematrices[J].IEEE Transaction on Pattern Analysis and Machine Intelligence,1985,7(6):738-743.
[13] Goshtasby A.Image registration by local approximation methods[J].Imageand Vision Computing,1988,6:255-261.
[14] Danip P,Chaudhuris.Automated assembling of images:Image montage preparation[J].Pattern Recognition,1995,28(3):431-445.
[15] 章毓晉.圖像工程:下冊[M].北京:清華大學出版社,2004.
[16] 李 峰,周源華.采用金字塔分解的最小二乘影像匹配算法[J].上海交通大學學報,1999,33(5):513-519.