CN105510913A

CN105510913A - 基于类光学像方改正的异源光学和sar遥感影像联合定位方法

Info

Publication number: CN105510913A
Application number: CN201510767975.1A
Authority: CN
Inventors: 吴颖丹; 明洋; 黄楚云; 吕辉; 叶志伟; 雷勇; 朱莹
Original assignee: Hubei University of Technology
Current assignee: Hubei University of Technology
Priority date: 2015-11-11
Filing date: 2015-11-11
Publication date: 2016-04-20
Anticipated expiration: 2035-11-11
Also published as: CN105510913B

Abstract

一种基于类光学像方改正的异源光学和SAR遥感影像联合定位方法，包括：1）利用严格几何模型建立立体格网，计算每个格网点在SAR遥感影像上对应的像点坐标，建立虚拟控制点格网，解算SAR遥感影像的RPC参数；2）获取光学和SAR遥感影像间连接点和地面控制点对应像点的像素坐标；3）建立对应的像方误差改正模型，利用基于类光学像方改正的异源光学和SAR遥感影像RFM区域网平差方法，实现精确联合定位。本发明采用类光学像方改正模型对其系统误差进行改正，在近似统一的平差模型下实现异源光学和SAR遥感影像的联合定位，可以取得与严格几何模型相当甚至更优的定位精度，极大降低对地面控制点的要求。

Description

基于类光学像方改正的异源光学和SAR遥感影像联合定位方法

技术领域

本发明属于遥感影像处理技术领域，涉及一种基于类光学像方改正的异源光学和SAR遥感影像联合定位方法。

背景技术

合成孔径雷达(SyntheticApertureRadar，SAR)，是一种是用微波探测地表目标的主动式成像传感器。星载SAR，具有全天候、全天时、大范围、强穿透性等独特优点，其获取的SAR遥感影像，可作为光学遥感影像的有益补充，已广泛地应用于地形测绘、国土资源调查、农作物估产、水文监测和灾害评估等领域。异源光学和SAR遥感影像精确联合定位，是诸多应用的重要前提和关键环节，对其研究具有重要的现实意义。

就几何定位而言，构成广义立体像对实现对地目标定位的充要条件是“不同位置拍摄”和“影像存在一定重叠”，与影像是否属于同类传感器无关。利用异源光学和SAR遥感影像进行联合定位，不仅可以降低遥感影像立体成像的要求，而且可以充分利用单景异源遥感影像的互补信息，为此，已有较多学者对其进行深入研究。目前，研究方法主要可分为如下三类：

1)先空间后方交会-再前方交会的方法，即首先利用大量均匀分布的地面控制点分别对每景影像进行定向，然后利用前方交会技术进行定位，如2004年Toutin学者利用其提出的几何模型对Landsat-7ETM+、SPOT4HRV、ASTERVNIR、RADARSAT、ERS等多源影像进行了联合定位试验、2008年邢帅等人尝试将光学与SAR卫星遥感影像联立构成“立体像对”进行复合式定位等。

2)基于严格几何模型的区域网平差方法，如2006年Toutin学者对Landsat-7ETM+、SPOT4HRV、ASTERVNIR、RADARSAT、ERS等影像，采用25个控制点进行平差试验，2008年邢帅对多星多源传感器SPOT、IKONOS、ERS及其RADARSAT影像采用严格几何模型进行了联合区域网平差实验研究等。

3)基于有理函数模型(RationalFunctionModel，RFM)的区域网平差方法，如秦绪文、张过等人直接将RFM模型用于多源遥感影像的区域网平差，并取得较好的定位结果。

然而，上述方法在实际应用中主要存在如下几点问题：

1)先空间后方交会-再前方交会以及基于严格几何模型的区域网平差方法，对地面控制点要求较高。采用第一种方法时，对每景影像进行定向，均需要野外布设大量、均匀分布的地面控制点；而采用第二种方法时，可一定程度减少所需的控制点数量，但对控制点的分布要求较高，一般需尽量四角布设。

2)基于严格几何模型的区域网平差方法，求解比较复杂。采用第一、二种方法时，当涉及不同类型传感器影像的联合定位时，需要利用不同传感器对应的严格几何模型建立误差方程进行求解，涉及较多复杂的物理参数，求解比较复杂。

3)基于RFM的区域网平差模型，虽具有与传感器无关、形式简单、便于计算等优点，但其像方改正模型主要针对于光学遥感影像。而对于SAR遥感影像，则直接沿用光学遥感影像的仿射变换误差补偿模型。尽管SAR遥感影像亦具有视场角窄、飞行高度高的特点，使其定位误差均具有一定系统性，但由于SAR侧视斜距投影的成像方式，造成地形起伏引起的像点位移较光学影像会严重几倍甚至几十倍，这在某种意义上将破坏误差系统性，从而使得用于光学影像的误差补偿模型不再完全适用于SAR遥感影像。

发明内容

本发明针对目前异源光学和SAR遥感影像联合定位方法存在的缺陷，提出了一种基于类光学像方改正的异源光学和SAR遥感影像联合定位方法，可在有限地面控制点条件下，实现异源光学和SAR遥感影像的精确几何定位，大大降低了野外工作量，具有显著的经济和社会效益。

为了达到上述目的，本发明采用如下技术方案：

一种基于类光学像方改正的异源光学和SAR遥感影像联合定位方法，其特征在于：所述基于类光学像方改正的异源光学和SAR遥感影像联合定位方法包括以下步骤：

1)解算异源光学和SAR遥感影像的RPC参数；

2)利用影像匹配或者人工量测的方法，获取光学和SAR遥感影像间连接点和地面控制点对应像点的像素坐标；

3)待完成步骤2)后，基于步骤1)所解算得到的异源光学和SAR遥感影像的RPC参数以及基于类光学像方改正的异源光学和SAR遥感影像区域网平差，实现精确联合定位。

作为优选，本发明所采用的步骤1)的具体实现方式是：

1.1)构建SAR遥感影像在物方均匀分布的虚拟控制点格网；

1.2)根据步骤1.1)所建立的虚拟控制点格网求解SAR遥感影像有理函数多项式模型；

1.3)对步骤1.2)所建立得到的有理函数多项式模型进行线性化，同时采用谱修正迭代算法计算得到异源光学和SAR遥感影像的RPC参数。

作为优选，本发明所采用的步骤1.1)的具体实现方式是：

1.1.1)根据SAR遥感影像的四个角点的像素坐标求解对应物方点的大地经纬度坐标，确定SAR遥感影像的覆盖范围；

1.1.2)读取SAR遥感影像的覆盖范围内的免费全球DEM数据，确定最大高程值以及最小高程值，并根据给定的高程分层数建立三维立体格网；

1.1.3)对步骤1.1.2)中所建立得到的三维立体格网中各格网点，由其大地经纬度坐标计算各格网点在SAR遥感影像上对应的像点坐标，并进一步对像点坐标进行从侧视斜距投影到光学透视投影的转换，建立物方均匀分布的虚拟控制点格网。

作为优选，本发明所采用的步骤1.1.1)中求解对应物方点的大地经纬度坐标的具体实现方式是：

由某像点像素坐标(x,y)和人工设定的投影面高程h，利用SAR严格几何模型并引入地球椭球模型，计算该像点对应物方点的地心旋转坐标系坐标(X,Y,Z)，通过坐标转换可得到对应的大地经纬度坐标(lat,lon,h)；

所述SAR严格几何模型是根据距离条件和多普勒频率条件来表达雷达图像瞬时构像的几何关系：

R_{S l a n t} = | R_{s} - R_{p} | = \sqrt{{(X - X_{S})}^{2} + {(Y - Y_{S})}^{2} + {(Z - Z_{S})}^{2}}

(X-X_s)V_SX+(Y-Y_s)V_SY+(Z-Z_s)V_SZ＝0

其中：

R_Slant为卫星天线中心S到地面点P的斜距；

R_p以及V_p分别为地面点P的物方坐标(X,Y,Z)以及速度矢量(V_px,V_py,V_pz)；

R_s以及V_s分别为卫星天线中心S成像瞬间对应的物方空间坐标(X_S,Y_S,Z_S)以及速度矢量(V_SX,V_SY,V_SZ)；

所述地球椭球模型的函数表达式是：

\frac{(X^{2} + Y^{2})}{{(r_{e} + h)}^{2}} + \frac{Z^{2}}{r_{p}^{2}} = 1

其中：

r_e,r_p分别为地球椭球的长、短半轴；

h为人工设定的投影面高程；

在计算过程中，首先从SAR遥感影像头文件中读取所需的模型参数，然后以物方点的地心旋转坐标系坐标(X,Y,Z)为未知数，对上述公式线性化，建立误差方程的矩阵形式为：

V＝BX-L,P

式中：

B = [\begin{matrix} B_{11} & B_{12} & B_{13} \\ B_{21} & B_{22} & B_{23} \\ B_{31} & B_{32} & B_{33} \end{matrix}]

为(X,Y,Z)的改正数对应的系数矩阵；

X为(X,Y,Z)的改正数向量；

L＝-[F₁F₂F₃]^T为常数向量；

P为权矩阵；

根据误差方程，按照最小二乘法迭代求解解算出(X,Y,Z)，经坐标转换即可得到大地经纬度坐标(lat,lon,h)。

作为优选，本发明所采用的步骤1.1.3)中由其大地经纬度坐标计算各格网点在SAR遥感影像上对应的像点坐标的具体实现方式是：

由某物方点的大地经纬度坐标(lat,lon,h)，转换为地心旋转坐标系坐标(X,Y,Z)，通过SAR严格几何模型的间接定位公式，计算该物方点在SAR遥感影像上对应像点的像素坐标(x,y)：

R_{S l a n t} = | R_{s} - R_{p} | = \sqrt{{(X - X_{S})}^{2} + {(Y - Y_{S})}^{2} + {(Z - Z_{S})}^{2}} - - - (1)

f_{D} = - \frac{2 (R_{s} - R_{p}) \cdot (V_{s} - V_{p})}{λ | R_{s} - R_{p} |} - - - (2)

其中：

R_Slant为卫星天线中心S到地面点P的斜距；

f_D为多普勒频率；

λ为卫星天线所接收的雷达的波长；

其中，斜距R_Slant为距离向坐标y的函数；

对于斜距影像，R_Slant＝R₀+yM_y，R₀为第一斜距，M_y为斜距影像的距离向分辨率，参数R₀和M_y可以从SAR遥感影像头文件中直接读取；y为像点的像素坐标(x,y)在斜距向的值；

对于SAR遥感影像，方位向上每一行对应不同的卫星位置和速度矢量，表示为成像时刻t的函数：

\{\begin{matrix} X_{s} = a_{0} + a_{1} t + a_{2} t^{2} \\ Y_{s} = b_{0} + b_{1} t + b_{2} t^{2} \\ Z_{s} = c_{0} + c_{1} t + c_{2} t^{2} \\ V_{s x} = a_{1} + 2 a_{2} t \\ V_{s y} = b_{1} + 2 b_{2} t \\ V_{s z} = c_{1} + 2 c_{2} t \end{matrix} - - - (3)

t＝t₀+(x-1)Δt

其中：

Δt为每行SAR遥感影像成像的时间间隔，Δt是根据SAR遥感影像的头文件提供的脉冲重复频率PRF得到，Δt＝1/PRF；

t为对应的成像时刻；

x为方位向坐标；

t₀为参考时刻；

a_i,b_i,c_i(i＝0,1,2)为多项式拟合系数，由SAR遥感影像的头文件中提供的若干个卫星状态向量拟合得到；

星载SAR遥感影像一般都进行了多普勒归零处理，且地物在地心旋转坐标系中速度为零，为此多普勒条件方程(2)可简化为：

(X-X_s)V_SX+(Y-Y_s)V_SY+(Z-Z_s)V_SZ＝0(4)

联立公式(1)和公式(4)即为距离多普勒数学模型；

对距离多普勒数学模型线性化，建立误差方程的矩阵形式为：

V＝AX-L,P(5)

其中：

A = [\begin{matrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{matrix}],

为未知数t、y改正数的系数；

X为未知数t、y改正数的向量；

L＝-[F₁F₂]^T为常数向量，其中，T表示矩阵转置，F1以及F2分别代表多普勒对应两个方程计算得到的值；

P为权矩阵；

根据误差方程(5)，按照最小二乘法迭代求解可解算出未知数(t,y)，进而由公式t＝t₀+(x-1)Δt可得到像点像素坐标(x,y)。

作为优选，本发明在步骤1.1.3)中对像点坐标进行从侧视斜距投影到光学透视投影的转换的具体实现方式是：

假设O点为卫星天线中心所在的位置，卫星天线中心距离整个测区平均高程面的高度均为H，O在地面上的投影为E点；A点为SAR遥感影像的近地点所对应的地面点，而B点为远地点所对应的地面点；OO'为SAR遥感影像成像面，DD'为等效焦距为f时的光学影像成像面；

则近地点A对应侧视角为：

θ₀＝arccos(H/R₀)

其中：

R₀为近地点斜距OA的距离；

H为卫星天线航高；

远地点B的斜距为：

OB＝R₀+M_xW

其中：

M_x为SAR遥感影像的斜距分辨率；

W为SAR遥感影像的宽度；

根据公式(7)，则有近地点A与远地点B之间的距离AB是：

A B = \sqrt{{OB}^{2} - H^{2}} - \sqrt{{OA}^{2} - H^{2}}

令近地点A与远地点B在等效焦距为f时的光学影像成像面DD'上的投影长度ab为W，则等效焦距f为：

f = \frac{H \times W}{A B}

设物方点P高程为h_p，物方点P在SAR遥感影像对应像点坐标为(x,y)，像点坐标(x,y)转换为透视中心投影坐标为(x',y')，则有：

\{\begin{matrix} x^{'} = \frac{f}{H - h_{p}} \sqrt{{(R_{0} + x \times M_{x})}^{2} - {(H - h_{p})}^{2}} - f t a n (θ_{0}) \\ y^{'} = y \end{matrix} .

作为优选，本发明所采用的步骤1.2)中有理函数多项式模型的具体定义是：

\{\begin{matrix} X = \frac{{Num}_{L} (P, L, H)}{{Den}_{L} (P, L, H)} \\ Y = \frac{{Num}_{s} (P, L, H)}{{Den}_{s} (P, L, H)} \end{matrix}

其中：

Num_L(P,L,H)、Den_L(P,L,H)、Num_s(P,L,H)以及Den_s(P,L,H)为三次多项式，三次多项式的系数为RPC参数；

(P,L,H)为正则化的地面坐标，(X,Y)为正则化的影像坐标，计算公式如下：

P = \frac{L a t - L A T_O F F}{L A T_S C A L E}

L = \frac{L o n - L O N G_O F F}{L O N G_S C A L E}

H = \frac{H e i g h t - H E I G H T_O F F}{H E I G H T_S C A L E}

X = \frac{S a m p l e - S A M P_O F F}{S A M P_S C A L E}

Y = \frac{L i n e - L I N E_O F F}{L I N E_S C A L E}

其中：

LAT_OFF、LAT_SCALE、LONG_OFF、LONG_SCALE、HEIGHT_OFF以及HEIGHT_SCALE为地面点坐标的正则化参数；

SAMP_OFF、SAMP_SCALE、LINE_OFF以及LINE_SCALE为影像坐标的正则化参数；

Lat表示纬度，Lon表示经度，Height表示高程，此三项即代表地面点的空间坐标；

Sample代表影像列坐标；Line代表影像行坐标。

作为优选，本发明所采用的步骤1.3)的具体实现方式是：

以RPC参数为未知数，对步骤1.2)所建立得到的有理函数多项式模型进行线性化，列立误差方程：

V＝BX-L,P

其中：

B = [\begin{matrix} \frac{\partial F_{x}}{\partial a_{i}} & \frac{\partial F_{x}}{\partial b_{i}} & \frac{\partial F_{x}}{\partial c_{i}} & \frac{\partial F_{x}}{\partial d_{i}} \\ \frac{\partial F_{y}}{\partial a_{i}} & \frac{\partial F_{y}}{\partial b_{i}} & \frac{\partial F_{y}}{\partial c_{i}} & \frac{\partial F_{y}}{\partial d_{i}} \end{matrix}], (i = 0, ..., 19)

L = [\begin{matrix} - F_{1}^{0} \\ - F_{2}^{0} \end{matrix}]

x＝[a_ib_ic_id_i]^T

P为权矩阵

采用谱修正迭代算法，计算得到SAR遥感影像的RPC参数。

作为优选，本发明所采用的步骤3)的具体实现方式是：

3.1)计算和给定各未知参数的初始值；

3.2)利用异源光学和SAR遥感影像的RPC参数，对步骤2)中所得到的异源光学和SAR遥感影像之间的连接点和地面控制点的每个连接点和地面控制点，以连接点和地面控制点对应的物方坐标，以及各图像对应的仿射变换系数为未知参数列立误差方程；

3.3)利用最小二乘原理，可解算出未知参数的改正数；

3.4)将步骤3.3)解算得到的未知参数的改正数与上次迭代的未知参数值相加，得到本次迭代更新后的未知参数值；当第一次迭代时，未知参数值为初始值；

3.5)利用更新后的未知参数值按步骤3.2)～3.4)重新进行计算，直至达到收敛条件，即完成异源光学和SAR遥感影像的联合定位。。

作为优选，本发明所采用的步骤3.2)的具体实现方式是：

对于光学遥感影像，直接采用像方仿射变换模型，即：

\{\begin{matrix} F_{x} = {px}_{0} + {px}_{1} \cdot S a m p l e + {px}_{2} \cdot L i n e + S a m p l e - x = 0 \\ F_{y} = {py}_{0} + {py}_{1} \cdot S a m p l e + {py}_{2} \cdot L i n e + L i n e - y = 0 \end{matrix} - - - (6)

其中：

(x,y)为物方点在光学遥感影像上对应像点的像素坐标；

(Sample,Line)为物方点利用所在光学遥感影像的RPC参数投影到光学遥感影像上的像素坐标；

(px_i,py_i)(i＝0,1,2)是待求光学遥感影像对应的仿射变换系数；其中，px₁、py₂的初始值为1，其他系数初始值均为0；

对于SAR遥感影像，采用光学遥感影像的像方仿射变换模型，即：

\{\begin{matrix} F_{x}^{'} = {px}_{0}^{'} + {px}_{1}^{'} \cdot s a m p l e + {px}_{2}^{'} \cdot l i n e + s a m p l e - x^{'} = 0 \\ F_{y}^{'} = {py}_{0}^{'} + {py}_{1}^{'} \cdot s a m p l e + {py}_{2}^{'} \cdot l i n e + l i n e - y^{'} = 0 \end{matrix} - - - (7)

其中：

(x′,y′)为对物方点在SAR遥感影像上对应像点的像素坐标进行透视中心投影转换后的像素坐标；

(sample,line)为利用所在SAR遥感影像的RPC参数计算的地面点投影到等效光学影像面上的像素坐标；

(px′_i,py′_i)(i＝0,1,2)是SAR遥感影像对应的仿射变换系数；

对于光学遥感影像上的每个连接点和地面控制点的像素坐标，按公式(6)进行线性化，列立误差方程：

V₁＝A₁t+B₁X-L₁,P₁(8)

对于SAR遥感影像上的每个连接点和地面控制点的像素坐标，按公式(7)进行线性化，列立误差方程：

V₂＝A₂t+B₂X-L₂,P₂(9)

对于地面控制点，将地面控制点的物方坐标视为带权观测值进行处理，还需列立如下误差方程：

V₃＝EX-L₃,P₃(10)

联立公式(8)～公式(10)，即为类光学像方改正的异源光学和SAR遥感影像区域网平差模型；

其中：

t为各影像对应的仿射变换系数改正数向量；

X为连接点以及地面控制点对应的物方坐标改正数向量；

A_i和B_i(i＝1,2)分别为未知数t和X对应的系数矩阵；

L_i(i＝1,2,3)为对应误差方程的常数向量；

P_i为(i＝1,2,3)为相应的权矩阵。

本发明的优点在于：

本发明主要分为如下两个部分：一是异源遥感影像RPC参数计算，二是基于类光学像方改正的异源光学和SAR遥感影像区域网平差。本发明首先对于SAR遥感影像，采用与地形无关的方法，利用严格几何模型建立立体格网，计算每个立体格网点在SAR遥感影像上对应的像点坐标，并进行从侧视斜距投影到光学透视投影的转换，根据每个立体格网点在SAR遥感影像上对应的像点坐标建立虚拟控制点格网，在此基础上，采用谱修正迭代算法，解算SAR遥感影像的RPC参数。其次，利用影像匹配或者人工量测的方法，获取光学以及SAR遥感影像间连接点和地面控制点对应像点的像素坐标。最后，针对不同的光学、SAR遥感影像，建立对应的像方误差改正模型，利用基于类光学像方改正的异源光学和SAR遥感影像RFM区域网平差方法，实现其精确联合定位。本发明针对SAR遥感图像，建立可顾及其系统误差特性的类光学像方改正模型，可以获得更优的定位精度。充分利用RFM模型传感器无关性，形式简单、解算方便等优点，建立了近似统一的异源光学和SAR遥感影像区域网平差模型。本方法具有RFM模型区域网平差方法受地面控制点分布影响较小的优点，对于地形复杂的困难地区，可大大降低地面控制点测量的难度。本发明为异源光学和SAR遥感影像的联合定位提供了一种新方法，创新性地提出了一种基于类光学像方改正的RFM模型区域平差方法，可充分顾及SAR遥感影像的误差特性，非常适用于缺少地面控制点的实际情况。实践表明，当地面控制点数量较少、分布不佳时，可以获得相比严格几何模型更优的定位精度。为此，本方法在借鉴RFM区域网平差模型方法的基础上，针对SAR遥感影像误差特性，提出了一种基于类光学像方改正的异源光学和SAR遥感影像联合定位方法，具有非常重要的现实意义。

附图说明

图1为SAR遥感影像像点从侧视斜距投影转换为透视中心投影的示意图。

具体实施方式

下面结合附图对本发明做进一步详细描述。

实施例1：

以下结合图1对本发明做进一步详细描述。一种基于类光学像方改正的异源光学和SAR遥感影像联合定位方法，各步骤详细阐述如下：

步骤1)异源光学和SAR遥感影像的RPC参数计算

目前，光学遥感影像均提供相应的RPC参数，可直接利用。下面将重点说明本方法中SAR遥感影像对应的RPC参数求解方法：

方式一)由某物方点的大地经纬度坐标(lat,lon,h)，转换为地心旋转坐标系坐标(X,Y,Z)，通过SAR严格几何模型的间接定位公式，计算该物方点在SAR遥感影像上对应像点的像素坐标(x,y)。

SAR严格几何模型，是根据距离条件和多普勒频率条件来表达雷达图像瞬时构像的几何关系，即：

R_{S l a n t} = | R_{s} - R_{p} | = \sqrt{{(X - X_{S})}^{2} + {(Y - Y_{S})}^{2} + {(Z - Z_{S})}^{2}} - - - (1)

f_{D} = - \frac{2 (R_{s} - R_{p}) \cdot (V_{s} - V_{p})}{λ | R_{s} - R_{p} |} - - - (2)

式中：

R_Slant为卫星天线中心S到地面点P的斜距；

f_D为多普勒频率；

λ为卫星天线所接收的雷达的波长；

其中，斜距R_Slant为距离向坐标y的函数；

对于SAR遥感影像，方位向上每一行对应不同的卫星位置和速度矢量，可表示为成像时刻t的函数：

\{\begin{matrix} X_{s} = a_{0} + a_{1} t + a_{2} t^{2} \\ Y_{s} = b_{0} + b_{1} t + b_{2} t^{2} \\ Z_{s} = c_{0} + c_{1} t + c_{2} t^{2} \\ V_{s x} = a_{1} + 2 a_{2} t \\ V_{s y} = b_{1} + 2 b_{2} t \\ V_{s z} = c_{1} + 2 c_{2} t \end{matrix} - - - (3)

t＝t₀+(x-1)Δt

其中：

Δt为每行SAR遥感影像成像的时间间隔，Δt是根据SAR遥感影像的头文件提供的脉冲重复频率PRF得到Δt＝1/PRF；

t为对应的成像时刻；

x为方位向坐标；

t₀为参考时刻；

a_i,b_i,c_i(i＝0,1,2)为多项式拟合系数，由SAR遥感影像的头文件中提供的若干个卫星状态向量拟合得到。

目前，星载SAR遥感影像一般都进行了多普勒归零处理，且地物在地心旋转坐标系中速度为零，为此多普勒条件方程(2)可简化为：

(X-X_s)V_SX+(Y-Y_s)V_SY+(Z-Z_s)V_SZ＝0(4)

联立公式(1)和公式(4)即为距离多普勒数学模型。

本发明实施例中，首先从SAR遥感影像头文件中读取所需的模型参数，然后将时间t和列号y作为未知数，对距离多普勒数学模型线性化，建立误差方程的矩阵形式为：

V＝AX-L,P(5)

式中：

A = [\begin{matrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{matrix}],

为未知数t、y改正数的系数；

X为未知数t、y改正数的向量；

P为权矩阵。

方式二：由某像点像素坐标(x,y)和人工设定的投影面高程h，利用SAR严格几何模型并引入地球椭球模型，计算该像点对应物方点的地心旋转坐标系坐标(X,Y,Z)，通过坐标转换可得到对应的大地经纬度坐标(lat,lon,h)。

所引入的地球椭球模型表示如下：

\frac{(X^{2} + Y^{2})}{{(r_{e} + h)}^{2}} + \frac{Z^{2}}{r_{p}^{2}} = 1 - - - (6)

式中：

r_e,r_p分别为地球椭球的长、短半轴；

h为人工设定的投影面高程。

实施例的计算过程中，首先从SAR遥感影像头文件中读取所需的模型参数，然后以物方点的地心旋转坐标系坐标(X,Y,Z)为未知数，对公式(1)、式(4)和式(6)线性化，建立误差方程的矩阵形式为：

V＝BX-L,P(7)

式中：

B = [\begin{matrix} B_{11} & B_{12} & B_{13} \\ B_{21} & B_{22} & B_{23} \\ B_{31} & B_{32} & B_{33} \end{matrix}]

为(X,Y,Z)的改正数对应的系数矩阵；

X为(X,Y,Z)的改正数向量；

L＝-[F₁F₂F₃]^T为常数向量；

P为权矩阵。

根据误差方程(7)，按照最小二乘法迭代求解可解算出(X,Y,Z)，经坐标转换即可得到大地经纬度坐标(lat,lon,h)。

1.2)根据方式一以及方式二所得结果，对物方点在SAR遥感影像上对应像点的像素坐标进行从侧视斜距投影到光学透视投影的转换，构建SAR遥感影像在物方均匀分布的虚拟控制点格网。

具体包括以下步骤：

影像覆盖范围的确定：根据影像四个角点的像素坐标，按方式二所述方法，求出对应物方点的大地经纬度坐标，从而确定影像覆盖范围。

立体格网的建立：读取影像覆盖范围内的免费全球DEM数据，确定出的最大和最小高程值，并根据给定的高程分层数建立三维立体格网。

虚拟控制点格网的建立：对三维立体格网中各格网点，由其大地经纬度坐标，根据方式一中所述方法计算各格网点在SAR遥感影像上对应的像点坐标，并进一步对像点坐标进行从侧视斜距投影到光学透视投影的转换，从而建立物方均匀分布的虚拟控制点格网。

下面将以图1为例对本发明的步骤1)进行详细说明：

图1给出了SAR遥感影像像点从侧视斜距投影转换为透视中心投影的示意图。O为SAR天线所在的位置，并假设天线距离整个测区平均高程面的高度为H，其在地面上的投影为E点。A点为SAR遥感影像的近地点所对应的地面点，而B点为远地点所对应的地面点。OO'为SAR遥感影像成像面，DD'为等效焦距为f时的光学影像成像面。设近地点斜距OA的距离为R₀，天线航高为H，故可以计算近地点对应侧视角：

θ₀＝arccos(H/R₀)(8)

远地点的斜距OB＝R₀+M_xW，其中M_x为SAR遥感影像的斜距分辨率，W为SAR遥感影像的宽度，则有：

A B = \sqrt{{OB}^{2} - H^{2}} - \sqrt{{OA}^{2} - H^{2}} - - - (9)

令近地点与远地点在光学影像成像面成像面上的投影长度ab为W，则等效焦距f为：

f = \frac{H \times W}{A B} - - - (10)

设物方点P高程为h_p，其在SAR遥感影像对应像点坐标为(x,y)，像点坐标(x,y)转换为透视中心投影坐标为(x',y')，则有：

\{\begin{matrix} x^{'} = \frac{f}{H - h_{p}} \sqrt{{(R_{0} + x \times M_{x})}^{2} - {(H - h_{p})}^{2}} - f t a n (θ_{0}) \\ y^{'} = y \end{matrix} - - - (11)

1.3)根据步骤1.2)建立的虚拟控制点格网求解SAR遥感影像有理函数多项式模型参数。

RFM模型的具体定义关系式如下：

\{\begin{matrix} X = \frac{{Num}_{L} (P, L, H)}{{Den}_{L} (P, L, H)} \\ Y = \frac{{Num}_{s} (P, L, H)}{{Den}_{s} (P, L, H)} \end{matrix} - - - (12)

式中：

Num_L(P,L,H)、Den_L(P,L,H)、Num_s(P,L,H)以及Den_s(P,L,H)为三次多项式，三次多项式的系数即称为RPC参数；

P = \frac{L a t - L A T_O F F}{L A T_S C A L E}

L = \frac{L o n - L O N G_O F F}{L O N G_S C A L E}

H = \frac{H e i g h t - H E I G H T_O F F}{H E I G H T_S C A L E}

X = \frac{S a m p l e - S A M P_O F F}{S A M P_S C A L E}

Y = \frac{L i n e - L I N E_O F F}{L I N E_S C A L E} - - - (13)

式中：

LAT_OFF、LAT_SCALE、LONG_OFF、LONG_SCALE、HEIGHT_OFF和HEIGHT_SCALE为地面点坐标的正则化参数；

SAMP_OFF、SAMP_SCALE、LINE_OFF和LINE_SCALE为影像坐标的正则化参数。Lat表示纬度，Lon表示经度，Height表示高程，此三项即可代表地面点的空间坐标。

Sample代表影像列坐标，Line代表影像行坐标。

本实施例采用最小二乘方法进行求解，以RPC参数为未知数，对有理函数多项式模型进行线性化，列立误差方程：

V＝BX-L,P(14)

式中

B = [\begin{matrix} \frac{\partial F_{x}}{\partial a_{i}} & \frac{\partial F_{x}}{\partial b_{i}} & \frac{\partial F_{x}}{\partial c_{i}} & \frac{\partial F_{x}}{\partial d_{i}} \\ \frac{\partial F_{y}}{\partial a_{i}} & \frac{\partial F_{y}}{\partial b_{i}} & \frac{\partial F_{y}}{\partial c_{i}} & \frac{\partial F_{y}}{\partial d_{i}} \end{matrix}], (i = 0, ..., 19)

L = [\begin{matrix} - F_{1}^{0} \\ - F_{2}^{0} \end{matrix}]

x＝[a_ib_ic_id_i]^T

P为权矩阵(15)

采用谱修正迭代算法，可计算出所述的RPC参数。

步骤2)基于类光学像方改正的异源光学和SAR遥感影像区域网平差

待人工量测所需要的异源光学和SAR遥感影像之间的连接点和地面控制点对应像点的像素坐标后，可以基于类光学像方改正的异源光学和SAR遥感影像区域网平差，实现其精确联合定位。

对于光学遥感影像，直接采用像方仿射变换模型，即：

\{\begin{matrix} F_{x} = {px}_{0} + {px}_{1} \cdot S a m p l e + {px}_{2} \cdot L i n e + S a m p l e - x = 0 \\ F_{y} = {py}_{0} + {py}_{1} \cdot S a m p l e + {py}_{2} \cdot L i n e + L i n e - y = 0 \end{matrix} - - - (16)

式中，(x,y)为物方点在光学遥感影像上对应像点的像素坐标，(Sample,Line)为物方点利用所在光学遥感影像的RPC参数投影到光学遥感影像上的像素坐标，(px_i,py_i)(i＝0,1,2)是待求光学遥感影像对应的仿射变换系数，px₁、py₂的初始值为1，其他系数初始值均为0。

对于SAR遥感影像，因将其侧视斜距投影的像素坐标进行了透视中心投影的转换，可采用光学遥感影像的像方仿射变换模型，即：

\{\begin{matrix} F_{x}^{'} = {px}_{0}^{'} + {px}_{1}^{'} \cdot s a m p l e + {px}_{2}^{'} \cdot l i n e + s a m p l e - x^{'} = 0 \\ F_{y}^{'} = {py}_{0}^{'} + {py}_{1}^{'} \cdot s a m p l e + {py}_{2}^{'} \cdot l i n e + l i n e - y^{'} = 0 \end{matrix} - - - (17)

式中，(x′,y′)为对物方点在SAR遥感影像上对应像点的像素坐标进行透视中心投影转换后的像素坐标，(sample,line)为利用所在SAR遥感影像的RPC参数计算的地面点投影到等效光学影像面上的像素坐标，(px′_i,py′_i)(i＝0,1,2)是SAR遥感影像对应的仿射变换系数。

本实施例中，首先人工量测所需要的异源光学和SAR遥感影像之间的连接点和地面控制点对应像点的像素坐标，然后利用光学、SAR遥感影像的RPC参数，对每个连接点和控制点，以其对应的物方坐标，以及各图像对应的仿射变换系数为未知数列立误差方程。

对于光学遥感影像上的每个连接点和地面控制点的像素坐标，按公式(16)进行线性化，列立误差方程：

V₁＝A₁t+B₁X-L₁,P₁(18)

对于SAR遥感影像上的每个连接点和地面控制点的像素坐标，按公式(17)进行线性化，列立误差方程：

V₂＝A₂t+B₂X-L₂,P₂(19)

对于地面控制点，将其物方坐标视为带权观测值进行处理，还需列立如下误差方程：

V₃＝EX-L₃,P₃(20)

联立式(18)～式(20)，即为类光学像方改正的异源光学和SAR遥感影像区域网平差模型。式中，t为各影像对应的仿射变换系数改正数向量，X为连接点、地面控制点对应的物方坐标改正数向量，A_i和B_i(i＝1,2)分别为未知数t和X对应的系数矩阵，L_i(i＝1,2,3)为对应误差方程的常数向量，P_i为(i＝1,2,3)为相应的权矩阵。利用最小二乘原理，可解算出未知参数的改正数。将与上次迭代的参数值相加，并利用更新后的参数值重新进行平差计算，直至达到收敛条件，即可完成异源光学和SAR遥感影像的联合定位。

Claims

1.一种基于类光学像方改正的异源光学和SAR遥感影像联合定位方法，其特征在于：所述基于类光学像方改正的异源光学和SAR遥感影像联合定位方法包括以下步骤：

1)解算异源光学和SAR遥感影像的RPC参数；

2)利用影像匹配或者人工量测的方法，获取光学以及SAR遥感影像间连接点和地面控制点对应像点的像素坐标；

2.根据权利要求1所述的基于类光学像方改正的异源光学和SAR遥感影像联合定位方法，其特征在于：所述步骤1)的具体实现方式是：

1.1)构建SAR遥感影像在物方均匀分布的虚拟控制点格网；

3.根据权利要求2所述的基于类光学像方改正的异源光学和SAR遥感影像联合定位方法，其特征在于：所述步骤1.1)的具体实现方式是：

4.根据权利要求3所述的基于类光学像方改正的异源光学和SAR遥感影像联合定位方法，其特征在于：所述步骤1.1.1)中求解对应物方点的大地经纬度坐标的具体实现方式是：

R_{S l a n t} = | R_{s} - R_{p} | = \sqrt{{(X - X_{S})}^{2} + {(Y - Y_{S})}^{2} + {(Z - Z_{S})}^{2}}

(X-X_s)V_SX+(Y-Y_s)V_SY+(Z-Z_s)V_SZ＝0

其中：

R_Slant为卫星天线中心S到地面点P的斜距；

所述地球椭球模型的函数表达式是：

\frac{(X^{2} + Y^{2})}{{(r_{e} + h)}^{2}} + \frac{Z^{2}}{r_{p}^{2}} = 1

其中：

r_e,r_p分别为地球椭球的长、短半轴；

h为人工设定的投影面高程；

V＝BX-L,P

式中：

B = [\begin{matrix} B_{11} & B_{12} & B_{13} \\ B_{21} & B_{22} & B_{23} \\ B_{31} & B_{32} & B_{33} \end{matrix}]

为(X,Y,Z)的改正数对应的系数矩阵；

X为(X,Y,Z)的改正数向量；

L＝-[F₁F₂F₃]^T为常数向量；

P为权矩阵；

5.根据权利要求4所述的基于类光学像方改正的异源光学和SAR遥感影像联合定位方法，其特征在于：所述步骤1.1.3)中由其大地经纬度坐标计算各格网点在SAR遥感影像上对应的像点坐标的具体实现方式是：

R_{S l a n t} = | R_{s} - R_{p} | = \sqrt{{(X - X_{S})}^{2} + {(Y - Y_{S})}^{2} + {(Z - Z_{S})}^{2}} - - - (1)

f_{D} = - \frac{2 (R_{s} - R_{p}) \cdot (V_{s} - V_{p})}{λ | R_{s} - R_{p} |} - - - (2)

其中：

R_Slant为卫星天线中心S到地面点P的斜距；

f_D为多普勒频率；

λ为卫星天线所接收的雷达的波长；

其中，斜距R_Slant为距离向坐标y的函数；

对于斜距影像，R_Slant＝R₀+yM_y，R₀为第一斜距，M_y为斜距影像的距离向分辨率，参数R₀和M_y可以从SAR遥感影像头文件中直接读取，y为像点的像素坐标(x,y)在斜距向的值；

\{\begin{matrix} X_{s} = a_{0} + a_{1} t + a_{2} t^{2} \\ Y_{s} = b_{0} + b_{1} t + b_{2} t^{2} \\ Z_{s} = c_{0} + c_{1} t + c_{2} t^{2} \\ V_{s x} = a_{1} + 2 a_{2} t \\ V_{s y} = b_{1} + 2 b_{2} t \\ V_{s z} = c_{1} + 2 c_{2} t \end{matrix} - - - (3)

t＝t₀+(x-1)Δt

其中：

t为对应的成像时刻；

x为方位向坐标；

t₀为参考时刻；

(X-X_s)V_SX+(Y-Y_s)V_SY+(Z-Z_s)V_SZ＝0(4)

联立公式(1)和公式(4)即为距离多普勒数学模型；

V＝AX-L,P(5)

其中：

A = [\begin{matrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{matrix}],

为未知数t、y改正数的系数；

X为未知数t、y改正数的向量；

P为权矩阵；

6.根据权利要求5所述的基于类光学像方改正的异源光学和SAR遥感影像联合定位方法，其特征在于：所述步骤1.1.3)中对像点坐标进行从侧视斜距投影到光学透视投影的转换的具体实现方式是：

则近地点A对应侧视角为：

θ₀＝arccos(H/R₀)

其中：

R₀为近地点斜距OA的距离；

H为卫星天线航高；

远地点B的斜距为：

OB＝R₀+M_xW

其中：

M_x为SAR遥感影像的斜距分辨率；

W为SAR遥感影像的宽度；

根据公式(7)，则有近地点A与远地点B之间的距离AB是：

A B = \sqrt{{OB}^{2} - H^{2}} - \sqrt{{OA}^{2} - H^{2}}

f = \frac{H \times W}{A B}

\{\begin{matrix} x^{'} = \frac{f}{H - h_{p}} \sqrt{{(R_{0} + x \times M_{x})}^{2} - {(H - h_{p})}^{2}} - f t a n (θ_{0}) \\ y^{'} = y \end{matrix} .

7.根据权利要求2或3或4或5或6所述的基于类光学像方改正的异源光学和SAR遥感影像联合定位方法，其特征在于：所述步骤1.2)中有理函数多项式模型的具体定义是：

\{\begin{matrix} X = \frac{{Num}_{L} (P, L, H)}{{Den}_{L} (P, L, H)} \\ Y = \frac{{Num}_{s} (P, L, H)}{{Den}_{s} (P, L, H)} \end{matrix}

其中：

P = \frac{L a t - L A T_O F F}{L A T_S C A L E}

L = \frac{L o n - L O N G_O F F}{L O N G_S C A L E}

H = \frac{H e i g h t - H E I G H T_O F F}{H E I G H T_S C A L E}

X = \frac{S a m p l e - S A M P_O F F}{S A M P_S C A L E}

Y = \frac{L i n e - L I N E_O F F}{L I N E_S C A L E}

其中：

Sample代表影像列坐标；Line代表影像行坐标。

8.根据权利要求7所述的基于类光学像方改正的异源光学和SAR遥感影像联合定位方法，其特征在于：所述步骤1.3)的具体实现方式是：

V＝BX-L,P

其中：

B = [\begin{matrix} \frac{\partial F_{x}}{\partial a_{i}} & \frac{\partial F_{x}}{\partial b_{i}} & \frac{\partial F_{x}}{\partial c_{i}} & \frac{\partial F_{x}}{\partial d_{i}} \\ \frac{\partial F_{y}}{\partial a_{i}} & \frac{\partial F_{y}}{\partial b_{i}} & \frac{\partial F_{y}}{\partial c_{i}} & \frac{\partial F_{y}}{\partial d_{i}} \end{matrix}] (i = 0, ..., 19)

L = [\begin{matrix} - F_{1}^{0} \\ - F_{2}^{0} \end{matrix}]

x＝[a_ib_ic_id_i]^T

P为权矩阵

采用谱修正迭代算法，计算得到SAR遥感影像的RPC参数。

9.根据权利要求8所述的基于类光学像方改正的异源光学和SAR遥感影像联合定位方法，其特征在于：所述步骤3)的具体实现方式是：

3.1)计算和给定各未知参数的初始值；

3.3)利用最小二乘原理，可解算出未知参数的改正数；3.4)将步骤3.3)解算得到的未知参数的改正数与上次迭代的未知参数值相加，得到本次迭代更新后的未知参数值；当第一次迭代时，未知参数值为初始值；

3.5)利用更新后的未知参数值按步骤3.2)～3.4)重新进行计算，直至达到收敛条件，即完成异源光学和SAR遥感影像的联合定位。

10.根据权利要求9所述的基于类光学像方改正的异源光学和SAR遥感影像联合定位方法，其特征在于：所述步骤3.2)的具体实现方式是：

对于光学遥感影像，直接采用像方仿射变换模型，即：

\{\begin{matrix} F_{x} = {px}_{0} + {px}_{1} \cdot S a m p l e + {px}_{2} \cdot L i n e + S a m p l e - x = 0 \\ F_{y} = {py}_{0} + {py}_{1} \cdot S a m p l e + {py}_{2} \cdot L i n e + L i n e - y = 0 \end{matrix} - - - (6)

其中：

(x,y)为物方点在光学遥感影像上对应像点的像素坐标；

\{\begin{matrix} F_{x}^{'} = {px}_{0}^{'} + {px}_{1}^{'} \cdot s a m p l e + {px}_{2}^{'} \cdot l i n e + s a m p l e - x^{'} = 0 \\ F_{y}^{'} = {py}_{0}^{'} + {py}_{1}^{'} \cdot s a m p l e + {py}_{2}^{'} \cdot l i n e + l i n e - y^{'} = 0 \end{matrix} - - - (7)

其中：

(px′_i,py′_i)(i＝0,1,2)是SAR遥感影像对应的仿射变换系数；

V₁＝A₁t+B₁X-L₁,P₁(8)

V₂＝A₂t+B₂X-L₂,P₂(9)

V₃＝EX-L₃,P₃(10)

其中：

t为各影像对应的仿射变换系数改正数向量；

X为连接点以及地面控制点对应的物方坐标改正数向量；

A_i和B_i(i＝1,2)分别为未知数t和X对应的系数矩阵；

L_i(i＝1,2,3)为对应误差方程的常数向量；

P_i为(i＝1,2,3)为相应的权矩阵。