方韓康,張波,陳衛(wèi)榮,吳樊,王超
(1 中國科學院空天信息創(chuàng)新研究院 中國科學院數(shù)字地球重點實驗室, 北京 100094;2 中國科學院大學資源與環(huán)境學院, 北京 100049; 3 中國資源衛(wèi)星應(yīng)用中心, 北京 100094)
高分三號衛(wèi)星是中國首顆C波段高分辨率合成孔徑雷達(synthetic aperture radar, SAR)衛(wèi)星。它具有12種不同的工作模式,在聚束模式下可達到1 m分辨率,為海洋監(jiān)測[1-2]、建筑提取[3-4]、建筑物損毀評估[5]、船舶檢測[6]、海面風場反演[7]等不同應(yīng)用提供了數(shù)據(jù)支持。高分三號發(fā)布的L1A級產(chǎn)品是過采樣復(fù)數(shù)據(jù)產(chǎn)品,因其含有相位、強度、極化信息,而被廣泛訂購使用。L1A數(shù)據(jù)集中包含*.meta.xml、*.rpc、*.incidence.xml等輔助數(shù)據(jù)文件。這些輔助數(shù)據(jù)文件集的提供,為用戶進行自定義的更高級別圖像產(chǎn)品形式的生產(chǎn)提供了方便。
SAR影像產(chǎn)品處理方式和等級不同,其像素值可由多種物理量表示:DN(digital number)值、幅度值、強度值、后向散射系數(shù)(sigmma nought)值[8]。DN值為無符號整型數(shù)值,常用于SAR影像復(fù)數(shù)據(jù)的短有效位存儲;幅度值與強度值對應(yīng)SAR影像像素復(fù)數(shù)據(jù)的模及模平方;后向散射系數(shù)值對應(yīng)SAR影像強度值經(jīng)輻射糾正后的數(shù)值表達。DN值、幅度值、強度值、后向散射系數(shù)值四者之間為指數(shù)或?qū)?shù)的變換關(guān)系,因此它們之間可相互轉(zhuǎn)換并為正相關(guān)關(guān)系。
經(jīng)輻射糾正后,SAR圖像像素值常用后向散射系數(shù)表示,其大小對應(yīng)被照射地物在該波段上對入射雷達波反射能力的強弱?;诤笙蛏⑸湎禂?shù)圖,同一傳感器在不同時刻獲取的圖像,以及同一時刻成像場景中不同位置處的地物后向散射特性均可做定量化對比,以便于基于后向散射系數(shù)實現(xiàn)生物量參數(shù)反演等。這些定量對比與參數(shù)反演的應(yīng)用就要求SAR圖像輻射糾正要有較高的精度與準確度。
針對高分三號衛(wèi)星影像單極化圖像的絕對輻射強度糾正、全極化數(shù)據(jù)的通道不平衡和串擾均有學者進行了實驗研究。Chang等[9]指出通常情況下高分三號衛(wèi)星通道強度不平衡優(yōu)于0.5 dB,通道相位差優(yōu)于10°,串擾優(yōu)于-35 dB。對少數(shù)殘差較大的圖像,通過進一步極化糾正,可以將極化通道不平衡降低到0.26 dB,通道相位差降低到±0.2°,串擾精度優(yōu)于-42 dB。Chen等[10]利用三面角角反射器在河北地區(qū)做了高分三號的同步散射強度測量和圖像質(zhì)量評價,依據(jù)實測結(jié)果發(fā)現(xiàn)高分三號衛(wèi)星的等效噪聲系數(shù)為-25 dB,稍微優(yōu)于Radarsat-2影像的-24.6 dB;高分三號衛(wèi)星絕對輻射定標值相較于Radarsat-2要低3~5 dB。高分三號不同模式產(chǎn)品對應(yīng)的標稱輻射糾正參數(shù)均包含在L1A數(shù)據(jù)集中,普通用戶可依據(jù)這些參數(shù)自行進行后向散射系數(shù)圖生成,但采用何種方式實現(xiàn)多視處理、克服L1A級產(chǎn)品量化過程中的空值與零值干擾,獲取糾正結(jié)果圖的高精度輸出仍值得研究。
星載SAR幾何糾正是消除衛(wèi)星成像時由于太陽同步軌道面傾斜、地球自轉(zhuǎn)以及局部地形等因素造成的影像幾何畸變的過程[11]。已經(jīng)發(fā)展了多種星載SAR影像幾何糾正方法,大致可分為嚴密的物理模型方法以及通用的幾何模型。張過等[12-13]通過選取地面控制點(ground control point, GCP)建立了通用的有理多項式(rational polynomial coefficient, RPC)模型,并驗證了RPC模型具有對傳統(tǒng)SAR嚴密物理模型即距離多普勒(range-Doppler, RD)模型的可替代性。中國高分三號衛(wèi)星數(shù)據(jù)在斜距方向與方位時間的標定精度有限。呂冠男等[14]選取京津冀地區(qū)高分三號影像數(shù)據(jù),探究影像的幾何定標后影像定位精度。實驗證明,高分三號影像無控制點幾何校正平面誤差范圍為109.428~116.290 m。在高分三號SAR影像L1A數(shù)據(jù)集中的*.rpc文件中提供了可進行幾何糾正的RPC模型參數(shù)。根據(jù)高分三號衛(wèi)星提供的用戶手冊,用戶可以只通過影像對應(yīng)的*.rpc文件,利用最新版的Envi5.5軟件對影像進行幾何校正。然而,Envi5.5軟件在對部分高分三號影像進行幾何糾正時,存在無法完成幾何糾正或輸出結(jié)果出錯的問題。
為實現(xiàn)從L1A級到更高等級產(chǎn)品處理結(jié)果的高可靠和高精度,提出一套完整的處理流程用以克服L1A級圖像產(chǎn)品過采樣問題、輻射量化過程中出現(xiàn)空值與零值帶來輻射糾正結(jié)果偏差的問題、以及RPC參數(shù)文件在Envi軟件中出現(xiàn)失效導致幾何糾正無法進行的問題等。并通過與同波段的Sentinel-1 SAR影像、Sentinel-2光學影像進行對比證明了本文方法的有效性。
高分三號衛(wèi)星L1A級數(shù)據(jù)產(chǎn)品中的柵格圖像是投影方向為斜距方向上的單視復(fù)數(shù)據(jù)(single look complex, SLC),其像素值為量化后的DN值。在遙感綜合應(yīng)用中,光學影像、DEM數(shù)據(jù)等均為投影面為地面的柵格圖。因此,在后續(xù)處理中也須對L1A級產(chǎn)品進行輻射強度與幾何投影的標準化。鑒于此,提出對L1A級產(chǎn)品進行輻射與幾何糾正的一體化處理流程如圖1所示。
在該處理流程的輻射糾正中,依據(jù)等效噪聲系數(shù)對SAR圖像中各像素值進行檢核,以剔除SAR圖像量化中空值與零值對定標結(jié)果的影響。在幾何校正中提出RPC反算方法,以應(yīng)對高分三號元數(shù)據(jù)文件提供的角點坐標與RPC參數(shù)無法自洽的問題。在后向投影采樣輸出中將傳統(tǒng)的Lee濾波融合到采樣過程中,以替代傳統(tǒng)的多視處理或雙線性插值處理,抑制SAR圖像斑點噪聲的影響。其各要點方法分解描述如下。
(1)
1.2.1 基于RPC模型參數(shù)的幾何糾正
高分三號SAR衛(wèi)星L1A級產(chǎn)品在元數(shù)據(jù)文件集中提供了以(*.rpc)命名的RPC參數(shù)文件。直接利用該參數(shù)文件,可以實現(xiàn)地面上一點P(B,L,H)到圖像上對應(yīng)像點坐標P′(r,c)的轉(zhuǎn)換。其中(B,L,H)代表該點的緯度、經(jīng)度、高程;(r,c)代表圖像上的行列坐標。RPC參數(shù)模型及其中參數(shù)的含義見文獻[12]。
為方便利用RPC參數(shù)實現(xiàn)幾何糾正過程,常采用如下所述的后向投影過程予以實現(xiàn):1)求出待輸出幾何糾正結(jié)果圖幅框4個角點的地理坐標,2)采用RPC模型建立從糾正結(jié)果圖P(B,L,H)到原圖P′(r,c)的后向投影關(guān)系,3)設(shè)定輸出結(jié)果中各像素間距(輸出步長)后,在原圖中進行采樣實現(xiàn)幾何糾正結(jié)果的輸出。
在上述糾正過程中的第一步尤為重要,圖幅框4個角點的像素坐標可通過RPC模型反算原圖的4個角點經(jīng)緯度坐標P1(BTL,LTL),P2(BTR,LTR),P3(BBL,LBL),P4(BBR,LBR),而后通過求取最大外接矩得到校正后影像范圍。4個角點的經(jīng)緯度坐標可通過圖1中RPC幾何校正模塊所示的反算迭代流程予以實現(xiàn),其過程概括如下:
1)設(shè)定計算初始值為影像中心經(jīng)緯度坐標(B0,L0);
4)將像素坐標差(Δr,Δc)轉(zhuǎn)換成為投影平面坐標差(Δx,Δy);
8)輸出角點經(jīng)緯度坐標Pcorner(B,L)為最終解。
在上述RPC模型反算過程中有以下兩個關(guān)鍵點說明如下:1)解算起始點參數(shù)的設(shè)置對反算結(jié)果的計算效率影響較大,起算點設(shè)置不合理會導致算法不收斂。RPC參數(shù)文件中提供的經(jīng)緯度偏移參數(shù)(Lat_Offset, Lon_Offset),其物理含義為景中心經(jīng)緯度,將其作為初始值可以提高迭代速度,且保證算法收斂。2)L1A級別產(chǎn)品的柵格影像存儲順序以SAR傳感器接收到回波的先后次序進行存儲,因此在實現(xiàn)從像素坐標差到大地坐標差的轉(zhuǎn)換時應(yīng)該顧及升降軌引起的像素坐標偏移量的方向變化,其計算公式如下:
(2)
其中:korbit_dir和klook_dir分別對應(yīng)軌道方向(方位向)和視向(距離向),其取值可為1或-1。升軌時取korbit_dir=1,klook_dir為右視時取-1、左視時取1;降軌時取korbit_dir=-1,klook_dir為右視時取1、左視時取-1。x_space、y_space為影像采樣距離向采樣間隔參數(shù)和方位向采樣間隔參數(shù),可在元數(shù)據(jù)文件中通過字段widthspace和heightspace查找。
1.2.2 基于Lee濾波的圖像插值
L1A級別圖像產(chǎn)品為過采樣產(chǎn)品,其方位向采樣間隔通常小于距離向采樣間隔。在進行幾何糾正的后向投影采樣計算中,需要距離向與方位向等間距輸出。本文顧及這一特點,采用Lee濾波直接代替在方位向上進行多視處理后再采用傳統(tǒng)的最近鄰插值或雙線性插值,可在得到的幾何校正結(jié)果圖像中表現(xiàn)出良好的斑噪抑制效果。
通過后向投影得到糾正結(jié)果圖上點P(B,L,H)對應(yīng)的原始影像點P′(r,c),取以P′點為中心的m×n大小的窗口,m為輸出圖方位向采樣間隔除以原圖方位向采樣間隔,n為輸出圖距離向采樣間隔除以原圖距離向采樣間隔。計算得到窗口內(nèi)像素均值、方差后,根據(jù)公式
(3)
(4)
計算濾波結(jié)果。其中x表示無噪聲影像的信號;z表示影像的強度;σv表示乘性噪聲的標準差,假定噪聲均值為1且其標準差σv與信號x無關(guān)。
為驗證本文提出L1A級影像處理流程的正確性,基于高分三號SAR衛(wèi)星影像和Sentinel-1/2影像開展了實驗驗證。選取高分三號和Sentinel-1 2019年5月拍攝的VV極化影像,對比地物的后向散射系數(shù)統(tǒng)計特征以驗證本文方法輻射校正流程的穩(wěn)健性。針對高分三號影像RPC參數(shù)與元數(shù)據(jù)參數(shù)無法自洽,Envi軟件不能正常工作的問題,選取北京地區(qū)和廣東地區(qū)各一景存在元數(shù)據(jù)無法自洽的問題數(shù)據(jù)進行處理,以驗證本文提出幾何糾正處理流程的可靠性。
1)輻射糾正實驗
選取高分三號FSII 10 m分辨率的VV極化影像數(shù)據(jù),以及成像間隔14 d、入射角間隔3.28°的Sentinel-1 IW模式 10 m分辨率的VV極化影像數(shù)據(jù)作為參照進行對比分析,實驗數(shù)據(jù)具體參數(shù)如表1所示。對零值處理、插值處理的結(jié)果進行了對比,實驗對比結(jié)果如圖2所示。
表1 輻射校正實驗數(shù)據(jù)參數(shù)列表Table 1 Imaging parameters in the experiment for radiometric correction
圖2 水域與植被區(qū)在輻射糾正中是否進行零值處理的結(jié)果對比Fig.2 Results comparison of traditional method and proposed method applied to water area and vegetation area
圖2的處理結(jié)果表明,本文提出的輻射校正改正公式剔除了低于等效噪聲系數(shù)的無效后向散射系數(shù)值,在圖2(a2),2(b2)中這些無效值引起的白色斑點噪聲在圖2(a3)、2(b3)中得到了抑制。對圖2(a2)水體區(qū)域,圖2(b2)農(nóng)田區(qū)域中的無效像素進行統(tǒng)計,它們分別占比4.49%、1.23%。由此可見,若對這些L1A產(chǎn)品中出現(xiàn)的量化噪聲不進行處理將會引起輻射強度估計的偏離。
采樣方法的選取影響影像的輻射結(jié)果,針對不同插值后處理采樣方法的對比實驗結(jié)果如圖3所示。
為了對比驗證本文提出的融合Lee濾波方法的重采樣方法對幾何校正結(jié)果影像的處理效果提升,選取影像水體部分的均質(zhì)區(qū)域和建筑部分的高對比度區(qū)域,統(tǒng)計影像像素直方圖,得到結(jié)果如圖3所示。從圖3上看,傳統(tǒng)的最近鄰插值、雙線性插值和本文提出采用的濾波多視處理方法得到的像素統(tǒng)計直方圖上有所差異。最近鄰插值方法結(jié)果的統(tǒng)計直方圖不連續(xù)。雙線性插值方法結(jié)果的地物標準差統(tǒng)計值較本文方法所得結(jié)果較大,在均質(zhì)區(qū)域4.64、高對比度區(qū)域6.56均大于Lee濾波方法結(jié)果的3.07和5.41。計算3種插值方法的等效視數(shù)(equalization number of look, ENL),最近鄰插值方法的ENL為1.75,雙線性插值方法的ENL為1.99,Lee濾波結(jié)果的ENL為2.16。由此可見,本文方法在實現(xiàn)對L1A圖像產(chǎn)品進行方位向插值處理的同時可降低SAR斑點噪聲的影響。
圖3 不同插值方法在影像均質(zhì)區(qū)域和高對比度區(qū)域結(jié)果的像素值統(tǒng)計直方圖Fig.3 Statistical histogram of different interpolation methods in homogeneous region and high contrast region
為了進一步說明本文的方法可以消除影像中無效像素值的影響,選取高分三號影像與Sentine-1影像中典型地物進行像素值統(tǒng)計對比。針對高分三號處理結(jié)果與Sentinel-1處理結(jié)果的輻射強度統(tǒng)計信息對比如表2所示。
表2 高分三號與Sentinel-1影像輻射糾正在典型地物上的統(tǒng)計結(jié)果對比Table 2 Comparison of radiometric statistic results between Gaofen-3 and sentinel-1 images on typical land covers
由表2可見,高分三號所得的同一地物后向散射強度值較Sentinel-1結(jié)果偏小,這與文獻[10]所得結(jié)論一致。但本文處理方法對輻射強度結(jié)果進行了無偏糾正,其有效性表現(xiàn)在以下兩個方面:1)填補了影像中的無效像素,剔除了影像中的統(tǒng)計偏差。本文方法所得水體區(qū)域的中值從-22.33 dB提升到-21.17 dB,農(nóng)田區(qū)域中值從-12.96 dB提升到-11.62 dB;地物統(tǒng)計水體區(qū)域的均值從-22.65 dB提升到-21.97 dB,農(nóng)田區(qū)域的均值從-14.45 dB提升到-14.22 dB;2)零值處理后的影像與Sentinel-1影像的地物統(tǒng)計結(jié)果差異減小,說明零值處理之后的影像消除了無效像素值的影響,修正了影像值的統(tǒng)計偏差。
2)幾何校正實驗
經(jīng)過對高分三號大量L1A級數(shù)據(jù)進行處理發(fā)現(xiàn)部分高分三號數(shù)據(jù)的元數(shù)據(jù)信息存在自洽性問題:1)高分三號影像元數(shù)據(jù)4個角點經(jīng)緯度坐標無法與RPC參數(shù)反算結(jié)果保持一致性。如表3所示,部分高分三號影像元數(shù)據(jù)文件提供的角點坐標經(jīng)過RPC模型獲得的影像像素坐標與真實像素坐標之間存在較大誤差;2)高分三號影像RPC參數(shù)內(nèi)部存在不自洽,導致Envi軟件在利用RPC參數(shù)進行幾何糾正時,部分影像無法完成計算或處理結(jié)果出錯。
從表3可以看到,高分三號不同升降軌方向或左右視方向部分存在元數(shù)據(jù)文件參數(shù)和RPC參數(shù)不自洽的問題。圖上最大誤差達到10 039.0像素,對應(yīng)緯度誤差為0.136°。因此,通過高分三號影像元數(shù)據(jù)文件中提供的角點坐標參數(shù)對影像直接進行RPC正算,無法得到正確的糾正結(jié)果。由此可見,利用RPC反算來獲取影像角點坐標可以保證RPC參數(shù)內(nèi)部的自洽性,對糾正流程的穩(wěn)健性構(gòu)成充分條件,保證處理流程的一致性。
表3 部分高分三號影像元數(shù)據(jù)4個角點坐標與RPC正算的不自洽問題統(tǒng)計表Table 3 The error about self-consistency between corner coordinates from metadata and RPC method in Gaofen-3
針對存在的第2類自洽性問題,即Envi軟件對部分高分三號影像進行幾何校正時無法得到正確結(jié)果,選取2景高分三號影像分別使用Envi軟件和本文方法進行實驗,影像的參數(shù)和Envi軟件錯誤說明如表4所示。Envi對第1景影像無法完成處理,對第2景的處理結(jié)果如圖4所示。本文方法得到的對應(yīng)處理結(jié)果見圖5(a)和5(b)。
表4 Envi軟件無法處理或處理錯誤的高分三號元數(shù)據(jù)參數(shù)說明表Table 4 Parameters of Gaofen-3 image failed in processing with Envi software
圖4 廣東地區(qū)(表4中所列)高分三號影像Envi軟件處理結(jié)果圖Fig.4 Error result of one Gaofen-3 geometric correction via Envi software
圖5 本文方法得到影像的幾何校正結(jié)果圖Fig.5 Geometric correction results by the proposed method (a) Beijing area (b) Guangdong area
從圖5可看出,對于Envi軟件僅依據(jù)RPC參數(shù)無法計算的影像和結(jié)果扭曲的影像,本文提出的方法均可得到正確的幾何校正結(jié)果。這是由于本文方法通過引入元數(shù)據(jù)文件中衛(wèi)星的軌道方向和視向參數(shù)(見式(2))來保證RPC參數(shù)的自洽性,使本文提出的幾何校正流程具有高魯棒性。
Sentinel-2影像95%以上像素定位誤差<11 m[16],為評價本文方法的影像糾正結(jié)果(即圖5)精度,選用Sentinel-2光學影像作為基準數(shù)據(jù)。分別在圖5(a)和5(b)影像中選取20個SAR影像和光學影像上的同名點作為檢查點。同名點之間的坐標誤差作為影像幾何校正結(jié)果的精度誤差。得到的同名點誤差分布如圖6所示。
圖6 高分三號影像(圖5)檢查點誤差分布圖Fig.6 Error distribution of selected check points in Gaofen-3 image
針對圖5(a)所得的5 m分辨率糾正結(jié)果,為避免高程影響,在北京平原地區(qū)選取SAR影像與光學影像中20個地物同名點為檢查點,同名點之間的平面坐標偏差作為幾何校正后的誤差。通過高分三號衛(wèi)星平臺的軌道傾角(10°),將平面坐標誤差旋轉(zhuǎn)到成像坐標系下(距離向、方位向),得到20個同名點在距離向和方位向上的誤差分布如圖6(a)所示。選取的北京平原地區(qū)高分三號FSI模式影像所有同名點的誤差均方根為118.98 m,通過與文獻[14]得到的高分三號影像京津冀地區(qū)FSI模式影像幾何校正結(jié)果誤差109.428~116.290 m相比,可以認為,利用本文提出的幾何校正方法來處理Envi無法處理的圖像,其所得結(jié)果也在合理范圍之內(nèi)。相較于同一地區(qū)的Radarsat-2影像無控制點RPC幾何校正后誤差呈現(xiàn)出均勻分布[17],而圖6(a)中同名點誤差的聚類中心卻位于(121.00,-2.10)處,可推定此高分三號影像在距離向上存在約121 m的系統(tǒng)差。針對圖5(b)廣東地區(qū)SS模式影像所得的25 m分辨率糾正結(jié)果,所選同名點的誤差均方根為862.89 m。分析北京地區(qū)和廣東地區(qū)影像糾正結(jié)果的差異可能由以下原因?qū)е拢?)2景影像區(qū)域內(nèi)的高程差不同。廣東地區(qū)影像大部分屬于丘陵地形,高程起伏較大,在選取的影像同名點分布范圍內(nèi)的影像高程均值為175.495 m,高程差為1 102 m。而北京地區(qū)控制點選在平原區(qū),影像高程為101.934 m,高程差為230 m。高程差異導致了斜視成像的SAR影像幾何校正精度的差異。2)分辨率差異。廣東地區(qū)SS模式影像分辨率為25 m,低于北京地區(qū)FSI影像的5 m分辨率。這導致在粗分辨率影像上選取同名檢查點較為困難,且定位誤差較大。綜合以上分析,本文提出的處理方法可避免高分三號影像數(shù)據(jù)RPC參數(shù)不自洽的問題,且獲得正確的幾何校正結(jié)果。
本文提出的高分三號L1A級影像精處理方法,可以有效地去除原始影像上低于等效噪聲系數(shù)的無效值引起的斑噪影響,提高影像視覺效果。從輻射強度統(tǒng)計量的變化上也說明本文提出的輻射校正流程消除了影像中無效值對影像地物輻射強度的統(tǒng)計偏差。本文提出的幾何校正流程具有高穩(wěn)健性,可以消除高分三號元數(shù)據(jù)信息不自洽問題,處理高分三號影像數(shù)據(jù)中存在的利用Envi軟件無法完成幾何校正或處理結(jié)果扭曲的問題。同時,本文提出的將Lee濾波方法替代常用的最近鄰和雙線性重采樣插值方法,可在實現(xiàn)L1A產(chǎn)品方位向過采樣的同時實現(xiàn)斑噪抑制效果。
綜上所述,本文提出的高分三號L1A級數(shù)據(jù)精處理方法對輻射校正處理結(jié)果和幾何校正處理結(jié)果具有糾正后精度提升和處理流程穩(wěn)健性的特點。但對于如何消除高分三號影像自身存在的絕對輻射值稍低和絕對幾何定位精度稍差的問題,還有待進一步研究。
感謝中國資源衛(wèi)星應(yīng)用中心提供的多模式、多分辨率高分三號L1A級SAR影像,用以驗證本文所提出的方法。