CN103218780B

CN103218780B - 基于逆rd定位模型的无控星载sar图像正射校正方法

Info

Publication number: CN103218780B
Application number: CN201310101077.3A
Authority: CN
Inventors: 刘佳音; 仇晓兰; 胡玉新; 雷斌; 付琨
Original assignee: Institute of Electronics of CAS
Current assignee: Jigang Defense Technology Co ltd; Aerospace Information Research Institute of CAS
Priority date: 2013-03-27
Filing date: 2013-03-27
Publication date: 2015-09-02
Anticipated expiration: 2033-03-27
Also published as: CN103218780A

Abstract

本发明提供一种基于逆RD定位模型的无控星载SAR图像正射校正方法，不仅解决了现有技术方案中需要迭代获取地面高程数值的难题，而且使SAR图像的正射校正流程易于并行编程，在保证几何校正精度基础上提高处理效率。首先利用像素对应的卫星轨道位置，构建斜距多普勒方程组，获得四角点的地理经纬度；然后地图投影，获得四角点的指东指北坐标和所需图像地图范围；其次进行逆地图投影获得像元的地理经纬度；按照公式转换为像元在地固坐标系统下的位置矢量；最后利用双线性插值法进行重采样所需图像像元的灰度值；直到所需图像地图范围内全部像元计算完毕，获得对应的灰度值，由此获得校正后图像。

Description

基于逆RD定位模型的无控星载SAR图像正射校正方法

技术领域

本发明是一种基于逆RD定位模型的无控星载SAR图像正射校正方法，涉及合成孔径雷达（SAR）图像处理及图像几何校正处理领域。

背景技术

随着高分辨率新型SAR传感器卫星（例如TerraSAR-X，RadarSAT-2等）的成功发射，其在国民经济建设中所发挥的作用越来越突出。但是由于SAR特殊的斜距成像机理，由地形起伏引起的几何失真大，定位误差严重，从而导致定量几何分析及影像解译困难，这在很大程度上限制了其在测绘、林业、地质灾害预防等领域中的应用。因此，对于实现自动化大批量星载SAR正射校正方法的研究有着广泛的应用需求和较大的技术难度。目前星载SAR正射校正方法分为有控校正和无控校正两大类。有控校正是在地面控制点的基础上，通过多项式、共线方程、有理函数等模型进行星载SAR图像的正射校正。该类方法在地面高程数据基础上还需要一定的地面控制点信息，不但限制了SAR图像正射校正的制作，而且由于控制点匹配等技术难题难以实现自动化大批量的星载SAR图像正射校正。

在无地面控制点情况下实现正射校正的方法也称为系统级正射校正，该类方法基于卫星轨道数据、SAR系统参数、成像处理参数以及地面高程数据对SAR图像实现正射校正。

目前，系统级正射校正方法普遍通过利用DEM模拟SAR影像来完成，首先根据DEM数据采用距离-多普勒模型和经验公式模拟出SAR图像，然后分别采用Harris算子和互信息匹配的方法提取模拟图像和实际SAR图像上的同名特征点，并根据这些特征点构建的不规则三角网（TIN）实现模拟SAR图像和实际SAR图像的精确配准，最后将星载SAR图像通过实际图像和模拟图像的精确配准关系以及模拟图像和DEM数据之间的对应关系校正到DEM所在的地理坐标中，实现SAR图像的正射校正[尤红建，《组合图像模拟和精配准的星载SAR图像正射校正》，测绘科学，2009年11月，第34卷第6期]。

系统级正射校正方法还利用迭代求解由斜距方程、多普勒方程以及地球模型构建的非线性方程组来实现星载SAR图像正射校正[J.C.Curlander,R.Kwork A post-processing Systemfor Automated Recitification and Registration of Space-borne SAR Imagery,Int.J.Remote Sensing,1987:621-638]。该方法首先根据卫星轨道数据、SAR系统参数以及成像处理结果构建斜距-多普勒-地球模型的非线性方程组，如下：

1.斜距方程：

R = | \overset{&RightArrow;}{R_{s}} - \overset{&RightArrow;}{R_{t}} |

2.多普勒中心频率方程：

f_{d} = \frac{- 2}{κR} (\overset{&RightArrow;}{R_{s}} - \overset{&RightArrow;}{R_{t}}) \cdot (\overset{&RightArrow;}{V_{s}} - \overset{&RightArrow;}{V_{t}})

3.地球模型方程：

\frac{x_{t}^{2} + y_{t}^{2}}{{(R_{e} + h)}^{2}} + \frac{z_{t}^{2}}{{(R_{p} + h)}^{2}} = 1

然后根据平均高程求解出像素位置并提取对应像素位置的高程数值，接着重新将新高程数值代入上述非线性方程组中，计算出新的像素位置，通过若干次迭代后，得到准确的像素位置。将SAR按照一定大小分成格网，格网每个角点重复上述迭代求解的过程。最后利用多项式模型或有理函数模型建立转换关系，并进行图像重采样，获得正射校正图像。

利用DEM数据模拟SAR图像的方法以SAR构象几何关系为基础模拟出SAR图像作为参考影像，然后根据较为成熟的有控校正方法进行星载SAR图像正射校正，但是存在以下缺点：

1.需要经验公式来模拟SAR散射特性，而经验公式缺乏普适性，难以用于工程化的批量生产；

2.利用DEM数据模拟SAR图像只能适应地形起伏剧烈的山区，对于既有平坦区域又有地形起伏的复杂情况的宽覆盖面积的星载SAR图像，则校正后定位精度较差。

利用迭代求解由斜距-多普勒-地球模型构建的非线性方程组的方法建立在SAR构象方程基础上，不需要任何经验数值，可以实现自动化生产，但是存在以下问题：

1.反复迭代求解过程难以进行并行化编程，运行效率低；

2.为降低运算量，格网内的点用多项式模型代替，不符合SAR严密几何模型，校正后几何畸变较大。

发明内容

本发明不仅解决了现有技术方案中需要迭代获取地面高程数值的难题，而且使SAR图像的正射校正流程易于并行编程，在保证几何校正精度基础上提高处理效率。

本发明的技术方案如下：

一种基于逆RD定位模型的无控星载SAR图像正射校正方法，包括以下步骤：

第一步：对卫星轨道数据以像素为自变量进行多项式拟合，其中多项式拟合模型如下：

\{\begin{matrix} x_{s} (m) = a_{0} + a_{1} m + a_{2} m^{2} + a_{3} m^{3} \\ y_{s} (m) = b_{0} + b_{1} m + b_{2} m^{2} + b_{3} m^{3} \\ z_{s} (m) = c_{0} + c_{1} m + c_{2} m^{2} + c_{3} m^{3} \end{matrix}

其中，m表示像素号，x_s,y_s,z_s表示该像素对应的卫星轨道位置x方向、y方向、z方向数值；

第二步：利用第一步得到的像素对应的卫星轨道位置，构建斜距多普勒方程组，并迭代求解，从而对SAR图像四角点进行像素定位，获得四角点的地理经纬度；

第三步：将第二步得到的四角点的地理经纬度进行地图投影，获得四角点的指东指北坐标；

第四步：根据第三步得到的四角点指东指北坐标，获得所需图像地图范围；

第五步：根据地图分辨率对所需图像像元进行逆地图投影，获得像元的地理经纬度；

第六步：根据所需图像像元的地理经纬度，在数字高程数据库（DEM）中提取对应的高程数值；

第七步：根据第五步得到的像元的地理经纬度和第六步得到的高程数值，按照下面的公式转换为像元在地固坐标系统下的位置矢量；

\{\begin{matrix} x_{i, j} = (\frac{R_{e}}{\sqrt{1 - e^{2} \sin^{2} ({lat}_{i, j})}} + h_{i, j}) \cos ({lat}_{i, j}) \cos ({lon}_{i, j}) \\ y_{i, j} = (\frac{R_{e}}{\sqrt{1 - e^{2} \sin^{2} ({lat}_{i, j})}} + h_{i, j}) \cos ({lat}_{i, j}) \sin ({lon}_{i, j}) \\ z_{i, j} = (\frac{R_{e}}{\sqrt{1 - e^{2} \sin^{2} ({lat}_{i, j})}} \cdot (1 - e^{2}) + h_{i, j}) \sin ({lat}_{i, j}) \end{matrix}

其中R_e为地球长半轴，e为地球第一偏心率，lat_i,j为像元对应的纬度，lon_i,j为像元对应的经度，h_i,j为像元对应的高程数值；

第八步：根据第一步得到的像素对应卫星轨道位置，构建像元对应的逆斜距多普勒方程组，方程如下，并采用牛顿法求解逆斜距多普勒方程组：

斜距方程：

F_{1} = \sqrt{{(x_{s} (m) - x_{t})}^{2} + {(y_{s} (m) - y_{t})}^{2} + {(z_{s} (m) - z_{t})}^{2}} - (R_{0} + \frac{c_{light}}{2 \cdot f_{s}} \cdot n)

多普勒中心频率方程：

F_{2} = 2 \cdot \frac{x_{s} (m) \cdot {vx}_{s} (m) + y_{s} (m) \cdot {vy}_{s} (m) + z_{s} (m) \cdot {vz}_{s} (m)}{\sqrt{{(x_{s} (m) - x_{t})}^{2} + {(y_{s} (m) - y_{t})}^{2} + {(z_{s} (m) - z_{t})}^{2}}} + λ \cdot f_{d} (n)

其中，R₀为初始斜距值，f_s为采样频率，f_d(n)为像素对应的多普勒中心频率；

第九步：利用双线性插值法进行重采样所需图像像元的灰度值；

第十步：重复第五步至第九步，直到所需图像地图范围内全部像元计算完毕，获得对应的灰度值，由此获得校正后图像。

本发明的有益效果：

本发明不需要进行图像模拟，避免了仿真SAR散射特性的瓶颈；同时利用逆斜距-多普勒模型，既符合SAR严密几何模型，降低了图像畸变，又无需迭代，易于计算机并行处理，极大提高了校正效率。

附图说明

图1本发明的正射校正方法流程图。

具体实施方式

选取一幅北京地区的星载SAR图像产品，结合附图说明本发明的具体实施方式：

首先从产品中提取出本发明所需的相关参数。包括以下参数以及若干组卫星位置速度矢量。

第一步：对卫星轨道数据以像素为自变量进行多项式拟合，获得卫星位置速度和图像方位像素的关系；

a₀	-1882349．705302733	b₀	5030621．015380397	c₀	4289948．518499875
						a₁	3211．325617072	b₁	-3811．343856856	c₁	5858．518584492
a₂	1．298738241	b₂	-3．181818766	c₂	-2．587766992
						a₃	-0．047247270	b₃	-0．015512607	c₃	-0．011412460

指东坐标最大值	468417.414424
		指东坐标最小值	434944.323800
指北坐标最大值	4404882.211167
		指北坐标最小值	4438471.615221

\{\begin{matrix} x_{i, j} = (\frac{R_{e}}{\sqrt{1 - e^{2} \sin^{2} ({lat}_{i, j})}} + h_{i, j}) \cos ({lat}_{i, j}) \cos ({lon}_{i, j}) \\ y_{i, j} = (\frac{R_{e}}{\sqrt{1 - e^{2} \sin^{2} ({lat}_{i, j})}} + h_{i, j}) \cos ({lat}_{i, j}) \sin ({lon}_{i, j}) \\ z_{i, j} = (\frac{R_{e}}{\sqrt{1 - e^{2} \sin^{2} ({lat}_{i, j})}} \cdot (1 - e^{2}) + h_{i, j}) \sin ({lat}_{i, j}) \end{matrix}

斜距方程：

F_{1} = \sqrt{{(x_{s} (m) - x_{t})}^{2} + {(y_{s} (m) - y_{t})}^{2} + {(z_{s} (m) - z_{t})}^{2}} - (R_{0} + \frac{c_{light}}{2 \cdot f_{s}} \cdot n)

多普勒中心频率方程：

F_{2} = 2 \cdot \frac{x_{s} (m) \cdot {vx}_{s} (m) + y_{s} (m) \cdot {vy}_{s} (m) + z_{s} (m) \cdot {vz}_{s} (m)}{\sqrt{{(x_{s} (m) - x_{t})}^{2} + {(y_{s} (m) - y_{t})}^{2} + {(z_{s} (m) - z_{t})}^{2}}} + λ \cdot f_{d} (n)

其中，R₀为初始斜距值，f_s为采样频率，f_d(n)为像素对应的多普勒中心频率。

牛顿法的具体求解内核如下：

F_{1} (m_{0}, n_{0}) + \frac{{&PartialD; F}_{1}}{&PartialD; m} dm + \frac{{&PartialD; F}_{1}}{&PartialD; n} dn = 0

F_{1} (m_{0}, n_{0}) + \frac{{&PartialD; F}_{2}}{&PartialD; m} dm + \frac{{&PartialD; F}_{2}}{&PartialD; n} dn = 0

\frac{{&PartialD; F}_{1}}{&PartialD; m} = \frac{1}{2} {({(x_{s} (m) - x_{t})}^{2} + {(y_{s} (m) - y_{t})}^{2} + {(z_{s} (m) - z_{t})}^{2})}^{- 0.5} \cdot (2 (x_{s} (m) - x_{t}) \frac{{&PartialD; x}_{s} (m)}{&PartialD; m}

2 + (y_{s} (m) - y_{t} \frac{&PartialD; y_{s} (m)}{&PartialD; m} + 2 (z_{s} (m) - z_{t}) \frac{&PartialD; z_{s} (m)}{&PartialD; m})

\frac{{&PartialD; F}_{z}}{&PartialD; m} = 2 \frac{{vx}_{z} (m) \frac{{&PartialD; x}_{z} (m)}{&PartialD; m} + x_{z} (m) \frac{{&PartialD; vx}_{z} (m)}{&PartialD; m} + {vy}_{z} (m) \frac{{&PartialD; y}_{z} (m)}{&PartialD; m} + y_{z} (m) \frac{{&PartialD; vy}_{z} (m)}{&PartialD; m} + {vz}_{z} (m) \frac{{&PartialD; x}_{z} (m)}{&PartialD; m} + z_{z} (m) \frac{&PartialD; {vz}_{z} (m)}{&PartialD; m}}{\sqrt{{(x_{z} (m) - x_{t})}^{z} + {(y_{z} (m) - y_{t})}^{z} + {(z_{z} (m) - z_{t})}^{z}}} - (x_{z} (m) \cdot {vx}_{z} (m) + y_{z} (m) \cdot {vy}_{z} (m)

+ z_{z} (m) \cdot {vz}_{z} (m)) \cdot {((x_{z} (m) - x_{t})}^{z} + {(y_{z} (m) - y_{t})}^{z} + {(z_{z} (m) - z_{t})^{z})}^{- 2.5}

\frac{{&PartialD; F}_{1}}{&PartialD; n} = - \frac{c_{light}}{2 \cdot f_{s}}

\frac{{&PartialD; F}_{2}}{&PartialD; n} = λ (d_{1} + 2 d_{2} \cdot n)

\frac{{&PartialD; x}_{s} (m)}{&PartialD; m} = a_{1} + {2 a}_{2} \cdot m + {3 a}_{3} \cdot m^{2}

\frac{{&PartialD; y}_{s} (m)}{&PartialD; m} = b_{1} + {2 b}_{2} \cdot m + {3 b}_{3} \cdot m^{2}

\frac{{&PartialD; z}_{s} (m)}{&PartialD; m} = c_{1} + {2 c}_{2} \cdot m + {3 c}_{3} \cdot m^{2}

\frac{{&PartialD; vx}_{s} (m)}{&PartialD; m} = ({2 a}_{2} + {6 a}_{3} \cdot m) \cdot PRT

\frac{{&PartialD; vy}_{s} (m)}{&PartialD; m} = (2 b_{2} + {6 b}_{3} \cdot m) \cdot PRT

\frac{{&PartialD; vz}_{s} (m)}{&PartialD; m} ({2 c}_{2} + {6 c}_{3} \cdot m) \cdot PRT

其中，PRT为脉冲重复时间。

第九步：将第八步求解的SAR图像像素号，利用双线性插值法进行重采样所需图像像元的灰度值；

下面分别对全景图像和以山区为主的部分图像，采用现有技术方案及本发明提供的方法进行正射校正，其中现有技术方案主要有模拟SAR图像方法和迭代求解斜距-多普勒-地球模型方程的方法。由于迭代求解斜距-多普勒-地球模型方程的方法无法对图像内全部像元逐点进行构建方程并迭代求解，因此分别在图像上分别对5*5,10*10,25*25以及100*100个像元进行构建方程并迭代求解，其余图像像元利用上述点采用多项式模型进行逼近。

表1对1000*1000的图像（3公里*3公里覆盖区域）进行正射校正结果

表2对10000*10000的图像（30公里*30公里覆盖区域）进行正射校正结果

表1和表2分别给出对以山区为主的小幅图像（即1000*1000的图像（3公里*3公里覆盖区域））和全景大幅图像（即10000*10000的图像（30公里*30公里覆盖区域））进行正射校正的结果。从结果可以看出，模拟SAR图像方法只能适应以山区为主的小幅图像，对于复杂地形的大幅图像，模拟的图像纹理与真实SAR图像区别较大，导致匹配失败，无法进行正射校正；迭代求解斜距-多普勒-地球模型方程的方法无法对图像内全部像元逐点进行构建方程并迭代求解，当选取的像元点少时（25个像元或100个像元）处理速度虽然高，但由于个别异常点会导致定位精度低，几何畸变大的问题，当选取的像元多时（625个像元或1万个像元）处理速度很慢；而本发明方法不但适用于小幅图像，而且对大幅图像也能进行逐点斜距-多普勒模型计算，通过采用多个CPU并行处理，本发明对于全景大幅图像（30公里*30公里覆盖区域）只需要150秒就可以处理完成，而且定位精度高（35..2米），几何畸变低（3.1米）。

如卫星轨道数据拟合采用的多项式，可以是三阶也可以是更高阶，只要采用多项式模型对卫星轨道数据进行拟合，均属于本发明的范畴之内。

如利用双线性插值实现像素重采样，也可以利用最近邻法或双三次卷积法等常见图像重采样方法，均属于本发明的范畴之内。

Claims

1.一种基于逆RD定位模型的无控星载SAR图像正射校正方法，其特征在于，包括以下步骤：

\{\begin{matrix} y_{s} (m) = a_{0} + a_{1} m + a_{2} m^{2} + a_{3} m^{3} \\ y_{s} (m) = b_{0} + b_{1} m + b_{2} m^{2} + b_{3} m^{3} \\ z_{s} (m) = c_{0} + c_{1} m + c_{2} m^{2} + c_{3} m^{3} \end{matrix}

其中，m表示像素号，x_s，y_s，z_s表示该像素对应的卫星轨道位置x方向、y方向、z方向数值；

第六步：根据所需图像像元的地理经纬度，在数字高程数据库(DEM)中提取对应的高程数值；

\{\begin{matrix} x_{i, j} = (\frac{R_{e}}{\sqrt{1 - e^{2} \sin^{2} ({lat}_{i, j})}} + h_{i, j}) \cos ({lat}_{i, j}) \cos ({lon}_{i, j}) \\ y_{i, j} = (\frac{R_{e}}{\sqrt{1 - e^{2} \sin^{2} ({lat}_{i, j})}} + h_{i, j}) \cos ({lat}_{i, j}) \sin ({lon}_{i, j}) \\ z_{i, j} = (\frac{R_{e}}{\sqrt{1 - e^{2} \sin^{2} ({lat}_{i, j})}} \cdot (1 - e^{2}) + h_{i, j}) \sin ({lat}_{i, j}) \end{matrix}

其中R_e为地球长半轴，e为地球第一偏心率，lat_i，j为像元对应的纬度，lon_i，j为像元对应的经度，h_i，j为像元对应的高程数值；

斜距方程：

F_{1} = \sqrt{{(x_{s} (m) - x_{t})}^{2} + {(y_{s} (m) - y_{t})}^{2} + {(z_{s} (m) - z_{t})}^{2}} - (R_{0} + \frac{(c_{light})}{2 \cdot f_{s}} \cdot n)

多普勒中心频率方程：

F_{2} = 2 \cdot \frac{x_{s} (m) \cdot {vx}_{s} (m) + y_{s} (m) \cdot {vy}_{s} (m) + z_{s} (m) \cdot {vz}_{s} (m)}{\sqrt{{(x_{s} (m) - x_{t})}^{2} + {(y_{s} (m) - y_{t})}^{2} + {(z_{s} (m) - z_{t})}^{2}}} + λ \cdot f_{d} (n)