新型SAR傳感器一級地距產品絕對輻射定標方法

杜偉娜， 徐愛功， 宋耀鑫， 孫華生

(1.遼寧工程技術大學,阜新 123000； 2.公安部交通管理科學研究所,無錫 214151)

?

新型SAR傳感器一級地距產品絕對輻射定標方法

杜偉娜1， 徐愛功1， 宋耀鑫2， 孫華生1

(1.遼寧工程技術大學,阜新 123000； 2.公安部交通管理科學研究所,無錫 214151)

針對目前新型SAR傳感器數據產品預處理軟件缺乏的現狀，全面細致地論述了ENVISAT ASAR,Radarsat2,Cosmoskymed,TerraSAR-X和Sentinel1等5種新型SAR傳感器的一級地距產品(level-1 detected ground range product,L-1 DGRP)數據的絕對輻射定標方法以及軟件化過程中所涉及的后向散射系數等參數獲取方法； 并以Sentinel1傳感器的L-1 DGRP數據為例，采用C++語言編程實現了絕對輻射定標過程； 最后將本文方法軟件處理結果和歐空局S1 ToolBox軟件處理結果進行對比，二者得到的后向散射系數值基本相同，證實了本文所介紹方法的正確性。

新型SAR； 絕對輻射定標； 一級地距產品(L-1 DGRP)； 后向散射系數； Sentinel1

0 引言

隨著國外高分辨率合成孔徑雷達(synthetic aperture Radar,SAR)數據對我國各領域的開放以及“龍計劃”項目[1]的深入開展，ENVISAT ASAR,Radarsat1/2,Cosmoskymed,TerraSAR-X,ALOS PALSAR和Sentinel1等高分辨率星載SAR影像數據受到廣大關注。然而，當前能夠支持這些星載SAR數據預處理的商業軟件還比較少； 因此，研究這些SAR傳感器各級產品數據的輻射定標算法并進行軟件實現，對于促進新型高分辨率合成孔徑成像雷達數據的定量應用具有重要意義。與SAR數據的其他各級產品相比，一級地距產品(level-1 detected ground range product,L-1 DGRP)數據在陸地資源環境遙感中的應用更為廣泛，也是衛星地面接收站和數據生產商主要提供的數據產品。L-1 DGRP數據記錄的是雷達微波后向散射信號的振幅，而非地球生物物理參數定量化研究所需的地物后向散射系數[2]。在地表參數反演中，輻射定標處理是不可或缺的，該過程是將傳感器接收的地物后向散射強度信息轉化為地物后向散射系數的唯一途徑。因此，本文根據SAR系統輻射定標的基本原理，介紹了ENVISAT ASAR,Radarsat2,Cosmoskymed,TerraSAR-X和Sentinel1等5種新型星載SAR傳感器的L-1 DGRP數據的絕對輻射定標方法； 并以Sentinel1數據產品的輻射定標為例，采用C++語言編程實現了其絕對輻射定標過程。

1 SAR輻射定標基本原理

星載合成孔徑雷達在探測地物目標時，雷達系統主動向目標發射無線電脈沖，其發射能量部分被目標吸收，其余大部分能量經目標散射后被雷達天線接收，形成回波信號。這一能量傳輸過程可用雷達方程定量表示[3-4]為

(1)

式中：Pr為雷達接收的信號功率；Pt為發射功率；Gr為總接收增益；G(θ,φ)為距離向和方位向的增益；λ為波長；S為地距；σ為目標的雷達截面積。

雷達截面積σ是在假設雷達所接收的回波信號功率密度與實際目標在同等條件下所接收的目標回波信號功率密度相同的理想條件下，定義的實際目標的等效雷達截面積，是表征地物目標散射特性的常用參數。

對于分布目標，雷達截面積σ可表示為后向散射系數的函數，即

σ=σ0Ac,

(2)

式中：σ0為單位面積的雷達截面積(即后向散射系數)，可用來表征目標對電磁波的散射能力；Ac為雷達地面分辨單元的面積。

星載SAR傳感器在發射之前，通常會利用角反射器在地面試驗場進行多次絕對輻射定標實驗； 并且大多數SAR傳感器會將輻射定標模型所需的參數簡化成角度、天線增益、定標常量等[5]，并存儲在數據產品的頭文件中供用戶使用。

2 典型SAR傳感器絕對輻射定標方法

與傳統SAR相比，新一代SAR傳感器普遍具有高空間分辨率、高靈敏輻射強度、多極化模式和多掃描帶寬等特點，表1列出5種典型SAR傳感器的參數[6-11]。

表1 5種新型SAR傳感器的參數

2.1 ENVISAT ASAR數據輻射定標

ENVISAT ASAR[6]是迄今為止在環境監測領域應用最為廣泛的星載合成孔徑雷達，為遙感應用領域的學者提供了大量高質量的數據。ASAR的L-1 DGRP數據包括： ASA_IMP_1P,ASA_IMM_1P,ASA_APP_1P, ASA_APM_1P, ASA_WSM_1P, ASA_IMG_1P和 ASA_APG_1P共7種，通常以“.N1”文件格式分發。圖像中任一像元(i,j)的后向散射系數σ0為

(3)

對后向散射系數求平均，以db為單位表示，其計算公式為

(4)

式中:U和V分別為平均窗口的行和列；N=UV。

以上定標模型涉及的參數均可從“.N1”格式的文件中直接或間接獲取。DN是數據產品的像元值，直接從SAR數據產品文件中的雷達數據記錄中讀取； 標定常數K從ASAR數據文件的元數據中讀取，該參數與極化方式有關，具體讀取方式見表2。

表2 ASAR數據絕對定標常數K

入射角αi,j需要通過計算得到，在注記數據集下的Geolocation Grid ADSRs中記錄了SAR的入射角信息，即影像方位向的行和距離向的列均存在對應的11組入射角，各像元(i,j)的入射角通過插值得到。由于影像在方位向不超過100 km時可認為入射角隨影像方位向行號變化的影響可以忽略不計，因此入射角只隨列號而變化。通常選擇一元二次方程作為插值函數，利用最小二乘法估計，即可計算得到列號j對應的像元入射角αi,j。

2.2 Radarsat-2數據輻射定標

2007年發射的Radarsat2[7]是Radarsat1的后續衛星，它引領了星載SAR的多項技術革命，如第一次出現了全極化模式。Radarsat2的L-1 DGRP數據以Product.xml文件組織，包含TIFF格式的影像數據和XML格式的增益偏移參數文件，其絕對輻射定標過程可以分為獲取雷達截面積影像和計算入射角。

1)獲取雷達截面積影像。Radarsat2影像雷達截面積σ為

σ=(DN2+A)/G ,

(5)

式中：DN為雷達影像像元值；A為增益偏移參數；G為第j像元(即第j列)的增益。其中偏移和增益參數可以在lutSigma.xml文件中獲取。

2)計算像元的入射角。雷達入射角與地球半徑和軌道高度有關，其幾何關系為

(6)

式中：h為軌道高度；R為地球半徑；s為斜距；Sj為圖像中第j列的地距。對于L-1DGRP，在升軌右視和降軌左視情況下(根據Product.xml文件中的passDirection和antennaPointing值判斷)，Sj的計算公式為

Sj=a+jLb+(jL)2c+(jL)3d+(jL)4e+(jL)5f ;

(7)

在降軌右視和升軌左視情況下，Sj的計算公式為

Sj=a+(J-j)Lb+[(J-j)L]2c+[(J-j)L]3d+[(J-j)L]4e+[(J-j)L]5f ,

(8)

式中：L為像元間隔；J為影像總列數；a,b,c,d,e,f分別為地距到斜距的6個轉換參數，可從Product.xml文件中的groundToSlantRangeCoefficients獲取。此外，Sj也可以通過最小入射角和最大入射角插值得到。

綜合以上過程，可得到雷達后向散射系數圖像σ0(以db為單位表示)，即

(9)

2.3 Cosmoskymed數據輻射定標

Cosmoskymed[8]是意大利研發的COSMO-skymed高分辨率雷達衛星星座的代稱，該星座共有4顆衛星，重訪能力高，具備全球范圍觀測能力。Cosmoskymed的L-1 DGRP數據(有DGM和GEC兩種)通常以HDF5格式(“.H5”)分發。其絕對輻射定標過程為

(10)

式中： |imginp(i,j)|2為計算強度影像表達式；αref為參考入射角；F為尺度因子；K為定標常數。以db為單位表示，則

σ0(i,j)db=10lg[σ0(i,j)] ,

(11)

以上參數均可在HDF5文件的元數據中找到，其中αref取自Reference Slant Range；F取自Rescaling Factor；K取自Calibration Constant。

為了減少系統噪聲的影響，通常可以采用一個窗口平均值替代，即

(12)

式中：r,c分別為窗口的行和列；wr,wc分別為窗口的行寬和列寬。

2.4 TerraSAR-X數據輻射定標

TerraSAR-X[9]是一顆成像分辨率和軌道精度都非常高的SAR衛星，在大范圍對地觀測和干涉測量領域具有明顯優勢。TerraSAR-X的L-1 DGRP數據也有2種： 多視地距探測產品(multi look ground range detected,MGD)和地理編碼橢球校正產品(geocoding of terrain correction,GEC)，通常以“.XML”格式組織，包含TIFF影像數據、快視數據和元數據等。TerraSAR-X的L-1 DGRP數據為絕對輻射定標處理提供了定標常量、等效噪聲分布和本地入射角參數，具體過程如式(13)―(16)所示，即

σ0=(ks|DN|2-B)sinθloc,

(13)

(14)

(15)

(16)

式(13)―(16)中：ks為輻射定標常量，對應XML格式元數據文件中calibrationConstant節點下calFactor的參數值；DN為像元的后向散射強度值，影像數據位于IMAGEDATA文件夾；B為強度等效噪聲，反映不同噪聲分布模型對信號的影響，所涉及的參數均可在元數據文件中獲取(由于B的影響很小，在有些情況下該參數可以忽略)；θloc為本地入射角，其計算所需的參數GGIM可以在georef.xml文件中獲取；deg和coffi分別為多項式擬合的階數和系數；τ為快時間；τref為參考快時間；τmin和τmax分別為快時間的最小值和最大值。

2.5 Sentinel1數據輻射定標

歐空局于2014年發射的新一代SAR衛星Sentinel1[11](亦稱“哨兵1號”)是ENVISAT衛星的后繼衛星，繼續擔負著對地環境觀測的任務，其數據可以在歐空局官方網站免費申請獲得。Sentinel1的L-1 DGRP數據代碼為GRD，以“manifest.safe”索引文件組織，包含影像數據(measurement文件夾)、快視數據(preview文件夾)、軌道及標定參數數據(annotation文件夾)以及說明數據(support文件夾)。對于Sentinel1的L-1 DGRP數據的絕對輻射定標，數據產品本身提供了一個定標矢量，將影像的強度值轉化為后向散射系數，具體轉換過程為

(17)

式中：DN為TIFF格式地距影像中對應像元的像素值；Ssigma為定標參數，該參數可通過“annotation”文件夾下的“Calibration”文件夾中的XML元數據文件sigmaNought域得到。元數據中Ssigma是一個查找表(LUT)，為行方向和列方向具有一定間隔的矢量數組； 對于任意像元的Ssigma，需要通過插值得到，插值方法往往選擇雙線性插值方法。

3 實驗與分析

本文以Sentinel1的L-1 DGRP數據為例，采用VC++編程語言和GDAL柵格數據讀寫庫，驗證了絕對輻射定標方法。實驗數據為2015年3月15日獲取的覆蓋環渤海地區Sentinel1的L-1 DGRP，影像大小為19 342行×25 267列，空間分辨率為5 m×20 m，干涉寬幅掃描(IW)成像模式，掃描帶寬為250 km。根據2.5節介紹的方法對Sentinel1的L-1 DGRP數據進行絕對輻射定標。首先讀取影像的像元值，Sentinel1影像采用16 bit記錄強度信息，因此像元值在0～65 535之間； 然后獲取每個像元定標參數值S，由于數據產品提供了30行×633列的定標參數矩陣以及對應的行列坐標矩陣，因此采用雙線性內插方法計算其他行列坐標像元的定標參數值； 最后根據式(17)計算得到每個像元的后向散射系數值，即得到輻射定標處理結果。

圖1(a)為用本文方法的軟件實現的Sentinel1的L-1 DGRP數據輻射定標結果，圖1(b)為用歐空局的S1-ToolBox軟件對相同數據進行輻射定標得到的結果，兩者在視覺上差別甚小。

(a) 本文方法 (b) 歐空局S1-ToolBox軟件

圖1 2種不同方法軟件輻射標定結果

Fig.1 Calibration results by using two kinds of software from different methods

為了進一步分析本文方法絕對輻射定標處理結果的精度，從以上2個輻射定標結果中隨機選取了12組樣本點的后向散射系數值進行對比，誤差分析結果如表3所示。

表3 2種方法標定精度對比

可以看出，2種方法處理結果相對誤差的最大值為0.003，最小值為0。考慮到計算機插值處理精度的問題，以及忽略有效數字的差別，本文方法的軟件處理得到的結果與歐空局S1-ToolBox軟件處理的結果影像是基本相同的，從而證實了本文所介紹方法的正確性。

4 結論

輻射定標是SAR影像地表參數定量反演中的關鍵步驟，本文在明確幾種新型SAR一級地距產品(L-1 DGRP)數據組織結構特點的基礎上，結合輻射定標基本原理，闡述了幾種新型SAR的L-1 DGRP輻射定標方法，結論如下：

1) 針對ENVISAT ASAR,Radarsat2,Cosmoskymed,TerraSAR-X和Sentinel1的L-1 DGRP數據，發展了利用產品的定標參數和入射角插值方法、進而利用輻射定標公式進行后向散射系數計算。

2) 以Sentinel1的L-1 DGRP數據為例，通過對2種方法的軟件處理結果對比研究，驗證了本文所述算法的正確性。

3)下一步的研究將充分考慮地形因子對新型SAR的L-1 DGRP數據輻射定標的影響，并發展更為嚴密的輻射定標模型。

志謝： 感謝歐洲空間局提供ESA S1 ToolBox軟件。

[1] Desnos Y L,Li Z Y,Gao Z H,et al.The dragon programme-status and achievements[C]//Proceedings of the Envisat Symposium 2007.Montreux,Switzerland:ESA,2007.

[2] Laur H,Bally P,Meadows P,et al.ERS SAR Calibration.Derivation of the Backscattering Coefficientσ0in ESA ERS SAR PRI Products[R].ESA/ESRIN ES-TN-RS-PM-HL09.2009:3-6.

[3] Cumming I G,Wong F H.Digital Processing of Synthetic Aperture Radar Data:Algorithms and Implementation[M].Boston,MA:Artech House Remote Sensing Library,2005:56-123.

[4] 袁孝康.星載合成孔徑雷達導論[M].北京:國防工業出版社,2003. Yuan X K.Introduction to Synthetic Aperture Radar[M].Beijing:National Defense Industry Press,2003.

[5] Schwerdt M,Brautigam B,Bachmann M,et al.TerraSAR-X calibration-first results[C]//Proceedings of the 2007 IEEE International Geoscience and Remote Sensing Symposium.Barcelona:IEEE,2007:3932-3935.

[6] Rosich B,Meadows P.Absolute Calibration of ASAR Level 1 Products Generated with PF-ASAR[R].ESA-ESRIN,ENVI-CLVL-EOPG-TN-03-0010,2004:5-8.

[7] Slade B.RADARSAT-2 Product Description[R].RN-SP-52-1238,2009:2-25.

[8] Agenzia Spaziale Italiana.COSMO-SkyMed SAR Products Handbook[M].Agenzia:ASI,2007:399-405.

[9] Fritz T,Werninghaus R.TerraSAR-X Ground Segment Level 1b Product Format Specification[R].Berlin:Clustert Applied Remote Sensing(CAF),German Aerospace Center(DLR),2007:156-170.

[10]Rosenqvist A,Shimada M,Suzuki S,et al.Operational performance of the ALOS global systematic acquisition strategy and observation plans for ALOS-2 PALSAR-2[J].Remote Sensing of Environment,2014,155:3-12.

[11]Geudtner D,Torres R,Snoeij P,et al.Sentinel1 system capabilities and applications[C]//Proceedings of the 2014 IEEE Geoscience and Remote Sensing Symposium.Quebec City,QC:IEEE,2014:1457-1460.

(責任編輯： 邢宇)

Absolute radiometric calibration of level-1 detected ground range products of new SAR sensors

DU Weina1， XU Aigong1， SONG Yaoxin2， SUN Huasheng1

(1.LiaoningTechnicalUniversity，Fuxin123000，China； 2.TrafficManagementResearchInstituteofMinistryofPublicSecurity，Wuxi214151,China)

For the current situation of the lack of the new SAR sensor data preprocessing software, this paper introduced in detail the methods of absolute radiometric calibration and the parameter acquisition for several new SAR sensor level-1 detected products, such as ENVISAT ASAR，Radarsat2，Cosmoskymed，TerraSAR-X and Sentinel1. In addition, the absolute radiometric calibration process was achieved by programming with the level-1 detected ground range products(L-1 DGRP) data of Sentinel1 sensor, and C++ programming language was used to achieve the absolute radiation of the calibration process. At last, the radiometric calibration results produced by the method developed in this paper and implemented in the authors’ software were compared with those by ESA S1 ToolBox, the freely distributed SAR data processing tool by European Space Agency, and it is shown that the two numerical back scattering systems are basically the same. The radiometric calibration method developed in this paper is proved to be correct by the program implementation.

new SAR sensors； absolute radiometric calibration； level-1 detected ground range products(L-1 DGRP)； back scattering coefficient； Sentinel1

10.6046/gtzyyg.2016.04.05

杜偉娜,徐愛功,宋耀鑫,等.新型SAR傳感器一級地距產品絕對輻射定標方法[J].國土資源遙感,2016,28(4)：30-34.(Du W N,Xu A G,Song Y X,et al.Absolute radiometric calibration of level-1 detected ground range products of new SAR sensors[J].Remote Sensing for Land and Resources,2016,28(4)：30-34.)

2015-05-13；

2015-07-08

國家自然科學基金項目“遙感數據的空間分辨率和波段數對土地覆蓋制圖的影響研究”(編號： 41201454)資助。

TP 751.1

A

1001-070X(2016)04-0030-05

杜偉娜(1988-),女,碩士研究生,主要研究方向為SAR圖像處理與應用。Email： RSwendu@163.com。