高分三号卫星是中国首颗C波段高分辨率合成孔径雷达(synthetic aperture radar, SAR)卫星。它具有12种不同的工作模式, 在聚束模式下可达到1 m分辨率, 为海洋监测^[1-2]、建筑提取^[3-4]、建筑物损毁评估^[5]、船舶检测^[6]、海面风场反演^[7]等不同应用提供了数据支持。高分三号发布的L1A级产品是过采样复数据产品, 因其含有相位、强度、极化信息, 而被广泛订购使用。L1A数据集中包含*.meta.xml、*.rpc、*.incidence.xml等辅助数据文件。这些辅助数据文件集的提供, 为用户进行自定义的更高级别图像产品形式的生产提供了方便。

SAR影像产品处理方式和等级不同, 其像素值可由多种物理量表示：DN(digital number)值、幅度值、强度值、后向散射系数(sigmma nought)值^[8]。DN值为无符号整型数值, 常用于SAR影像复数据的短有效位存储; 幅度值与强度值对应SAR影像像素复数据的模及模平方; 后向散射系数值对应SAR影像强度值经辐射纠正后的数值表达。DN值、幅度值、强度值、后向散射系数值四者之间为指数或对数的变换关系, 因此它们之间可相互转换并为正相关关系。

经辐射纠正后, SAR图像像素值常用后向散射系数表示, 其大小对应被照射地物在该波段上对入射雷达波反射能力的强弱。基于后向散射系数图, 同一传感器在不同时刻获取的图像, 以及同一时刻成像场景中不同位置处的地物后向散射特性均可做定量化对比, 以便于基于后向散射系数实现生物量参数反演等。这些定量对比与参数反演的应用就要求SAR图像辐射纠正要有较高的精度与准确度。

针对高分三号卫星影像单极化图像的绝对辐射强度纠正、全极化数据的通道不平衡和串扰均有学者进行了实验研究。Chang等^[9]指出通常情况下高分三号卫星通道强度不平衡优于0.5 dB, 通道相位差优于10°, 串扰优于-35 dB。对少数残差较大的图像, 通过进一步极化纠正, 可以将极化通道不平衡降低到0.26 dB, 通道相位差降低到±0.2°, 串扰精度优于-42 dB。Chen等^[10]利用三面角角反射器在河北地区做了高分三号的同步散射强度测量和图像质量评价, 依据实测结果发现高分三号卫星的等效噪声系数为-25 dB, 稍微优于Radarsat-2影像的-24.6 dB; 高分三号卫星绝对辐射定标值相较于Radarsat-2要低3~5 dB。高分三号不同模式产品对应的标称辐射纠正参数均包含在L1A数据集中, 普通用户可依据这些参数自行进行后向散射系数图生成, 但采用何种方式实现多视处理、克服L1A级产品量化过程中的空值与零值干扰, 获取纠正结果图的高精度输出仍值得研究。

星载SAR几何纠正是消除卫星成像时由于太阳同步轨道面倾斜、地球自转以及局部地形等因素造成的影像几何畸变的过程^[11]。已经发展了多种星载SAR影像几何纠正方法, 大致可分为严密的物理模型方法以及通用的几何模型。张过等^[12-13]通过选取地面控制点(ground control point, GCP)建立了通用的有理多项式(rational polynomial coefficient, RPC)模型, 并验证了RPC模型具有对传统SAR严密物理模型即距离多普勒(range-Doppler, RD)模型的可替代性。中国高分三号卫星数据在斜距方向与方位时间的标定精度有限。吕冠男等^[14]选取京津冀地区高分三号影像数据, 探究影像的几何定标后影像定位精度。实验证明, 高分三号影像无控制点几何校正平面误差范围为109.428~116.290 m。在高分三号SAR影像L1A数据集中的*.rpc文件中提供了可进行几何纠正的RPC模型参数。根据高分三号卫星提供的用户手册, 用户可以只通过影像对应的*.rpc文件, 利用最新版的Envi5.5软件对影像进行几何校正。然而, Envi5.5软件在对部分高分三号影像进行几何纠正时, 存在无法完成几何纠正或输出结果出错的问题。

为实现从L1A级到更高等级产品处理结果的高可靠和高精度, 提出一套完整的处理流程用以克服L1A级图像产品过采样问题、辐射量化过程中出现空值与零值带来辐射纠正结果偏差的问题、以及RPC参数文件在Envi软件中出现失效导致几何纠正无法进行的问题等。并通过与同波段的Sentinel-1 SAR影像、Sentinel-2光学影像进行对比证明了本文方法的有效性。

1 高分三号L1A级产品精处理方法

高分三号卫星L1A级数据产品中的栅格图像是投影方向为斜距方向上的单视复数据(single look complex, SLC), 其像素值为量化后的DN值。在遥感综合应用中, 光学影像、DEM数据等均为投影面为地面的栅格图。因此, 在后续处理中也须对L1A级产品进行辐射强度与几何投影的标准化。鉴于此, 提出对L1A级产品进行辐射与几何纠正的一体化处理流程如图 1所示。

在该处理流程的辐射纠正中, 依据等效噪声系数对SAR图像中各像素值进行检核, 以剔除SAR图像量化中空值与零值对定标结果的影响。在几何校正中提出RPC反算方法, 以应对高分三号元数据文件提供的角点坐标与RPC参数无法自洽的问题。在后向投影采样输出中将传统的Lee滤波融合到采样过程中, 以替代传统的多视处理或双线性插值处理, 抑制SAR图像斑点噪声的影响。其各要点方法分解描述如下。

1.1 高分三号SAR影像辐射纠正

对高分三号SAR影像L1A级别产品进行辐射纠正, 也就是每个像素值对应后向散射系数σ_dB⁰的计算可通过下式实现：

$\begin{cases}\sigma_{\mathrm{dB}}^0=10 \log _{10}\left(P^{I *}\left(\frac{\text { QualifyValue }}{32767}\right)^2\right)-K_{\mathrm{dB}}, & \sigma_{\mathrm{dB}}^0>N E_{\sigma_{\mathrm{dB}}^0}, \\ \sigma_{\mathrm{dB}}^0=N E_{\sigma_{\mathrm{dB}}^0}, & \sigma_{\mathrm{dB}}^0 \leqslant N E_{\sigma_{\mathrm{dB}}^0}\end{cases}$

(1)

其中: P^I=I²+Q², I为L1A级产品像素值的实部, Q为1A级产品像素值的虚部; QualifyValue为该景图像量化前的最大值, K_dB为定标常数, 这2个常数值可在元数据文件(*.meta.xml)中分别通过QualifyValue和CalibrationConst字段查询获取。

在式(1)中P^I值须大于零, 否则无数学解。在一幅SAR图像中, σ_dB⁰的最小值对应于SAR传感器的等效噪声系数(NE_σdB⁰), 它是SAR传感器辐射灵敏度的衡量指标^[15]。因此, 在辐射纠正过程中, 当读取到P^I为零或空(Null)值时, 辐射纠正的结果值应当以NE_σdB⁰值进行代替, 且保证辐射纠正结果大于等于σ_dB⁰。若以其他非物理意义的最小值替代NE_σ⁰dB, 将会引起对匀质区域(如水域、森林、草地等)后向散射特性统计量的偏差。

1.2 高分三号SAR影像几何纠正 1.2.1 基于RPC模型参数的几何纠正

高分三号SAR卫星L1A级产品在元数据文件集中提供了以(*.rpc)命名的RPC参数文件。直接利用该参数文件, 可以实现地面上一点P(B, L, H)到图像上对应像点坐标P′(r, c)的转换。其中(B, L, H)代表该点的纬度、经度、高程; (r, c)代表图像上的行列坐标。RPC参数模型及其中参数的含义见文献[12]。

为方便利用RPC参数实现几何纠正过程, 常采用如下所述的后向投影过程予以实现：1)求出待输出几何纠正结果图幅框4个角点的地理坐标, 2)采用RPC模型建立从纠正结果图P(B, L, H)到原图P′(r, c)的后向投影关系, 3)设定输出结果中各像素间距(输出步长)后, 在原图中进行采样实现几何纠正结果的输出。

在上述纠正过程中的第一步尤为重要, 图幅框4个角点的像素坐标可通过RPC模型反算原图的4个角点经纬度坐标P₁(B_TL, L_TL), P₂(B_TR, L_TR), P₃(B_BL, L_BL), P₄(B_BR, L_BR), 而后通过求取最大外接矩得到校正后影像范围。4个角点的经纬度坐标可通过图 1中RPC几何校正模块所示的反算迭代流程予以实现, 其过程概括如下：

1) 设定计算初始值为影像中心经纬度坐标(B₀, L₀);

2) 将起算点(B₀, L₀)赋值到迭代点(B′_i, L′_i);

3) RPC正算计算迭代点(B′_i, L′_i)对应的影像像素坐标(r, c), 与待求点像素坐标P_corner(r, c)作差得到待求点与迭代点像素坐标差(Δr, Δc);

4) 将像素坐标差(Δr, Δc)转换成为投影平面坐标差(Δx, Δy);

5) 计算迭代点(B′_i, L′_i)的投影平面坐标(x′_i, y′_i), 迭代点投影平面坐标(x′_i, y′_i)加投影平面坐标差(Δx, Δy)得到迭代后投影平面坐标(x′_i+1, y′_i+1)=(x′_i, y′_i)+(Δx, Δy);

6) 高斯反解计算迭代后迭代点经纬度坐标(B′_i+1, L′_i+1), 计算迭代前后影像经纬度坐标差(ΔB, ΔL);

7) 若坐标差(ΔB, ΔL)小于阈值, 即ΔB < 4.8×10^-9且ΔL < 4.8×10^-9, 或迭代次数达到上限则转入8), 否则迭代点赋值为(B′_i+1, L′_i+1)转入3);

8) 输出角点经纬度坐标P_corner(B, L)为最终解。

在上述RPC模型反算过程中有以下两个关键点说明如下：1)解算起始点参数的设置对反算结果的计算效率影响较大, 起算点设置不合理会导致算法不收敛。RPC参数文件中提供的经纬度偏移参数(Lat_Offset, Lon_Offset), 其物理含义为景中心经纬度, 将其作为初始值可以提高迭代速度, 且保证算法收敛。2)L1A级别产品的栅格影像存储顺序以SAR传感器接收到回波的先后次序进行存储, 因此在实现从像素坐标差到大地坐标差的转换时应该顾及升降轨引起的像素坐标偏移量的方向变化, 其计算公式如下：

$\left\{\begin{array}{l}\Delta x=\Delta r \times k_{\text {orbit_dir }} \times k_{\text {look_dir }} \times x_{-} \text {space , } \\ \Delta y=\Delta c \times k_{\text {orbit_dir }} \times k_{\text {look_dir }} \times y_{-} \text {space. }\end{array}\right.$

(2)

其中：k_{orbit_dir}和k_{look_dir}分别对应轨道方向(方位向)和视向(距离向), 其取值可为1或-1。升轨时取k_{orbit_dir}=1, k_{look_dir}为右视时取-1、左视时取1;降轨时取k_{orbit_dir}=-1, k_{look_dir}为右视时取1、左视时取-1。x_space、y_space为影像采样距离向采样间隔参数和方位向采样间隔参数, 可在元数据文件中通过字段widthspace和heightspace查找。

1.2.2 基于Lee滤波的图像插值

L1A级别图像产品为过采样产品, 其方位向采样间隔通常小于距离向采样间隔。在进行几何纠正的后向投影采样计算中, 需要距离向与方位向等间距输出。本文顾及这一特点, 采用Lee滤波直接代替在方位向上进行多视处理后再采用传统的最近邻插值或双线性插值, 可在得到的几何校正结果图像中表现出良好的斑噪抑制效果。

通过后向投影得到纠正结果图上点P(B, L, H)对应的原始影像点P′(r, c), 取以P′点为中心的m×n大小的窗口, m为输出图方位向采样间隔除以原图方位向采样间隔, n为输出图距离向采样间隔除以原图距离向采样间隔。计算得到窗口内像素均值、方差后, 根据公式

$\hat{x}=\bar{x}+k(z-\bar{x})$,

(3)

$k=\frac{\operatorname{var}(x)}{\bar{x}^2 \sigma_v^2+\operatorname{var}(x)}$,

(4)

计算滤波结果。其中x表示无噪声影像的信号; z表示影像的强度; σ_v表示乘性噪声的标准差, 假定噪声均值为1且其标准差σ_v与信号x无关。

2 实验结果及分析

为验证本文提出L1A级影像处理流程的正确性, 基于高分三号SAR卫星影像和Sentinel-1/2影像开展了实验验证。选取高分三号和Sentinel-1 2019年5月拍摄的VV极化影像, 对比地物的后向散射系数统计特征以验证本文方法辐射校正流程的稳健性。针对高分三号影像RPC参数与元数据参数无法自洽, Envi软件不能正常工作的问题, 选取北京地区和广东地区各一景存在元数据无法自洽的问题数据进行处理, 以验证本文提出几何纠正处理流程的可靠性。

1) 辐射纠正实验

选取高分三号FSII 10 m分辨率的VV极化影像数据, 以及成像间隔14 d、入射角间隔3.28°的Sentinel-1 IW模式10 m分辨率的VV极化影像数据作为参照进行对比分析, 实验数据具体参数如表 1所示。对零值处理、插值处理的结果进行了对比, 实验对比结果如图 2所示。

图 2的处理结果表明, 本文提出的辐射校正改正公式剔除了低于等效噪声系数的无效后向散射系数值, 在图 2(a2), 2(b2)中这些无效值引起的白色斑点噪声在图 2(a3)、 2(b3)中得到了抑制。对图 2(a2)水体区域, 图 2(b2)农田区域中的无效像素进行统计, 它们分别占比4.49 %、1.23 %。由此可见, 若对这些L1A产品中出现的量化噪声不进行处理将会引起辐射强度估计的偏离。

采样方法的选取影响影像的辐射结果, 针对不同插值后处理采样方法的对比实验结果如图 3所示。

为了对比验证本文提出的融合Lee滤波方法的重采样方法对几何校正结果影像的处理效果提升, 选取影像水体部分的均质区域和建筑部分的高对比度区域, 统计影像像素直方图, 得到结果如图 3所示。从图 3上看, 传统的最近邻插值、双线性插值和本文提出采用的滤波多视处理方法得到的像素统计直方图上有所差异。最近邻插值方法结果的统计直方图不连续。双线性插值方法结果的地物标准差统计值较本文方法所得结果较大, 在均质区域4.64、高对比度区域6.56均大于Lee滤波方法结果的3.07和5.41。计算3种插值方法的等效视数(equalization number of look, ENL), 最近邻插值方法的ENL为1.75, 双线性插值方法的ENL为1.99, Lee滤波结果的ENL为2.16。由此可见, 本文方法在实现对L1A图像产品进行方位向插值处理的同时可降低SAR斑点噪声的影响。

为了进一步说明本文的方法可以消除影像中无效像素值的影响, 选取高分三号影像与Sentine-1影像中典型地物进行像素值统计对比。针对高分三号处理结果与Sentinel-1处理结果的辐射强度统计信息对比如表 2所示。

由表 2可见, 高分三号所得的同一地物后向散射强度值较Sentinel-1结果偏小, 这与文献[10]所得结论一致。但本文处理方法对辐射强度结果进行了无偏纠正, 其有效性表现在以下两个方面：1)填补了影像中的无效像素, 剔除了影像中的统计偏差。本文方法所得水体区域的中值从-22.33 dB提升到-21.17 dB, 农田区域中值从-12.96 dB提升到-11.62 dB; 地物统计水体区域的均值从-22.65 dB提升到-21.97 dB, 农田区域的均值从-14.45 dB提升到-14.22 dB; 2)零值处理后的影像与Sentinel-1影像的地物统计结果差异减小, 说明零值处理之后的影像消除了无效像素值的影响, 修正了影像值的统计偏差。

2) 几何校正实验

经过对高分三号大量L1A级数据进行处理发现部分高分三号数据的元数据信息存在自洽性问题：1)高分三号影像元数据4个角点经纬度坐标无法与RPC参数反算结果保持一致性。如表 3所示, 部分高分三号影像元数据文件提供的角点坐标经过RPC模型获得的影像像素坐标与真实像素坐标之间存在较大误差; 2)高分三号影像RPC参数内部存在不自洽, 导致Envi软件在利用RPC参数进行几何纠正时, 部分影像无法完成计算或处理结果出错。

从表 3可以看到, 高分三号不同升降轨方向或左右视方向部分存在元数据文件参数和RPC参数不自洽的问题。图上最大误差达到10 039.0像素, 对应纬度误差为0.136°。因此, 通过高分三号影像元数据文件中提供的角点坐标参数对影像直接进行RPC正算, 无法得到正确的纠正结果。由此可见, 利用RPC反算来获取影像角点坐标可以保证RPC参数内部的自洽性, 对纠正流程的稳健性构成充分条件, 保证处理流程的一致性。

针对存在的第2类自洽性问题, 即Envi软件对部分高分三号影像进行几何校正时无法得到正确结果, 选取2景高分三号影像分别使用Envi软件和本文方法进行实验, 影像的参数和Envi软件错误说明如表 4所示。Envi对第1景影像无法完成处理, 对第2景的处理结果如图 4所示。本文方法得到的对应处理结果见图 5(a)和5(b)。

从图 5可看出, 对于Envi软件仅依据RPC参数无法计算的影像和结果扭曲的影像, 本文提出的方法均可得到正确的几何校正结果。这是由于本文方法通过引入元数据文件中卫星的轨道方向和视向参数(见式(2))来保证RPC参数的自洽性, 使本文提出的几何校正流程具有高鲁棒性。

Sentinel-2影像95 % 以上像素定位误差 < 11 m^[16], 为评价本文方法的影像纠正结果(即图 5)精度, 选用Sentinel-2光学影像作为基准数据。分别在图 5(a)和5(b)影像中选取20个SAR影像和光学影像上的同名点作为检查点。同名点之间的坐标误差作为影像几何校正结果的精度误差。得到的同名点误差分布如图 6所示。

针对图 5(a)所得的5 m分辨率纠正结果, 为避免高程影响, 在北京平原地区选取SAR影像与光学影像中20个地物同名点为检查点, 同名点之间的平面坐标偏差作为几何校正后的误差。通过高分三号卫星平台的轨道倾角(10°), 将平面坐标误差旋转到成像坐标系下(距离向、方位向), 得到20个同名点在距离向和方位向上的误差分布如图 6(a)所示。选取的北京平原地区高分三号FSI模式影像所有同名点的误差均方根为118.98 m, 通过与文献[14]得到的高分三号影像京津冀地区FSI模式影像几何校正结果误差109.428~116.290 m相比, 可以认为, 利用本文提出的几何校正方法来处理Envi无法处理的图像, 其所得结果也在合理范围之内。相较于同一地区的Radarsat-2影像无控制点RPC几何校正后误差呈现出均匀分布^[17], 而图 6(a)中同名点误差的聚类中心却位于(121.00, -2.10)处, 可推定此高分三号影像在距离向上存在约121 m的系统差。针对图 5(b)广东地区SS模式影像所得的25 m分辨率纠正结果, 所选同名点的误差均方根为862.89 m。分析北京地区和广东地区影像纠正结果的差异可能由以下原因导致：1)2景影像区域内的高程差不同。广东地区影像大部分属于丘陵地形, 高程起伏较大, 在选取的影像同名点分布范围内的影像高程均值为175.495 m, 高程差为1 102 m。而北京地区控制点选在平原区, 影像高程为101.934 m, 高程差为230 m。高程差异导致了斜视成像的SAR影像几何校正精度的差异。2)分辨率差异。广东地区SS模式影像分辨率为25 m, 低于北京地区FSI影像的5 m分辨率。这导致在粗分辨率影像上选取同名检查点较为困难, 且定位误差较大。综合以上分析, 本文提出的处理方法可避免高分三号影像数据RPC参数不自洽的问题, 且获得正确的几何校正结果。

3 结论

本文提出的高分三号L1A级影像精处理方法, 可以有效地去除原始影像上低于等效噪声系数的无效值引起的斑噪影响, 提高影像视觉效果。从辐射强度统计量的变化上也说明本文提出的辐射校正流程消除了影像中无效值对影像地物辐射强度的统计偏差。本文提出的几何校正流程具有高稳健性, 可以消除高分三号元数据信息不自洽问题, 处理高分三号影像数据中存在的利用Envi软件无法完成几何校正或处理结果扭曲的问题。同时, 本文提出的将Lee滤波方法替代常用的最近邻和双线性重采样插值方法, 可在实现L1A产品方位向过采样的同时实现斑噪抑制效果。

综上所述, 本文提出的高分三号L1A级数据精处理方法对辐射校正处理结果和几何校正处理结果具有纠正后精度提升和处理流程稳健性的特点。但对于如何消除高分三号影像自身存在的绝对辐射值稍低和绝对几何定位精度稍差的问题, 还有待进一步研究。

感谢中国资源卫星应用中心提供的多模式、多分辨率高分三号L1A级SAR影像, 用以验证本文所提出的方法。