先进的对地观测卫星二号(ALOS-2)于2014年5月24日由日本宇航局发射,是继先进的对地观测卫星(ALOS)之后的新一代L波段雷达卫星。作为上一代ALOS卫星的后续星,ALOS-2以更好的性能延续了ALOS 卫星在区域制图、灾害监测和资源调查等方面的优势,可继续发挥重大的作用^[1]。特别是为了响应更广泛的社会需求,其搭载的PALSAR-2合成孔径雷达系统获取的影像细节更清晰、覆盖范围更广,成像模式也更为丰富。Arikawa等^[2]指出,ALOS-2在通道数量、量测精度和导航精度方面的性能指标较ALOS有很大的提升,其数据产品的质量评价对进一步推广应用有重要的意义^[3]。但是,由于斜距成像的成像原理,影像中存在严重的几何变形,对PALSAR-2影像进行正射校正及精度评价显得尤为重要,且目前针对同一类型的2代同波段传感器的校正精度差异比较等相关研究成果尚少。

就校正方法而言,基于成像几何模型结合数字高程模型(digital elevation model,DEM)和距离-多普勒(range doppler,RD)定位模型对合成孔径雷达(synthetic aperture Radar,SAR)影像进行正射校正方法是应用较广泛的一类方法^[4,5,6]。该方法无需地面控制点,适用于对不易人工获取地面控制点的区域进行正射校正^[7]。但是由于用于计算的卫星轨道参数^[8,9]、DEM高程^[10]及成像参数存在误差,经校正后的图像仍然存在误差^[11]。校正过程中关键部分是如何计算准确的方位向与距离向偏移,为此一些学者曾提出了基于线性逼近和迭代来解算RD模型的方法^[12],但这类方法计算量大,且精度不是很高。Shimada给出了一种简单而且精确的正射校正方法^[13],但是要不断判断轨道精度,在轨道精度不满足要求的时候,需要通过模拟影像与真实影像的配准进行反复计算。

本文提出一种基于轨道参数修正和RD模型简化解算的PALSAR-2影像校正方法,利用模拟SAR与真实SAR影像的配准,修正轨道参数,再利用修正后的轨道参数与RD模型,完成校正工作。为了对同一类型的2代同波段传感器的校正精度差异进行比较,评估新一代L波段传感器的性能,将该方法同时应用于PALSAR-2与PALSAR影像,并对2种影像校正结果精度进行比较; 并与没有经过轨道参数修正的PALSAR-2影像校正结果进行了对比,验证方法的有效性。

PALSAR-2成像几何模型为侧视模式,根据地面点与天线中心距离大小按顺序记录地面点位置。因此,地面点高度变化引起图像几何变形非常明显。图1展示了成像过程中的地面点偏移情况。

图1中,SAR为卫星位置; P为地面点; h为高程; θ_inci为正常入射角; R_S为卫星与P之间的距离; P₀为P点在椭球面上的投影点; R'_S为P₀与卫星之间的距离; 目标P用多普勒频率f_D1表示,多普勒中心频率f_D和f_D1有一微小差值Δf_D(因为高程非零),由于非零像素高,P点沿距离向偏移到椭球面上的P₁点; Δr_g为地距偏移量; Δr为斜距偏移量; 由于多普勒频移Δf_D的影响,使得P₁点沿方位向偏移到P₂点; Δx为偏移量。像点位移的存在,使得图像存在叠掩、阴影、透视收缩和近距离压缩等变形,这也使得正射校正尤为重要。

基于DEM进行PALSAR-2影像模拟,对真实PALSAR-2影像进行校正的主要步骤如下:

1)读取SAR影像头文件,获取卫星运动参考坐标系统,并根据SAR头文件中提供的N个时刻的卫星位置矢量和速度矢量对卫星轨道进行拟合,生成轨道参数。

2)读取DEM影像并进行坐标系转换,转换为SAR头文件中读取卫星运动的参考坐标系。

3)根据SAR影像成像参数、拟合后的卫星轨道参数以及坐标系转换后的DEM,建立RD定位模型,进行影像模拟,生成模拟SAR影像。

4)模拟SAR影像与真实SAR影像进行特征点检测与匹配。

5)根据匹配得到的控制点对卫星轨道参数进行修正。

6)根据修正后的卫星轨道数据,重新建立RD模型。

7)求解RD模型,生成正射校正影像。

8)精度验证与分析。

具体流程如图2所示。

本文主要对卫星轨道描述、RD定位模型与简化计算以及轨道参数修正进行详细介绍,其他步骤采用常规方法进行。

SAR数据头文件中都会提供N个时刻的卫星位置矢量和速度矢量的记录,可以采用二阶及以上的多项式方法拟合出卫星随时间变化的轨道曲线方程,随着阶数增加,计算复杂度也会增加,且二阶和三阶多项式定位差异仅在cm级^[14],故本文使用二阶多项式对轨道进行拟合,天线相位中心位置关于时刻t的二阶多项式模型为

$\begin{bmatrix}X_{S} \\Y_{S} \\Z_{S}\end{bmatrix}=\begin{bmatrix}a_{0} & a_{1} & a_{2} \\b_{0} & b_{1} & b_{2} \\c_{0} & c_{1} & c_{2} \\\end{bmatrix}\begin{bmatrix}1 \\t \\t^{2} \\\end{bmatrix}$ ,(1)

式中: (X_S,Y_S,Z_S) 为卫星天线相位中心的三维位置; (a_k,b_k,c_k)为多项式系数,k=0,1,2。

式(1)对时间t求导,可得任一像点成像瞬间卫星天线相位中心的瞬时速度(V_X,V_Y,V_Z),即

$\begin{bmatrix}V_{X} \\V_{Y} \\V_{Z}\end{bmatrix}=\begin{bmatrix}a_{1} & a_{2} \\b_{1} & b_{2} \\c_{1} & c_{2} \\\end{bmatrix}\begin{bmatrix}1 \\t \\\end{bmatrix}$ ,(2)

利用头文件中提供的位置矢量和速度矢量,依据最小二乘原则,求取初始多项式系数,即本文所选取的轨道参数。

采用RD定位模型进行定位,实验中卫星轨道记录在地心旋转坐标系中,地面点相对于坐标轴是静止的。运动速度为0,建立RD模型为

$\begin{cases}R_{S}=|r_{P}-r_{S}| \\f_{D1}=\frac{2f_{0}}{cR_{S}}V_{S}(r_{P}-r_{S})\end{cases}$, (3)

式中: R_S可以通过雷达参数求取; r_P为地面点P的位置矢量; r_S和v_S分别为卫星位置矢量和参数矢量,均可以通过雷达参数求取; f₀为载波频率,均可以从头文件中读取; c为光速。

为了简化计算,把式(3)第二个公式中r_P用r_P-r_P1+r_P1代替,可得

f_D1= 2f0cRSv_S(r_P1-r_S)+ 2f0cRSv_S(r_P-r_P1)=f_D+Δf_D , (4)

Δf_D为相对地面点P相对卫星运动引起的多普勒频移,频移导致的方位向偏移量为

Δx=- ΔfDfDDv_g , (5)

式中: f_DD为多普勒调频率; v_g为卫星在地面上沿方位向的速度。因此,地面点P沿距离向偏移到P₁,再沿方位向偏移到P₂(图1)时,地面点高程引起距离向偏移,偏移量为

Δr_g≅ htanθinci。 (6)

卫星运动速度随时间变化导致多普勒频移也是随时间变化的,其变化规律可以利用三次多项式拟合。通常情况下,SAR头文件中都提供了多普勒中心频率变化的拟合参数,一种典型的计算方法为

f_D=d₀+d₁(t-t₀)+d₂(t-t0)2+d₃(t-t0)3, (7)

t₀=2R₀/c , (8)

式中: t₀和R₀分别为初始时间和初始斜距; d₀,d₁,d₂和d₃分别为头文件中给出的方程参数; t为当前时间,可以通过雷达参数转换为采样窗口方位向数据地址。

P₁点坐标可以通过下式求得,即

$\begin{cases}f_{D}=\frac{2f_{0}}{cR_{S}}V_{S}(r_{P1}-r_{S}) \\\frac{x^2_{P1}}{R^2_a}+\frac{y^2_{P1}}{R^2_a}+\frac{z^2_{P1}}{R^2_b}=1 \\R_{b}=R_{a}\sqrt{1-e^2sin^{2}\varphi} \\e^2=\frac{R^{2}_{a}-R^2_{b}}{R^2_{a}}\end{cases}$, (9)

式中: x_P1,y_P1和z_P1分别为矢量r_P1的3个元素; R_a为赤道半径; R_b为极半径; e为椭球第一偏心率; φ为大地纬度。

对于SAR影像上任意一个像元(i,j),满足如下关系

$\begin{cases}R_{S}=R_{0}+j\delta_{r} \\T=iM/PRF+T_{0}\end{cases}$, (10)

式中: δ_r为距离向斜距像元大小; T为像元(i,j)的成像时间; PRF为SAR的脉冲重复频率; M为方位向视数; T₀为初始成像时间。

根据修正后的轨道参数重新拟合卫星轨道,对于PALSAR-2影像任意一像元(i,j),根据式(10)以及R_S,R'_S与h之间的关系(图1),建立高程h关于(i,j)的函数,遍历PALSAR-2影像,求解出每一像元所对应的高程h,生成高程影像。将h代入式(6),即可求解出PALSAR-2每一像元所对应的距离向偏移量Δr_g,生成距离向偏移影像。联合式(4),(7)和(9),求解出地面点P相对卫星运动引起的多普勒频移Δf_D,代入式(5)即可求解出PALSAR-2每一像元所对应的方位向偏移量Δx,生成方位向偏移影像。根据这些影像完成校正工作。

构建函数如下:

F1=RS-rP-rSF2=fD1-2f0cRSvS(rP-rS)。 (11)

上述函数均和轨道参数相关,使用泰勒级数展开为线性形式,即

$\begin{cases}F^0_1=\sum_{k=0}^{2}\frac{\partial{F_1}}{\partial{a_k}}da_{k}+\sum_{k=0}^{2}\frac{\partial{F_1}}{\partial{b_k}}db_{k}+\sum_{k=0}^{2}\frac{\partial{F_1}}{\partial{c_k}}dc_{k} \\F^0_2=\sum_{k=0}^{2}\frac{\partial{F_2}}{\partial{a_k}}da_{k}+\sum_{k=0}^{2}\frac{\partial{F_2}}{\partial{b_k}}db_{k}+\sum_{k=0}^{2}\frac{\partial{F_2}}{\partial{c_k}}dc_{k}\end{cases}$, (12)

式中: ∂F1∂ak, ∂F1∂bk和 ∂F1∂ck为F₁对轨道参数的偏导; ∂F2∂ak, ∂F2∂bk和 ∂F2∂ck为F₂对轨道参数的偏导; F10和 F20分别为F₁和F₂的初始值。

对模拟SAR影像与真实SAR影像进行匹配,将匹配得到的控制点坐标代入,迭代求解出修正量 ∑k=02da_k, ∑k=02db_k和 ∑k=02dc_k,对轨道参数进行修正,当迭代次数大于最大迭代次数或单位权方差变化量小于阈值时终止迭代。

选取日本富士山地区为主要研究区域,地理位置为E138.30°~138.99°,N35.10°~35.70°,最高海拔为3 759 m,最低海拔为0 m,地形复杂多变。实验数据包括研究区的L波段ALOS PALSAR影像、ALOS-2 PALSAR-2影像,详细信息见表1,以及相应区域的SRTM 30 m空间分辨率的 DEM(来源于 USGS 网站)。

2.2.1 PALSAR-2影像预处理

采用的PALSAR-2影像为CEOS格式数据,需预先用SARscape将原始CEOS格式数据转为单视复数据(single look complex,SLC),并进行多视处理,生成多视处理之后的幅度图。此处多视处理的目的有2个: ①SLC图像数据中包含大量的斑点噪声,为了降噪,需进行多视处理,降低了空间分辨率,但同时提高辐射分辨率; ②SLC图像的距离向和方位向空间分辨率往往不一致,通过多视处理可以把两者的空间分辨率调整到相近的值。

PALSAR-2 SLC影像空间分辨率为5.54 m(距离向(地距))×2.84 m(方位向),经过距离向4视,方位向8视处理后,得到的多视图像空间分辨率为22.16 m(距离向(地距))×22.72 m(方位向),如图3所示。

2.2.2 DEM预处理

采用的DEM是SRTM 30 m空间分辨率DEM数据,DEM经纬度范围包含研究区,需预先确定PALSAR影像所对应研究区的经纬度范围,并对DEM进行裁剪,裁剪出研究区所对应的DEM影像。

由于精度需求,对DEM进行过采样。DEM过采样率f确定方法为

f≥ δrδrg2ΔSdemΔSout, (13)

式中: δ_r为斜距空间分辨率; δ_rg为地距空间分辨率; ΔS_out为模拟SAR影像斜距空间分辨率,可以和δ_r相同,也可以不同; ΔS_dem为DEM的空间分辨率。

研究区影像经多视处理后可得到,δ_rg=22.16 m,δ_r=11.44 m,ΔS_dem=30 m,ΔS_out=11.44 m。根据式(13),可以计算出过采样率f≥1.91。因此本文对DEM进行2倍过采样,处理后的DEM影像如图4所示。

根据PALSAR-2头文件中位置参数,建立卫星轨道拟合公式,生成初始轨道参数,采用常规的方法进行影像模拟,生成PALSAR-2模拟影像(图5)。

对比图5和图3可以看出,基于DEM的PALSAR-2模拟影像可以很好地反映出山地地形中的地形起伏特征,纹理信息很明显。

采用基于特征点的图像配准方法,对真实PALSAR-2影像和模拟PALSAR-2影像进行SURF特征点自动检测与匹配,将匹配的控制点代入式(12)对轨道参数进行修正,再根据修正后的轨道参数及RD模型简化计算的方法进行正射校正,得到正射校正后的PALSAR-2影像,如图6所示。

对比图6和图3可以看出,正射校正后,透视收缩等变形得到很好改善,表明该方法较为有效,具有可行性。为了对同一类型的2代同波段传感器的校正精度差异进行比较,评估新一代L波段传感器的性能,同时对PALSAR影像也进行相同实验,对实验用PALSAR影像进行校正,影像如图7所示。

为了评价该方法的校正精度,对PALSAR-2影像进行了没有经过轨道修正的正射校正,并对2种方法校正结果进行比较。PALSAR-2未经轨道修正校正的影像如图8所示。

采用正射校正后的空间分辨率为4.78 m的光学卫星影像图作为参考基准,在研究区均匀选取15个明显地物点作为检验点,如线路交叉点或拐点,在校正的PALSAR-2校正影像上和卫星影像图上都能准确识别。分别读取每一个检验点在PALSAR-2校正影像和卫星影像上的坐标,并计算两者的坐标差△x和△y(见表2)。同时,对PALSAR-2未经轨道修正校正影像也进行精度评价,如表2所示。根据表2计算均方根误差RMSE,即

RMSE= ∑i=1n(Pri-Pci)2n, (14)

式中: n为检验点的个数; P_ri为真实坐标; P_ci为校正坐标。通过计算可知,PALSAR-2校正的RMSE为8.32 m(0.55个像元),同样,PALSAR的RMSE

为10.91 m(0.73个像元),PALSAR-2未经轨道修正方法校正的RMSE为14.03 m(0.94个像元)。结果表明,新一代L波段传感器影像的校正精度较高,这也进一步证实了新一代L波段传感器有更强的性能指标; 而且相比之下,本文给出的校正方法较未经轨道修正方法校正精度有明显提高。

基于轨道参数修正和RD模型简化解算的PALSAR-2影像校正方法通过利用模拟SAR与真实SAR影像的配准,修正轨道参数,再利用修正后的轨道参数与RD定位模型,完成校正工作。将该方法同时应用于PALSAR-2与PALSAR影像,并对2种影像校正结果精度进行比较,结果表明新一代L波段传感器影像的校正精度较高,这也进一步证实了新一代L波段传感器有更强的性能指标。同时,将该方法与基于DEM和RD定位模型的PALSAR-2影像直接地理定位结果进行了比较,结果表明本文给出的校正方法可操作性强,相比于没有经过轨道参数修正的校正结果精度有明显提高。

由于PALSAR-2的空间分辨率远高于DEM的空间分辨率,这在一定程度上限制了PALSAR-2的校正精度。将来,更高空间分辨率的DEM可以在很大程度上提高PALSAR-2的校正精度。同时,由于只采用单景影像,叠掩和阴影等没有得到很好改善,在今后的研究中,可以采用多景影像来有效改善叠掩和阴影等变形。