CN107341778B - 基于卫星控制点库和dem的sar影像正射纠正方法 - Google Patents

Abstract

本发明公开了一种基于卫星控制点库和DEM的SAR影像正射纠正方法，其利用外部已有DEM数据和资源三号控制点数据，通过SAR影像模拟技术和影像精确匹配，得到资源三号控制点在SAR影像中的像点坐标，随后通过高精度的资源三号控制点参与SAR正射纠正，得到SAR正射纠正后的DOM影像。本发明可以自动地利用资源三号控制点库参与SAR影像正射纠正，减少人工干预带来的效率低下或人工选点带来的刺点误差；同时能有效用于SAR的正射纠正，降低获取控制点的成本，SAR影像纠正效果好。

Description

基于卫星控制点库和DEM的SAR影像正射纠正方法

技术领域

本发明涉及雷达图像几何处理技术领域，特别涉及一种基于卫星控制点库和DEM的SAR影像正射纠正方法。

背景技术

资源三号测绘卫星于2012年1月9号成功发射，是我国第一颗民用三线阵立体测图卫星，主要用于1:50000立体测图及更大比例尺基础地理产品的生产和更新以及开展国土资源调查与监测。为满足国产资源三号测绘卫星影像数据处理、应急测绘保障等规模化应用，首次建立了覆盖全国范围的满足1:50000精度的影像控制点数据库，包含361025个高精度控制点数据，它们已全面应用于我国1:50000比例尺光学数字正射影像生产和1:25000比例尺的全国正射影像更新。由于光学影像与SAR影像的成像几何不同，散射机理不同，导致纹理和强度信息相差极大，因此使直接将资源三号卫星影像控制点转刺到SAR影像中，会带来数个像素的转刺误差。而如果直接采用严密成像模型直接定位，则会由于成像几何误差的影响带来更大的定位误差，直接应用于SAR影像正射纠正很难满足我国测绘比例尺的精度要求。

一般在无实测地面控制点数据时，SAR影像的正射纠正是利用外部参考DEM数据模拟SAR影像进行几何校正，其原理是：根据真实SAR影像的轨道状态矢量和成像参数以及外部参考DEM数据，通过SAR模拟技术模拟SAR影像的成像几何，得到模拟SAR影像。由于模拟SAR影像与真实SAR影像使用相同的轨道状态矢量和成像参数，它们之间的几何差异可以通过影像匹配确定，进而得到模拟SAR影像与真实SAR影像一一对应的像点映射关系。而模拟SAR影像与参考DEM之间同样存在一一映射的像点关系，可以在DEM上选取一部分点作为外部参考控制点数据，用于SAR影像的正射纠正处理。一般工程应用均使用的是WGS84参考椭球系统，WGS84椭球高与EGM96正高系统之间需要通过大地水准面差距进行转换。此外，国内学者对SRTM DEM的高程精度进行验证，利用三峡、华北、云贵地区总计206个GPS观测数据进行DEM高程精度验证，结果表明高程均方根误差为10.31m。资源三号卫星影像控制点数据的平面和高程精度优于厘米级，其高程精度要高于从SRTM DEM上选取的控制点数据，即使用资源三号卫星控制点数据可以进一步提高SAR影像正射纠正精度。现有技术中一般对SAR影像采用严密的成像模型进行定位处理和模型参数优化，采用的模型参数较多且计算复杂。

因此，需要提出一种能够克服上述缺点的使用资源三号卫星控制点数据和DEM辅助下的用于SAR影像正射纠正的方法。

发明内容

本发明提供了一种基于卫星控制点库和DEM的SAR影像正射纠正方法，其特征在于，所述方法具体包括以下步骤：

步骤S1，采集原始数据，所述原始数据优选的包括资源三号卫星影像控制点数据、外部参考DEM数据、原始SAR影像以及原始SAR影像成像参数；

步骤S2，对所述原始SAR影像进行多视处理，得到多视SAR影像；

步骤S3，基于所述外部参考DEM数据和多视SAR影像，计算外接矩形数字高程模型(RefDEM)数据、DEM模拟SAR影像、DEM模拟SAR影像和多视SAR影像雷达坐标空间一一映射的坐标关系；

步骤S4，根据步骤S3.3获取的DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系和步骤S3.4获取的DEM模拟SAR影像，计算得到与多视SAR影像雷达坐标空间一致的SAR模拟影像；

步骤S5，对所述SAR模拟影像和所述多视SAR影像进行灰度匹配，计算得到新的DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系；

步骤S6，基于所述新的DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系，计算资源三号卫星影像控制点数据在多视SAR影像中的像点坐标集；

步骤S7，将所述资源三号卫星影像控制点数据在多视SAR影像中的像点坐标集转换为所述资源三号卫星影像控制点数据在原始SAR影像中的坐标集；

步骤S8，利用资源三号卫星影像控制点数据在原始SAR影像中的坐标集，进行SAR影像正射纠正处理，计算得到距离向和方位向多项式的改正系数；

步骤S9，根据所述距离向和方位向多项式的改正系数，对所述原始SAR影像进行SAR正射纠正，得到SAR正射纠正后的DOM影像。

优选的，步骤2具体包括：根据下式对所述原始SAR影像进行多视处理，得到多视SAR影像：

其中，Pix(i,j)表示多视SAR影像中的像元(i,j)对应的像元值；nAzML表示原始SAR影像的方位向多视数，nRgML表示原始SAR影像的距离向多视数；Pix(n,m)表示原始SAR影像中的像元(n,m)对应的像元值，n表示原始SAR影像中的像元(n,m)的方位向坐标，取值范围为n＝1,...,nAzML，m表示原始SAR影像中的像元(n,m)的距离向坐标，取值范围为m＝1,...,nRgML。

优选的，所述步骤3优选的包括以下子步骤：

步骤S3.1，基于所述外部参考DEM数据和多视SAR影像，计算外接矩形数字高程模型(RefDEM)数据，其中，所述外接矩形数字高程模型(RefDEM)数据包括所述外接矩形数字高程模型(RfDEM)在外部参考DEM数据中的宽、高、起始点地理坐标、纬度方向分辨率、经度方向分辨率和高程值；

步骤S3.2，根据步骤S3.1获取的外接矩形数字高程模型(RefDEM)的宽nWidth_DEM和高nHeight_DEM，创建两个大小均为(nHeight_DEM,nWidth_DEM)的实数空矩阵，用于分别存储外接矩形数字高程模型(RefDEM)数据和DEM模拟SAR影像数据；创建一个大小为(nHeight_DEM,nWidth_DEM)的复数空矩阵，用于存储DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系数据；

步骤S3.3，基于所述外部参考DEM数据和多视SAR影像，根据SAR严密成像模型(R-D模型)，计算外接矩形数字高程模型(RefDEM)数据与多视SAR影像雷达坐标空间一一映射的坐标关系；

步骤S3.4，根据SAR影像模拟公式，根据外接矩形数字高程模型(RefDEM)数据模拟SAR像元(i,j)处的影像值Pix_Sim(i,j)，最终所有模拟SAR像元值Pix_Sim(i,j)的集合构成DEM模拟SAR影像。

优选的，所述步骤S3.1可以包括：

步骤S3.1.1，利用外部参考DEM(ExtDEM)数据和多视SAR影像的轨道状态、成像参数信息，根据SAR严密成像模型(R-D模型)计算得到多视SAR影像的四个角点坐标在外部参考DEM(ExtDEM)数据中的像素坐标；

步骤S3.1.2，根据步骤S3.1.1获取的多视SAR影像的四个角点在外部参考DEM(ExtDEM)数据中的像素坐标(i_DEM,j_DEM)，计算外接矩形数字高程模型(RefDEM)的四个角点坐标在外部参考DEM(ExtDEM)数据中的坐标，以及外接矩形数字高程模型(RefDEM)在外部参考DEM数据中的宽nWidth_DEM和高nHeight_DEM；

步骤S3.1.3，根据步骤S3.1.2得到的外接矩形数字高程模型(RefDEM)的左上角C_lt(i,j)角点坐标，结合外部参考DEM(ExtDEM)数据提供的起点地理坐标(B_ext,0,L_ext,0)、纬度方向的分辨率B_ext,res、经度方向的分辨率L_ext,res，根据下式计算外接矩形数字高程模型(RefDEM)的起始点地理坐标(B₁,L₁)；

B₁＝B_ext,0+B_ext,res*i

L₁＝L_ext,0+L_ext,res*j

其中，(B₁,L₁)为外接矩形数字高程模型(RefDEM)的起始点地理坐标；(B_ext,0,L_ext,0)为外部参考DEM(ExtDEM)数据提供的起点地理坐标，B_ext,res为外部参考DEM(ExtDEM)数据的纬度方向的分辨率，L_ext,res为外部参考DEM(ExtDEM)数据的经度方向的分辨率；i为外接矩形数字高程模型(RefDEM)在外部参考DEM(ExtDEM)数据中的左上角的行方向坐标，j为外接矩形数字高程模型(RefDEM)在外部参考DEM(ExtDEM)数据中的左上角的列方向坐标；

步骤S3.1.4，根据步骤S3.1.3获取的外接矩形数字高程模型(RefDEM)的起始点地理坐标(B₁,L₁)，计算外接矩形数字高程模型(RefDEM)的高程值。

优选的，所述步骤S3.1可以包括以下子步骤：

步骤S3.1.1.1，根据多视SAR影像的角点的方位向坐标i_AZ,计算卫星方位向成像时间t_T，并根据卫星方位向成像时间tT计算卫星位置矢量(X_s,Y_s,Z_s)和速度矢量

步骤S3.1.1.2，对多视SAR影像的中心点经纬度

进行坐标转换，将其从大地坐标系转换到地心空间直角坐标系，作为待求点的地心坐标初值(X_T ⁰,Y_T ⁰,Z_T ⁰)；

步骤S3.1.1.3，基于步骤S3.1.1.1求得的卫星位置矢量(X_s,Y_s,Z_s)和速度矢量

和步骤S3.1.1.2求得的待求点的地心坐标初值(X_T ⁰,Y_T ⁰,Z_T ⁰)，根据SAR严密成像模型，解算所述待求点的改正数；

步骤S3.1.1.4，根据下式，利用待求点的改正数更新待求点的地心坐标：

X_T ¹＝X_T ⁰+ΔX_T，Y_T ¹＝Y_T ⁰+ΔY_T，Z_T ¹＝Z_T ⁰+ΔZ_T；

其中，

是更新后的待求点的地心坐标，

是待求点的地心坐标初值，ΔX_T，ΔY_T，ΔZ_T是待求点的改正坐标；

步骤S3.1.1.5，利用更新后的待求点的地心坐标(X_T ¹,Y_T ¹,Z_T ¹)进行迭代解算，重复S3.1.1.3-S3.1.1.4，直至改正数(ΔX_T,ΔY_T,ΔZ_T)的绝对值都小于设定阈值ε迭代结束，确定待求点的地心坐标(X_T,Y_T,Z_T)；

步骤S3.1.1.6，对通过步骤S3.1.1.5获取的待求点的地心坐标(X_T,Y_T,Z_T)进行坐标转换，将其从地心空间直角坐标系转换到大地坐标系，得到并输出待求点的大地坐标(B_T,L_T,H_T)；

步骤S3.1.1.7，根据下述公式，计算得到多视SAR影像的角点坐标在外部参考DEM(ExtDEM)数据中的像素坐标(i_DEM,j_DEM)：

i_DEM＝(B_ext,0-B_i)/B_ext,res

j_DEM＝(L_ext,0-L_i)/L_ext,res

其中，(i_DEM,j_DEM)为多视SAR影像的角点坐标在外部参考DEM(ExtDEM)数据中的像素坐标；(B_ext,0,L_ext,0)为外部参考DEM(ExtDEM)数据的起点地理坐标；B_ext,res为纬度方向分辨率，L_ext,res为经度方向分辨率；B_i为多视SAR影像的角点坐标对应的大地坐标纬度值，L_i为多视SAR影像角点坐标对应的大地坐标经度值。

优选的，步骤S3.1.1.3优选的包括：

所述SAR严密成像模型为：

其中，Vx_S·(X_S-X_P)+Vy_S·(Y_S-Y_P)+Vz_S·(Z_S-Z_P)＝-f_d·λ·(R₀+ΔR·x)/2为多普勒方程，

为距离方程，(X_S,Y_S,Z_S)为卫星位置矢量，(X_P,Y_P,Z_P)为角点的地心空间直角坐标，(V_XS,V_YS,V_ZS)为卫星速度矢量，x为角点的距离向坐标，R₀为近地点斜距，ΔR为距离向分辨率，f_d为多普勒中心频率，λ为雷达波长；

根据上述公式(2)，构建观测方程，计算多视SAR影像的角点的误差方程系数，形成误差方程系数矩阵，再根据误差方程系数矩阵建立法方程：B^TBX-B^TL＝0，根据下式解算待求点的改正数：

X＝(B^TB)^-1B^TL

其中，X＝(ΔX_T,ΔY_T,ΔZ_T)，ΔX_T，ΔY_T，ΔZ_T是待求点的改正数坐标；B为误差方程系数矩阵；L为常数项矩阵；；

其中，，在步骤S3.1.4中，根据下述公式计算外接矩形数字高程模型(RefDEM)的高程值：

B_RefDEM＝B₁+B_ref,res*m,m＝0,1,…,nHeight_DEM-1

L_RefDEM＝L₁+L_ref,res*n,n＝0,1,…,n Width_DEM-1

i_DEM＝(B_ext,0-B_RefDEM)/B_ext,res

j_DEM＝(L_ext,0-L_RefDEM)/L_ext,res

H_DEM＝Pix_DEM(i_DEM,j_DEM)

其中，(B_RefDEM,L_RefDEM)为外接矩形数字高程模型(RefDEM)的某个点的经纬度坐标，B_RefDEM为外接矩形数字高程模型(RefDEM)的某个点的纬度坐标，L_RefDEM为外接矩形数字高程模型(RefDEM)的某个点的经度坐标；(B₁,L₁)为外接矩形数字高程模型(RefDEM)的起始点地理坐标，B_ref,res为外接矩形数字高程模型(RefDEM)的纬度方向的分辨率，L_ref,res为外接矩形数字高程模型(RefDEM)的经度方向的分辨率，nWidth_DEM和nHeight_DEM分别为外接矩形数字高程模型(RefDEM)的宽和高，m为高度方向的坐标索引，其中m＝0,1,…,nHeight_DEM-1；n为宽度方向的坐标索引，其中n＝0,1,…,nWidth_DEM-1；(B_ext,0,L_ext,0)为外部参考DEM(ExtDEM)数据的起点地理坐标，B_ext,res为外部参考DEM(ExtDEM)在纬度方向的分辨率，L_ext,res为外部参考DEM(ExtDEM)在经度方向的分辨率；(i_DEM,j_DEM)为外接矩形数字高程模型(RefDEM)的某个点的经纬度坐标在外部参考DEM(ExtDEM)数据中的行列坐标；Pix_DEM(i_DEM,j_DEM)表示行列坐标为(i_DEM,j_DEM)的点在外部参考DEM(ExtDEM)数据中的高程值，H_DEM表示外接矩形数字高程模型(RefDEM)的高程值。

优选的，步骤S3.3优选的包括以下子步骤：

步骤S3.3.1，基于公式(2)中的多普勒方程迭代计算外接矩形数字高程模型(RefDEM)数据中的每个像素点的初始方位向坐标y；

S3.3.2，根据公式(2)的距离方程计算外接矩形数字高程模型(RefDEM)数据中的每个像素点的初始距离向坐标x；

S3.3.3，组合步骤S3.3.1和步骤S3.3.2得到的外接矩形数字高程模型(RefDEM)数据中的每个像素点的初始方位向坐标y和初始距离向坐标x，得到外接矩形数字高程模型(RefDEM)数据中的每个像素点在多视SAR影像雷达坐标空间中的像点坐标(x,y)，所述像点坐标(x,y)为外接矩形数字高程模型(RefDEM)数据与多视SAR影像雷达坐标空间一一映射的坐标关系。

优选的，

步骤S4具体包括以下三个步骤：

步骤S4.1，获得DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系(即坐标映射数据1)的每一个像元pix(i,j)的复数像元值(i_Value,j_Value)，如满足如下条件式，则保留该像元pix(i,j)，参与下一步计算：

其中，i_value为某一个像元pix(i,j)对应的复数像元值的高度方向坐标，j_value为某一个像元pix(i,j)对应的复数像元值的宽度方向坐标，nMLHeight为多视SAR影像的高，nMLWidth为多视SAR影像的宽；

步骤S4.2，根据步骤S4.1获取的像元pix(i,j)和对应的复数像元值(i_Value,j_Value)，在DEM模拟SAR影像中取(i,j)像元值填充到复数像元值(i_Value,j_Value)对应的SAR模拟影像坐标中；

步骤S4.3，对DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系(即坐标映射数据1)的每一个像元均采用步骤S4.1和步骤S4.2，得到SAR模拟影像；

步骤S5具体包括以下三个子步骤：

步骤S5.1，采用匹配方法，在SAR模拟影像上规则选取一定数量的点，对搜索窗口数为nWin范围内的像素进行匹配，得到窗口内每个像元的匹配测度ρ，取最大的匹配测度ρ，以及对应的方位向偏移量记为off_Az，距离向偏移量记为off_Rg。匹配方法采用常用的NCC匹配方法(归一化互相关系数)，其公式如下：

其中，r_i表示多视SAR影像(参考影像)窗口内的像素值；m_i表示SAR模拟影像(待匹配影像)窗口内的像素值；

表示多视SAR影像窗口内的像素均值；

表示模拟SAR影像窗口内的像素均值；σ_r表示多视SAR影像窗口内的像素标准差；σ_m表示模拟SAR影像窗口内的像素标准差；n表示匹配窗口内的像素个数，n＝nWin²，nWin为搜索窗口大小；

步骤S5.2，对步骤S5.1获取的方位向的偏移量off_Az、距离向的偏移量off_Rg，进行多项式拟合，得到方位向和距离向的偏移量拟合系数；多项式拟合公式为：

y(x_Az)＝a₀+a₁x_Az

y(x_Rg)＝b₀+b₁x_Rg

其中，y(x_Az)、y(x_Rg)为方位向和距离向的偏移量矩阵，矩阵大小为k个元素；x_Az、x_Rg为对应的元素索引，取值为x_Az＝1,2,...,k，x_Rg＝1,2,...,k；a₀、a₁为方位向的多项式常数项和一次项系数；b₀、b₁为距离向的多项式常数项和一次项系数；

步骤S5.3，根据步骤S5.2获取的方位向和距离向的偏移量拟合系数，更新坐标映射数据1得到新的坐标映射数据2。更新的公式如下：

对坐标映射数据1的每一个像元(i,j)，得到其复数像元值(i_Value,1,j_Value,1)，根据方位向和距离向的偏移量，有

i_Value,2＝i_Value,1+a₀+a₁·i_Value,1

j_Value,2＝j_Value,1+b₀+b₁·j_Value,1

其中，i_Value,2、j_Value,2为更新后的坐标映射数据2的方位向和距离向像元值；i_Value,1、j_Value,1为坐标映射数据1的方位向和距离向像元值；a₀、a₁为方位向的多项式常数项和一次项系数；b₀、b₁为距离向的多项式常数项和一次项系数。

优选的，步骤S6具体包括以下子步骤：

步骤S6.1，计算ZY3数据在外接矩形数字高程模型(RefDEM)中的坐标(i_DEM,j_DEM)，计算公式如下所示：

i_DEM＝(B_i-B_ref,1)/B_ref,res

j_DEM＝(L_i-L_ref,1)/L_ref,res

其中，(B_ref,1,L_ref,1)为外接矩形数字高程模型(RefDEM)的起始点地理坐标；(B_i,L_i)为资源三号卫星影像控制点数据ZY3的经纬度地理坐标；B_ref,res为纬度方向外接矩形数字高程模型(RefDEM)分辨率，L_ref,res为经度方向外接矩形数字高程模型(RefDEM)分辨率；(i_DEM,j_DEM)为ZY3数据在外接矩形RefDEM中的坐标；

步骤S6.2，计算(i_DEM,j_DEM)在SAR模拟影像中的像素坐标(i_MSAR,j_MSAR)；

步骤S6.3，将所有ZY3数据重复步骤S6.1和S6.2，获得ZY3在SAR模拟影像中的像点坐标(i_MSAR,j_MSAR)，并基于所述像点坐标(i_MSAR,j_MSAR)和ZY3的地理坐标和高程构成ZY3控制点数据在多视SAR影像中的坐标集(B_i,L_i,H_i,i_MSAR,j_MSAR)；其中，(B_i,L_i,H_i)为ZY3的地理坐标和高程；(i_MSAR,j_MSAR)为ZY3在多视SAR影像中的像点坐标。

优选的，步骤S9具体包括以下子步骤：

步骤S9.1，根据原始SAR影像的四个角点像素坐标、原始SAR成像参数和输入的DOM数据的纬度方向分辨率、经度方向分辨率，确定待校正的DOM数据的宽和高，以及DOM起始点地理坐标；

步骤S9.2，根据步骤S9.1的地理参数，利用步骤S3.4构建的RD模型，逐像素对待校正的DOM影像进行计算，得到待校正的DOM像素在原始SAR影像中的像素坐标(x',y')；

步骤S9.3，对步骤S9.2得到的像素坐标(x',y')进行改正，得到改正后的原始SAR像素坐标(x,y)，获取原始SAR像素坐标(x,y)对应的像元值Pix_DOM(x,y)，将所述像元值Pix_DOM(x,y)保存到步骤S9.2所述的待校正的DOM像素坐标pix(x_DOM,y_DOM)中，遍历整个待校正的DOM像素，获得DOM正射纠正影像。

相对于现有技术，本发明的主要优点在于：

1、通过本方法可以自动地利用资源三号控制点库参与SAR影像正射纠正，减少人工干预带来的效率低下或人工选点带来的刺点误差；

2、充分利用了通过光学手段获取的高精度资源三号控制点库，使其能有效用于SAR的正射纠正，降低获取控制点的成本；

3、通过外部DEM模拟SAR影像并与多视SAR影像进行匹配，确定了高精度资源三号控制点库在SAR影像上的坐标集，比直接在外部DEM上选取的点参与SAR影像纠正效果要好。

附图说明

图1为基于卫星控制点库和DEM的SAR正射纠正方法流程图。

图2(a)是外部参考DEM，图2(b)是多视SAR影像。

图3(a)为外接矩形RefDEM；图3(b)为坐标映射数据；图3(c)为DEM模拟SAR影像。

图4(a)为SAR模拟影像；图4(b)为多视SAR影像。

具体实施方式

下面参照附图对本发明进行更全面的描述，其中说明本发明的示例性实施例。

本发明的主要目的是利用外部已有DEM数据和资源三号控制点数据，通过SAR影像模拟技术和影像精确匹配，得到资源三号控制点在SAR影像中的像点坐标，随后通过高精度的资源三号控制点(此时包含像点坐标和经纬度高程坐标)参与SAR正射纠正，得到SAR正射纠正后的DOM影像。

如图1所示，本发明提出的基于卫星控制点库和DEM(Digital Elevation Model,DEM；数字高程模型)的SAR正射纠正方法，利用资源三号卫星影像控制点数据和外部参考DEM数据，首先对原始SAR影像进行多视处理得到多视SAR影像，其次通过SAR模拟技术，模拟一副SAR影像，运用匹配策略对SAR模拟影像和多视SAR影像进行匹配，确定它们之间的几何偏移关系，随后通过几何偏移关系计算得到高精度资源三号卫星影像控制点数据在多视SAR影像中的像点坐标，进而得到资源三号卫星影像控制点数据在原始SAR影像中的像点坐标，最后利用资源三号卫星影像控制点的高程数据替换外部参考DEM对应的高程数据，得到高精度的资源三号控制点(此时包含像点坐标和经纬度高程坐标)，参与SAR影像正射纠正，其可实现资源三号卫星影像控制点数据应用于SAR影像正射纠正处理。该方法具体包括以下步骤：

步骤S1，采集原始数据，所述原始数据包括资源三号卫星影像控制点数据、外部参考DEM数据、原始SAR影像以及原始SAR影像成像参数。优选的，外部参考DEM数据可以是SRTM(Shuttle Radar Topography Mission,SRTM；美国航天飞机雷达测图任务)DEM数据。

具体数据如下：1)资源三号卫星影像控制点数据采用高精度的GPS设备进行外业测量，从资源三号卫星影像控制点数据库中查找并选取实验区控制点数据，如表1所示的资源三号卫星影像控制点数据，以ZY3表示，控制点坐标包括：第一列表示控制点的名称，第二列表示纬度，第三列表示经度，第四列表示大地高程；2)外部参考DEM数字高程模型数据，以ExtDEM表示，它的地理范围是完全包括SAR影像的地理范围；3)原始SAR影像，为了描述方便，将步骤1获取的SAR影像称为“原始SAR影像”；4)原始SAR影像对应的成像参数如表2所示。

控制点名称	纬度/度	经度/度	大地高程/米
				F03	39.80890	116.09411	75.428
F18	39.74416	116.01044	50.419
				F05	39.79243	116.16860	38.754
C24	39.88796	116.01093	220.683
				C22	39.89176	116.02806	230.632
W15	40.12384	116.22403	30.164
				C02	40.12662	116.22181	33.304
W17	40.06218	116.28779	36.213

表1 资源三号卫星影像控制点数据ZY3

表2 原始SAR影像对应的成像参数

由于资源三号卫星影像控制点数据精度更高，因此本发明用资源三号卫星影像控制点数据替换具有相同地理坐标的外部参考DEM高程数据，参与SAR正射纠正处理，能得到更好的SAR影像正射纠正精度，表3为资源三号卫星影像控制点数据与SRTM DEM高程精度对比结果。

控制点名称	纬度/度	经度/度	实测高程/米	SRTM高程/米	误差值
						F03	39.80890	116.09411	75.428	76.727	-1.299
F18	39.74416	116.01044	50.419	53.834	-3.415
						F05	39.79243	116.16860	38.754	37.279	1.475
C24	39.88796	116.01093	220.683	234.741	-14.058
						C22	39.89176	116.02806	230.632	234.092	-3.460
W15	40.12384	116.22403	30.164	37.105	-6.941
						C02	40.12662	116.22181	33.304	34.631	-1.327
W17	40.06218	116.28779	36.213	38.467	-2.254

表3 资源三号卫星影像控制点数据与SRTM DEM高程精度比较

步骤S2，对所述原始SAR影像进行多视处理，得到多视SAR影像。

具体的，在步骤S2中，利用步骤S1获得的原始SAR影像和原始SAR影像成像参数，对原始SAR影像进行多视处理，得到多视SAR影像。

考虑到原始SAR影像受噪声影响比较大，为了增加SAR模拟影像与SAR影像匹配的成功率，首先对原始SAR影像进行多视处理，达到降低噪声对匹配结果的影响。将原始SAR影像的方位向(行方向)多视数设为nAzML，距离向(列方向)多视数设为nRgML。假设一幅原始SAR影像的宽为nWidth，高为nHeight，那么进行多视处理之后，其多视SAR影像的宽nMLWidth为nWidth/nRgML，多视SAR影像的高nMLHeight为nHeight/nAzML。根据下式对所述原始SAR影像进行多视处理，得到多视SAR影像：

其中，Pix(i,j)表示多视SAR影像中的像元(i,j)对应的像元值；nAzML表示原始SAR影像的方位向多视数，nRgML表示原始SAR影像的距离向多视数；Pix(n,m)表示原始SAR影像中的像元(n,m)对应的像元值，n表示原始SAR影像中的像元(n,m)的方位向坐标，取值范围为n＝1,...,nAzML，m表示原始SAR影像中的像元(n,m)的距离向坐标，取值范围为m＝1,...,nRgML。对所有像元均采用公式(1)进行计算，得到多视SAR影像。

优选的，方位向多视数nAzML＝4，距离向多视数nRgML＝4；多视SAR影像的宽nMLWidth＝3461，高nMLHeight＝7176。

步骤S3，基于所述外部参考DEM数据和多视SAR影像，计算外接矩形数字高程模型(RefDEM)数据、DEM模拟SAR影像、DEM模拟SAR影像和多视SAR影像雷达坐标空间一一映射的坐标关系。

具体的，利用步骤S1提供的外部参考DEM数据和步骤S2提供的多视SAR影像的轨道状态和成像参数信息，计算得到外接矩形数字高程模型(RefDEM)数据，通过SAR模拟技术得到DEM模拟SAR影像，以及计算得到DEM模拟SAR影像和多视SAR影像之间一一映射的坐标关系。

SAR模拟技术，首先包括几何模拟，它是以数字高程模型(DEM)、卫星状态矢量参数和成像参数为基本输入信息，基于距离-多普勒方程(如公式2所示)建立几何构像模型，模拟出影像的几何位置；其次，后向散射模拟，本发明基于半经验后向散射模型进行SAR影像的灰度模拟。完成模拟SAR影像之后，采用影像匹配技术，自动完成SAR模拟影像与多视SAR影像上同名目标点的寻找，是SAR正射纠正处理中常用的一种方法。

RefDEM(即，上文所述的外接矩形数字高程模型)数据是模拟SAR影像对应的外接矩形数字高程模型数据；DEM模拟SAR影像是基于外接矩形数字高程模型(RefDEM)数据模拟计算得到的模拟SAR影像，其大小和影像方向与外接矩形数字高程模型(RefDEM)数据相一致。

步骤S3包括：首先根据外部参考DEM数据和多视SAR影像的轨道状态、成像参数信息，结合多视SAR影像的四个角点像元坐标，计算外接矩形数字高程模型(RefDEM)在外部参考DEM数据中的宽nWidth_DEM和高nHeight_DEM，以及相应外接矩形数字高程模型(RefDEM)数据的起始点地理坐标(B₁,L₁)；其次，创建两个大小为(nHeight_DEM,nWidth_DEM)的实数空矩阵，用于存储外接矩形数字高程模型(RefDEM)数据和DEM模拟SAR影像数据，创建一个大小为(nHeight_DEM,nWidth_DEM)的复数空矩阵，用于存储DEM模拟SAR影像和多视SAR影像雷达坐标空间一一映射的坐标关系数据；随后，根据外接矩形数字高程模型(RefDEM)数据的起始点地理坐标(B₁,L₁)，在外部参考DEM数据中取高程值，数据值存储于外接矩形数字高程模型(RefDEM)数据中；根据SAR严密成像模型(R-D模型)，结合外接矩形数字高程模型(RefDEM)数据的起始点地理坐标(B₁,L₁)和相应高程数据，计算其在多视SAR影像雷达坐标空间中的像点坐标，存储于DEM模拟SAR影像和多视SAR影像雷达坐标空间一一映射的坐标关系数据中；随后，根据SAR影像模拟公式(即半经验后向散射模型)，模拟SAR影像灰度值，存储于DEM模拟SAR数据中；最后输出三个数据文件，参与后续影像匹配、资源三号控制点坐标获取和SAR正射纠正处理。

具体来说，步骤S3具体包括以下子步骤：

步骤S3.1，基于所述外部参考DEM数据和多视SAR影像，计算外接矩形数字高程模型(RefDEM)数据，其中，所述外接矩形数字高程模型(RefDEM)数据包括外接矩形数字高程模型(RefDEM)在外部参考DEM数据中的宽、高、起始点地理坐标、纬度方向分辨率、经度方向分辨率和高程值。

具体的，在步骤S3.1中，利用外部参考DEM(ExtDEM)数据和多视SAR影像的轨道状态矢量、成像参数信息，计算外接矩形数字高程模型(RefDEM)在外部参考DEM数据中的宽nWidth_DEM和高nHeight_DEM，以及外接矩形数字高程模型(RefDEM)的起始点地理坐标(B₁,L₁)。

如图2(a)所示，其中，21表示外部参考DEM数据，22表示多视SAR影像的实际范围，23表示计算得到的外接矩形数字高程模型数据。图2(a)是图2(b)所使用到的外部参考DEM数据，图2(a)的范围要完全包含图2(b)的范围，图2(b)对应的是图2(a)中的“多视SAR影像的实际范围22”所示的框线。

外部参考DEM数据21的范围很大，所以要根据多视SAR影像参数文件(SAR影像的实际范围如图2(a)的22所示)，进行裁剪，只保留与SAR影像范围一致的DEM数据，即计算得到外接矩形数字高程模型(RefDEM)(如图2(a)的23所示)。

如图2(a)所示，要确定外接矩形DEM数据在外部参考DEM数据中的范围，首先是根据图2(a)中多视SAR影像的实际范围22(“斜着的白色”矩形)的4个角点坐标，反算这4个角点坐标在外部参考DEM中的4个新坐标，进而根据这4个新坐标确定外接矩形数字高程模型23(“正着的白色”矩形)的四个新角点坐标，最后根据四个新角点坐标得到宽nWidth_DEM和高nHeight_DEM以及外接矩形数字高程模型23的左上角起始点地理坐标(B₁,L₁)。

步骤S3.1具体包括以下子步骤：

步骤S3.1.1，利用外部参考DEM(ExtDEM)数据和多视SAR影像的轨道状态矢量、成像参数信息，根据SAR严密成像模型(R-D模型)计算得到多视SAR影像的四个角点坐标在外部参考DEM(ExtDEM)数据中的像素坐标；

具体的，利用外部参考DEM(ExtDEM)数据和多视SAR影像的轨道状态矢量、成像参数信息，根据SAR严密成像模型(R-D模型)进行直接地理定位得到多视SAR影像数据的每个角点的地理坐标，所述SAR严密成像模型公式如下式(2)所示，计算多视SAR影像的四个角点坐标在外部参考DEM(ExtDEM)数据中的坐标，得到多视SAR影像数据的四个角点坐标在外部参考DEM(ExtDEM)数据中的四个角点坐标为左上角坐标C_lt'(i,j)、右上角坐标C_rt'(i,j)、左下角坐标C_ld'(i,j)、右下角坐标C_rd'(i,j)；

其中，Vx_S·(X_S-X_P)+Vy_S·(Y_S-Y_P)+Vz_S·(Z_S-Z_P)＝-f_d·λ·(R₀+ΔR·x)/2为多普勒方程，

为距离方程，(X_S,Y_S,Z_S)为卫星位置矢量，(X_P,Y_P,Z_P)为角点的地心空间直角坐标，(V_XS,V_YS,V_ZS)为卫星速度矢量，x为角点的距离向坐标，R₀为近地点斜距，ΔR为距离向分辨率，f_d为多普勒中心频率，λ为雷达波长。在求解的过程中，均是在WGS84椭球系下进行。

步骤S3.1.1具体包括以下子步骤：

步骤S3.1.1.1，根据多视SAR影像的角点(所述角点是指多视SAR影像的左上角点、右上角点、左下角点、右下角点这四个角点中的一个，需要说明的是，这四个角点坐标的计算方法相同)的方位向坐标i_AZ,计算卫星方位向成像时间t_T，并根据卫星方位向成像时间t_T计算卫星位置矢量(X_s,Y_s,Z_s)和速度矢量

具体的，根据下式计算卫星方位向成像时间t_T：t_T＝t₀+t_Azi·i_Az，其中，i_Az表示方位向坐标，t₀表示方位向起始时间，t_Azi表示方位向时间间隔。

其中，方位向起始时间为t₀，方位向时间间隔t_Azi(即，多视SAR影像的成像参数信息)，均由表2成像参数信息提供。优选的，令多视SAR影像的左上角的方位向和距离向坐标为(0,0)，右上角为(0，nMLWidth)，左下角为(nMLHeight，0)，右下角(nMLHeight，nMLWidth)。以多视SAR影像的左上角点(0,0)为例进行说明：多视SAR影像的四个角点方位向和距离向坐标(i_Az,j_Rg)分别为左上角坐标(0,0)、右上角坐标(0,3461)、左下角坐标(7176,0)、右下角坐标(7176,3461)。

根据卫星方位向成像时间t_T计算卫星位置矢量(X_s,Y_s,Z_s)具体包括：

拟合卫星位置矢量与所述卫星方位向成像时间t_T之间的二阶多项式，所述二阶多项式如下公式所示：

其中，(X_S,Y_S,Z_S)为待拟合的卫星位置矢量，

为位置拟合二阶多项式系数矩阵，t_T为卫星方位向成像时间；

以X方向为例，使用表2提供的19个状态位置矢量(即，上文提到的多视SAR影像的轨道状态信息)可构建19个方程，即

X＝(a₀ a₁ a₂)T

其中，X＝(-2628711.7465 -2629649.8296 -2630583.9876 -2631514.2213 -2632440.5295 -2633362.9117 -2634281.3685 -2635195.8984 -2636106.5024 -2637013.1796 -2637915.9289 -2638814.7513 -2639709.6459 -2640600.6118 -2641487.6501 -2642370.7589 -2643249.9397 -2644125.1910 -2644996.5124)'；T＝(t₀t₁ … t₁₂)′＝(80787.910594 80788.910594 80789.910594 80790.910594 80791.91059480792.910594 80793.910594 80794.910594 80795.910594 80796.910594 80797.91059480798.910594 80799.910594 80800.910594 80801.910594 80802.910594 80803.91059480804.910594 80805.910594)，状态矢量的时间间隔为1s。解算可获得多项式系数，

(a_x0 a_x1 a_x2)＝(-2627769.7393 -943.9692 1.9619)

以上计算得到了多项式系数，进而根据多项式系数和卫星方位向成像时间计算任意位置的卫星位置矢量。

根据卫星方位向成像时间t_T计算卫星速度矢量

具体包括：

拟合卫星速度矢量与所述卫星方位向成像时间t_T之间的二阶多项式，所述二阶多项式如公式(4)所示：

其中，

为待拟合的卫星速度矢量，t_T为卫星方位向成像时间，

为速度拟合二阶多项式系数矩阵。

步骤S3.1.1.2，对多视SAR影像的中心点经纬度

进行坐标转换，将其从大地坐标系转换到地心空间直角坐标系，作为待求点的地心坐标初值(X_T ⁰,Y_T ⁰,Z_T ⁰)，所述待求点是指多视SAR影像的左上角点、右上角点、左下角点、右下角点四个角点中的一个(需要说明的是，这四个角点坐标的计算方法相同)；

步骤S3.1.1.3，基于步骤S3.1.1.1求得的卫星位置矢量(X_s,Y_s,Z_s)和速度矢量

和步骤S3.1.1.2求得的待求点的地心坐标初值(X_T ⁰,Y_T ⁰,Z_T ⁰)，根据SAR严密成像模型，解算所述待求点的改正数；

具体的，如上所述，所述SAR严密成像模型为：

其中，Vx_S·(X_S-X_P)+Vy_S·(Y_S-Y_P)+Vz_S·(Z_S-Z_P)＝-f_d·λ·(R₀+ΔR·x)/2为多普勒方程，

为距离方程，(X_S,Y_S,Z_S)为卫星位置矢量，(X_P,Y_P,Z_P)为角点的地心空间直角坐标，

为卫星速度矢量，x为角点的距离向坐标，R₀为近地点斜距，ΔR为距离向分辨率，f_d为多普勒中心频率，λ为雷达波长。

根据上述公式(2)，构建观测方程，计算多视SAR影像的角点的误差方程系数，形成误差方程系数矩阵，再根据误差方程系数矩阵建立法方程(该过程为本领域公知的技术，在此不再赘述)：B^TBX-B^TL＝0，根据下式解算待求点的改正数：

X＝(B^TB)^-1B^TL

其中，X＝(ΔX_T,ΔY_T,ΔZ_T)，ΔX_T，ΔY_T，ΔZ_T是待求点的改正数坐标；B为误差方程系数矩阵(在该步骤的计算中，上述公式(2)中所涉及的卫星位置矢量和速度矢量由步骤S3.1.1.1计算得到)；L为常数项矩阵，由步骤S3.1.1.2计算得到。

步骤S3.1.1.4，根据下式，利用待求点的改正数更新待求点的地心坐标：

X_T ¹＝X_T ⁰+ΔX_T，Y_T ¹＝Y_T ⁰+ΔY_T，Z_T ¹＝Z_T ⁰+ΔZ_T；

其中，

是更新后的待求点的地心坐标，

是待求点的地心坐标初值，ΔX_T，ΔY_T，ΔZ_T是待求点的改正坐标。

步骤S3.1.1.5，利用更新后的待求点的地心坐标(X_T ¹,Y_T ¹,Z_T ¹)进行迭代解算，重复S3.1.1.3-S3.1.1.4，直至改正数(ΔX_T,ΔY_T,ΔZ_T)的绝对值都小于设定阈值ε迭代结束，确定待求点的地心坐标(X_T,Y_T,Z_T)；根据本发明的优选实施例，阈值ε为10^-8。

步骤S3.1.1.6，对通过步骤S3.1.1.5获取的待求点的地心坐标(X_T,Y_T,Z_T)进行坐标转换，将其从地心空间直角坐标系转换到大地坐标系，得到并输出待求点的大地坐标(B_T,L_T,H_T)；

该大地坐标(B_T,L_T,H_T)是多视SAR影像的直接定位结果，作为多视SAR影像待求点的大地坐标(即经纬度坐标)。

步骤S3.1.1.7，根据下述公式(5)，计算得到多视SAR影像的角点坐标在外部参考DEM(ExtDEM)数据中的像素坐标(i_DEM,j_DEM)；

具体的，根据外部参考DEM(ExtDEM)数据提供的起点地理坐标(B_ext,0,L_ext,0)、纬度方向的分辨率B_ext,res、经度方向的分辨率L_ext,res，再结合步骤S3.1.1.6输出的多视SAR影像的角点的大地坐标，利用下述公式(5)进行计算，得到多视SAR影像角点坐标在外部参考DEM(ExtDEM)数据中的像素坐标(i_DEM,j_DEM)；

其中，(i_DEM,j_DEM)为多视SAR影像的角点坐标在外部参考DEM(ExtDEM)数据中的像素坐标；(B_ext,0,L_ext,0)为外部参考DEM(ExtDEM)数据的起点地理坐标；B_ext,res为外部参考DEM(ExtDEM)数据的纬度方向分辨率，L_ext,res为外部参考DEM(ExtDEM)数据的经度方向分辨率；B_i为多视SAR影像的角点坐标对应的大地坐标纬度值，L_i为多视SAR影像角点坐标对应的大地坐标经度值，由步骤S3.1.1.6获得。

需要说明的是，计算多视SAR影像的左上角点、右上角点、左下角点、右下角点这四个角点的大地坐标的过程相同，所以通过步骤S3.1.1.1-S3.1.1.7可以计算得到多视SAR影像的左上角点、右上角点、左下角点、右下角点这四个角点的大地坐标。

根据本发明的具体实施例，外部参考DEM(ExtDEM)数据的起点地理坐标(B_ext,0,L_ext,0)＝(41.00042,114.999583)、纬度方向分辨率为B_ext,res＝-8.3333333e-5、经度方向分辨率为L_ext,res＝8.3333333e-5。由此计算得到的多视SAR影像的四个角点坐标在外部参考DEM(ExtDEM)数据中的坐标分别为左上角C_lt'(i,j)＝(9315,13865)、右上角C_rt'(i,j)＝(12374,19341)、左下角C_ld'(i,j)＝(17310,9605)、右下角C_rd'(i,j)＝(19620,16874)。

步骤S3.1.2，根据步骤S3.1.1获取的多视SAR影像的四个角点在外部参考DEM(ExtDEM)数据中的像素坐标(i_DEM,j_DEM)，计算外接矩形数字高程模型

(RefDEM)的四个角点(即，外接矩形数字高程模型的左上角点、右上角点、左下角点、右下角点)坐标在外部参考DEM(ExtDEM)数据中的坐标，以及外接矩形数字高程模型(RefDEM)在外部参考DEM数据中的宽nWidth_DEM和高nHeight_DEM；

在步骤S3.1.1中得到多视SAR影像四个角点在外部参考DEM(ExtDEM)数据中的像素坐标(i_DEM,j_DEM)，即图2(a)中的附图标记22所示的框图的四个角点在外部参考DEM(ExtDEM)数据中的(i_DEM,j_DEM)。可以看到附图标记22所示的框不是一个规则的矩形，为了计算方便，需要寻找(i_DEM,j_DEM)中的行和列的最大最小值，最终确定图2中的附图标记23所示的框图的像素坐标和该框图的宽、高。

要确定外接矩形数字高程模型(RefDEM)在外部参考DEM(ExtDEM)数据中的边界范围，必须对多视SAR影像的四个角点在外部参考DEM(ExtDEM)数据中的像素坐标(i_DEM,j_DEM)进行判断，寻找它们中的行方向的最大、最小值，列方向的最大、最小值，而行方向的最大值与最小值之差即为宽nWidth_DEM，列方向的最大值与最小值之差即为高nHeight_DEM，进而得到外接矩形数字高程模型(RefDEM)在外部参考DEM数据中的宽nWidth_DEM和高nHeight_DEM；另外，在对多视SAR影像的四个角点在外部参考DEM(ExtDEM)数据中的像素坐标(i_DEM,j_DEM)进行判断，寻找它们中的行方向的最大、最小值，列方向的最大、最小值后，可获得外接矩形数字高程模型(RefDEM)的四个角点坐标在外部参考DEM(ExtDEM)数据中分别为左上角坐标C_lt(i,j)(即行最小，列最小)、右上角坐标C_rt(i,j)(即行最小，列最大)、左下角坐标C_ld(i,j)(即行最大，列最小)、右下角坐标C_rd(i,j)(即行最大，列最大)；

举例来说，根据步骤S3.1.1.7提供的多视SAR影像的四个角点坐标在外部参考DEM(ExtDEM)数据中的坐标分别为左上角C_lt'(i,j)＝(9315,13865)、右上角C_rt'(i,j)＝(12374,19341)、左下角C_ld'(i,j)＝(17310,9605)、右下角C_rd'(i,j)＝(19620,16874)，它们中的行方向最大、最小值分别为19341和9605，它们中的列方向最大、最小值分别为19620和9315，因此外接矩形数字高程模型(RefDEM)的四个角点坐标在外部参考DEM(ExtDEM)数据中分别为左上角C_lt(i,j)＝(9315,9605)(即行最小，列最小)、右上角C_rt(i,j)＝(9315,19341)(即行最小，列最大)、左下角C_ld(i,j)＝(19620,9605)(即行最大，列最小)、右下角C_rd(i,j)＝(19620,19341)(即行最大，列最大)。

步骤S3.1.3，根据步骤S3.1.2得到的外接矩形数字高程模型(RefDEM)的左上角C_lt(i,j)角点坐标，结合外部参考DEM(ExtDEM)数据提供的起点地理坐标(B_ext,0,L_ext,0)(SAR处理领域，一般以左上角为图像的起始点)、纬度方向的分辨率B_ext,res、经度方向的分辨率L_ext,res，由下式(6)计算外接矩形数字高程模型(RefDEM)的起始点地理坐标(B₁,L₁)；

其中，(B₁,L₁)为外接矩形数字高程模型(RefDEM)的起始点地理坐标；

(B_ext,0,L_ext,0)为外部参考DEM(ExtDEM)数据提供的起点地理坐标，B_ext,res为外部参考DEM(ExtDEM)数据的纬度方向的分辨率，L_ext,res为外部参考DEM(ExtDEM)数据的经度方向的分辨率；i为外接矩形数字高程模型(RefDEM)在外部参考DEM(ExtDEM)数据中的左上角的行方向坐标，j为外接矩形数字高程模型(RefDEM)在外部参考DEM(ExtDEM)数据中的左上角的列方向坐标(该行、列坐标通过步骤S3.1.2求得)。

此外，外接矩形数字高程模型(RefDEM)的纬度方向的分辨率(B_ref,res)和经度方向的分辨率(L_ref,res)分别与外部参考DEM(ExtDEM)数据的纬度方向分辨率和经度方向的分辨率相同，即：

B_ref,res＝B_ext,res，

L_ref,res＝L_ext,res。

根据本发明的具体实施例，外接矩形数字高程模型(RefDEM)的宽和高分别为nWidth_DEM＝9736和nHeight_DEM＝10305。外接矩形数字高程模型(RefDEM)在外部参考ExtDEM数据中的范围，如图2中的23所示。相应的外接矩形数字高程模型(RefDEM)的起始点地理坐标为(B₁,L₁)＝(40.224,115.799)。

步骤S3.1.4，根据步骤S3.1.3获取的外接矩形数字高程模型(RefDEM)的起始点地理坐标(B₁,L₁)，计算外接矩形数字高程模型(RefDEM)的高程值；

具体的，根据步骤S3.1.3获取的外接矩形数字高程模型(RefDEM)的起始点地理坐标(B₁,L₁)、外接矩形数字高程模型(RefDEM)的纬度方向的分辨率B_ref,res、外接矩形数字高程模型(RefDEM)的经度方向的分辨率L_ref,res，在外部参考DEM(ExtDEM)数据中取高程值，并将取得的高程值存储于外接矩形数字高程模型(RefDEM)中。

具体的，根据下述公式(7)计算外接矩形数字高程模型(RefDEM)的高程值：

其中，(B_RefDEM,L_RefDEM)为外接矩形数字高程模型(RefDEM)的某个点的经纬度坐标，B_RefDEM为外接矩形数字高程模型(RefDEM)的某个点的纬度坐标，L_RefDEM为外接矩形数字高程模型(RefDEM)的某个点的经度坐标；(B₁,L₁)为外接矩形数字高程模型(RefDEM)的起始点地理坐标，B_ref,res为外接矩形数字高程模型(RefDEM)的纬度方向的分辨率，L_ref,res为外接矩形数字高程模型(RefDEM)的经度方向的分辨率，nWidth_DEM和nHeight_DEM分别为外接矩形数字高程模型(RefDEM)的宽和高，m为高度方向的坐标索引，其中m＝0,1,…,nHeight_DEM-1；n为宽度方向的坐标索引，其中n＝0,1,…,nWidth_DEM-1；(B_ext,0,L_ext,0)为外部参考DEM(ExtDEM)数据的起点地理坐标，B_ext,res为外部参考DEM(ExtDEM)在纬度方向的分辨率，L_ext,res为外部参考DEM(ExtDEM)在经度方向的分辨率；(i_DEM,j_DEM)为外接矩形数字高程模型(RefDEM)的某个点的经纬度坐标在外部参考DEM(ExtDEM)数据中的行列坐标，即，i_DEM为计算得到的外接矩形数字高程模型(RefDEM)的像素坐标在外部参考DEM(ExtDEM)数据中的行坐标，j_DEM为计算得到的外接矩形数字高程模型(RefDEM)的像素坐标在外部参考DEM(ExtDEM)数据中的列坐标；Pix_DEM(i_DEM,j_DEM)表示行列坐标为(i_DEM,j_DEM)的点在外部参考DEM(ExtDEM)数据中的高程值，H_DEM表示外接矩形数字高程模型(RefDEM)的高程值。

根据本发明的具体实施例，外部参考DEM(ExtDEM)数据的起点地理坐标(B_ext,0,L_ext,0)＝(41.00042,114.999583)、纬度方向分辨率为B_ext,res＝-8.3333333e-5、经度方向、分辨率为L_ext,res＝8.3333333e-5。

步骤S3.2，根据步骤S3.1获取的外接矩形数字高程模型(RefDEM)的宽nWidth_DEM和高nHeight_DEM，创建两个大小均为(nHeight_DEM,nWidth_DEM)的实数空矩阵，用于分别存储外接矩形数字高程模型(RefDEM)数据和DEM模拟SAR影像数据；创建一个大小为(nHeight_DEM,nWidth_DEM)的复数空矩阵，用于存储DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系数据。

在图3a，3b，3c中的数据对应的宽和高是根据步骤S3.1获取的；为了后续计算的方便，在步骤S3.2中先创建三个空的文件，创建两个大小均为(nHeight_DEM,nWidth_DEM)的实数空矩阵，分别存储外接矩形数字高程模型(RefDEM)数据(最终结果为图3a，对应的步骤为S3.1)和DEM模拟SAR影像数据(最终结果为图3c，对应的步骤为S3.4)；

创建一个大小为(nHeight_DEM,nWidth_DEM)的复数空矩阵，用于存储DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系(即坐标映射数据1)(最终结果为图3b，对应的步骤为S3.3)。

此外还可以看到，图3c的结果为DEM模拟SAR影像，它的方向与图3a中的RefDEM的方向是一致的；而为了和原始SAR影像进行匹配，必须将图3c转换为与SAR影像一致的方向上，这个处理过程在步骤S4进行。

步骤S3.3，基于所述外部参考DEM数据和多视SAR影像，根据SAR严密成像模型(R-D模型)，计算外接矩形数字高程模型(RefDEM)数据与多视SAR影像雷达坐标空间一一映射的坐标关系；

具体的，根据SAR严密成像模型(R-D模型)，结合由步骤S3.1获得的外接矩形数字高程模型(RefDEM)数据的起始点地理坐标和高程，以及多视SAR影像的成像参数中的卫星状态矢量参数，计算外接矩形数字高程模型(RefDEM)数据中的每个像素点在多视SAR影像雷达坐标空间中的像点坐标(x,y)，所述像点坐标(x,y)为外接矩形数字高程模型(RefDEM)数据与多视SAR影像雷达坐标空间一一映射的坐标关系，通过这个映射关系，即可在两个数据(不同坐标系)中相互找到对应关系。所述SAR严密成像模型如公式(2)所示。

DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系数据如图3(b)所示。因为，DEM模拟SAR影像是根据外接矩形数字高程模型(RefDEM)数据和SAR模拟技术得到的。这个坐标映射数据起着桥梁衔接的作用，通过坐标映射数据，可以实现DEM模拟SAR影像与多视SAR影像雷达坐标空间任何坐标的查找。

在得到DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系数据(如图3(b))之后，即可根据坐标映射数据得到SAR模拟影像(如图4(a)所示，对应步骤S4)。

步骤S3.3具体包括以下子步骤：

S3.3.1，基于公式(2)中的多普勒方程迭代计算外接矩形数字高程模型(RefDEM)数据中的每个像素点的初始方位向坐标y；

步骤S3.3.1具体包括以下步骤：

S3.3.1.1，设定外接矩形数字高程模型(RefDEM)数据中的每个像素点的像素坐标初值为多视SAR影像的中心点(i_center,j_center)(其中j_center＝nMLWidth/2，i_center＝nMLHeight/2，nMLWidth和nMLHeight分别为多视SAR影像的宽和高)，根据多视SAR影像中心点的方位向坐标i_center计算卫星方位向成像时间t_y，计算公式如下式所示：

t_y＝t₀+t_Azi·i_center

其中，t_y为卫星方位向成像时间；i_center为视SAR影像中心点的方位向坐标；t₀为卫星方位向成像起始时间，t_Azi为方位向时间间隔，t₀和t_Azi均由表2提供。

S3.3.1.2，根据所述卫星方位向成像时间t_y计算卫星位置矢量(X_S,Y_S,Z_S)(具体计算过程在步骤S3.1.1.1中已描述，不再赘述)；

S3.3.1.3，根据所述卫星方位向成像时间t_y计算速度矢量

(具体计算过程在步骤S3.1.1.1中已描述，不再赘述)；

S3.3.1.4，利用所述卫星位置矢量(X_S,Y_S,Z_S)、卫星速度矢量

和外接矩形数字高程模型(RefDEM)数据中的每个像素点的地心空间直角坐标(X_P,Y_P,Z_P)，根据多普勒方程计算得到多普勒值f_d(多普勒方程中的雷达波长、近斜距和距离向分辨率参数由表2提供)，同时计算得到多普勒变化率f_d'，结合多普勒测量值f_dop，根据下式计算出需要调整的时间变量Δt：

Δt＝-(f_d-f_dop)/f′_d；

其中，多普勒变化率f′_d是通过数值微分的方法近似求出，假设要计算当前方位向坐标y对应的时间t_y的f′_d，设dt＝0.01s，利用多普勒方程分别计算t_y和t_y+dt处的多普勒值f_d(t_y)和f_d(t_y+dt)，则

f′_d＝(f_d(t_y+dt)-f_d(t_y))/dt

此外，外接矩形数字高程模型(RefDEM)数据中的每个像素点的地心空间直角坐标(X_P,Y_P,Z_P)，它是根据公式(7)计算得到的(B_RefDEM,L_RefDEM,H_DEM)，其中B_RefDEM为该点的纬度，L_RefDEM为该点的经度，H_DEM为该点的高程。随后对(B_RefDEM,L_RefDEM,H_DEM)进行坐标转换，转换到地心空间直角坐标系中得到(X_P,Y_P,Z_P)。

S3.3.1.5，根据下式，利用所述时间变量Δt对卫星成像时间t_y进行修改，得到新的成像时刻t_y+1：

t_y+1＝t_y+Δt；

S3.3.1.6，重复步骤S3.3.1.2-S3.3.1.5，进行迭代计算，直到Δt绝对值小于阈值ε，迭代结束，输出卫星方位向成像时间t_y；

S3.3.1.7，由解算的卫星方位向成像时间t_y，根据下式计算每个像素点的初始方位向坐标y：y＝(t_y-t₀)/Δt_y，其中，t₀为卫星方位向成像起始时间，Δt_y为方位向时间间隔，这两个参量均由表2提供。

S3.3.2，根据公式(2)的距离方程计算外接矩形数字高程模型(RefDEM)数据中的每个像素点的初始距离向坐标x。

由距离方程

结合步骤S3.3.1计算得到的每个像素点的位置矢量(X_P,Y_P,Z_P)，以及由表2提供的近地点斜距R₀、距离向分辨率ΔR，计算每个像素点的初始距离向坐标x。

S3.3.3，组合步骤S3.3.1和步骤S3.3.2得到的外接矩形数字高程模型(RefDEM)数据中的每个像素点的初始方位向坐标y和初始距离向坐标x，得到外接矩形数字高程模型(RefDEM)数据中的每个像素点在多视SAR影像雷达坐标空间中的像点坐标(x,y)，所述像点坐标(x,y)为外接矩形数字高程模型(RefDEM)数据与多视SAR影像雷达坐标空间一一映射的坐标关系。通过这个映射关系，即可在两个数据(坐标系)中相互找到对应关系。

步骤S3.4，根据SAR影像模拟公式(即半经验后向散射模型)，根据外接矩形数字高程模型(RefDEM)数据模拟SAR像元(i,j)处的影像值Pix_Sim(i,j)，最终所有模拟SAR像元值Pix_Sim(i,j)的集合构成DEM模拟SAR影像。SAR影像模拟公式为：

其中，Pix_Sim(i,j)为模拟SAR像元(i,j)处的影像值，σ为雷达散射截面，K为常数，η为当地入射角大小，其中，

R_ts为雷达入射方向矢量；

为地面目标点所在表面的法线矢量。

Pix_Sim(i,j)为模拟SAR像元(i,j)处的影像值，只是模拟SAR图像中的一个像素值，所有Pix_Sim(i,j)组成的像素集合为DEM模拟SAR影像，如图3(C)所示。.

优选的，K＝[1,10]区间内的数。

输出三个数据文件，如图3所示。图3-(a)中保存外接矩形数字高程模型(RefDEM)高程数据；因外接矩形DEM数据方向与多视SAR影像(雷达坐标)方向不一致，需要经过坐标映射才能将它们之间的像元一一对应起来，因此3-(b)中保存的是外接矩形DEM数据与多视SAR雷达坐标空间的对应关系，根据坐标映射数据即可实现相互转换；图3-(c)中保存的是基于DEM和公式(8)得到的所有像元值构成DEM模拟SAR影像，同样可以根据坐标映射数据，将DEM模拟SAR影像与多视SAR影像雷达坐标空间一一对应起来。

步骤S4，根据步骤S3.3获取的DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系和步骤S3.4获取的DEM模拟SAR影像，计算得到与多视SAR影像雷达坐标空间一致的SAR模拟影像，结果如图4a所示。

如图2(b)和3(c)所示，DEM模拟SAR影像与多视SAR影像之间的数据方向是不一致的，需要通过DEM模拟SAR影像和多视SAR影像雷达坐标空间一一映射的坐标关系(即坐标映射数据)实现数据转换，得到与多视SAR影像(雷达坐标)方向一致的SAR模拟影像。遍历坐标映射数据，取其复数像元值，若其复数像元值满足多视SAR影像的宽nMLWidth和高nMLHeight，则将对应的DEM模拟SAR影像的像元值赋值给SAR模拟影像。

步骤S4具体包括以下三个步骤：

步骤S4.1，获得DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系(即坐标映射数据1)的每一个像元pix(i,j)的复数像元值(i_Value,j_Value)，如满足如下条件式，则保留该像元pix(i,j)，参与下一步计算。

其中，i_value为某一个像元pix(i,j)对应的复数像元值的高度方向坐标，j_value为某一个像元pix(i,j)对应的复数像元值的宽度方向坐标，nMLHeight为多视SAR影像的高，nMLWidth为多视SAR影像的宽。

步骤S4.2，根据步骤S4.1获取的像元pix(i,j)和对应的复数像元值(i_Value,j_Value)，在DEM模拟SAR影像中取(i,j)像元值填充到复数像元值(i_Value,j_Value)对应的SAR模拟影像坐标中。

步骤S4.3，对DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系(即坐标映射数据1)的每一个像元均采用步骤S4.1和步骤S4.2，得到SAR模拟影像，结果如图4-(a)所示。图4-(b)为多视SAR影像，现在它们的影像方向与大小是一致的，可以参与下一步的影像匹配。

步骤S5，对所述SAR模拟影像和所述多视SAR影像进行灰度匹配，计算得到新的DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系；

具体的，根据步骤S4得到的SAR模拟影像(简称SSAR)和多视SAR影像(简称MSAR)，对SAR模拟影像与多视SAR影像在一定范围的窗口内进行灰度匹配，得到方位向和距离向的偏移量，并基于所述方位向和距离向的偏移量更新DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系(即坐标映射数据1)得到新的DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系(即坐标映射数据2)。

步骤S5，具体包括以下三个子步骤：

步骤S5.1，采用匹配方法，在SAR模拟影像上规则选取一定数量的点，对搜索窗口数为nWin范围内的像素进行匹配，得到窗口内每个像元的匹配测度ρ，取最大的匹配测度ρ，以及对应的方位向偏移量记为off_Az，距离向偏移量记为off_Rg。匹配方法采用常用的NCC匹配方法(归一化互相关系数)，其公式如下：

其中，r_i表示多视SAR影像(参考影像)窗口内的像素值；m_i表示SAR模拟影像(待匹配影像)窗口内的像素值；

表示多视SAR影像窗口内的像素均值；

表示SAR模拟影像窗口内的像素均值；σ_r表示多视SAR影像窗口内的像素标准差；σ_m表示SAR模拟影像窗口内的像素标准差；n表示匹配窗口内的像素个数，n＝nWin²，nWin为搜索窗口大小。

优选的，nWin＝512。

步骤S5.2，对步骤S5.1获取的方位向的偏移量off_Az、距离向的偏移量off_Rg，进行多项式拟合，得到方位向和距离向的偏移量拟合系数。多项式拟合公式为：

其中，y(x_Az)、y(x_Rg)为方位向和距离向的偏移量矩阵，矩阵大小为k个元素；x_Az、x_Rg为对应的元素索引，取值为x_Az＝1,2,...,k，x_Rg＝1,2,...,k；a₀、a₁为方位向的多项式常数项和一次项系数；b₀、b₁为距离向的多项式常数项和一次项系数。

优选的，方位向偏移量系数为常数项a₀＝0.77653、一次项a₁＝1.19188e-005；距离向偏移量系数为常数项b₀＝0.09955、一次项b₁＝2.58910e-004。

步骤S5.3，根据步骤S5.2获取的方位向和距离向的偏移量拟合系数，更新坐标映射数据1得到新的坐标映射数据2。更新的公式如下：

对坐标映射数据1的每一个像元(i,j)，得到其复数像元值(i_Value,1,j_Value,1)，根据方位向和距离向的偏移量，有

其中，i_Value,2、j_Value,2为更新后的坐标映射数据2的方位向和距离向像元值；i_Value,1、j_Value,1为坐标映射数据1的方位向和距离向像元值；a₀、a₁为方位向的多项式常数项和一次项系数；b₀、b₁为距离向的多项式常数项和一次项系数。

步骤S6，基于所述新的DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系，计算资源三号卫星影像控制点数据在多视SAR影像中的像点坐标集；

具体的，根据步骤S3.3获取的外接矩形数字高程模型(RefDEM)数据的起始角点地理坐标(B_ref,1,L_ref,1)，结合资源三号卫星影像控制点数据ZY3，逐点计算ZY3数据在外接矩形数字高程模型(RefDEM)中的坐标(i_DEM,j_DEM)，其次利用步骤S5.3获取的坐标映射数据2，得到ZY3数据在多视SAR影像中的像点坐标(i_MSAR,j_MSAR)，构成ZY3控制点数据在多视SAR影像中的像点坐标集(B_i,L_i,H_i,i_MSAR,j_MSAR)。

步骤S6具体包括以下三个子步骤：

步骤S6.1，根据步骤S3.3获取的外接矩形数字高程模型(RefDEM)数据的起始角点地理坐标(B₁,L₁)，结合资源三号卫星影像控制点数据ZY3，见表1所示，逐点计算ZY3数据在外接矩形数字高程模型(RefDEM)中的坐标(i_DEM,j_DEM)，计算公式如下所示：

其中，(B_ref,1,L_ref,1)为外接矩形数字高程模型(RefDEM)的起始点地理坐标；(B_i,L_i)为资源三号卫星影像控制点数据ZY3的经纬度地理坐标，如表1所示；B_ref,res为纬度方向外接矩形数字高程模型(RefDEM)分辨率，L_ref,res为经度方向外接矩形数字高程模型(RefDEM)分辨率；(i_DEM,j_DEM)为ZY3数据在外接矩形RefDEM中的坐标。

优选的，外接矩形数字高程模型(RefDEM)的起始点地理坐标为(B_ref,1,L_ref,1)＝(40.224,115.799)度；(B_i,L_i)为表1中纬度和经度信息；外接矩形数字高程模型(RefDEM)的纬度方向分辨率为B_ref,res＝-8.3333e-005度；外接矩形数字高程模型(RefDEM)的经度方向分辨率为L_ref,res＝8.3333e-005度。

控制点名称	纬度/度	经度/度	i<sub>DEM</sub>	j<sub>DEM</sub>
					F03	39.80890	116.09411	4981.22	3541.334
F18	39.74416	116.01044	5758.10	2537.29
					F05	39.79243	116.16860	5178.86	4435.218
C24	39.88796	116.01093	4032.50	2543.17
					C22	39.89176	116.02806	3986.90	2748.731
W15	40.12384	116.22403	1201.92	5100.38
					C02	40.12662	116.22181	1168.56	5073.74
W17	40.06218	116.28779	1941.85	5865.503

表4 ZY3数据在外接矩形数字高程模型(RefDEM)中的坐标

步骤S6.2，利用步骤S6.1获取的ZY3数据在外接矩形数字高程模型(RefDEM)中的坐标(i_DEM,j_DEM)，结合步骤S5.3获取的坐标映射数据2，计算(i_DEM,j_DEM)在SAR模拟影像中的像素坐标(i_MSAR,j_MSAR)。因为根据步骤S3获取的外接矩形数字高程模型(RefDEM)和坐标映射数据2，它们之间的像素关系是一一对应的，如图3所示，即仅仅需要把外接矩形数字高程模型(RefDEM)的像素坐标在坐标映射数据2中将Pix(i_DEM,j_DEM)的像素值输出即可，得到(i_MSAR,j_MSAR)。

优选的，以表1中的F03点为例，计算得到F03点在多视SAR影像中的像点坐标为(i_MSAR,j_MSAR)＝(1285.46,5329.784)。

步骤S6.3，将所有ZY3数据重复步骤S6.1和S6.2，获得ZY3在SAR模拟影像中的像点坐标(i_MSAR,j_MSAR)，同时由于SAR模拟影像已与多视SAR影像进行匹配，它们之间的像素关系是一一对应的，即ZY3在多视SAR影像中的像点坐标(i_MSAR,j_MSAR)是可以直接使用的，不需要进行任何转换，并基于所述像点坐标(i_MSAR,j_MSAR)和ZY3的地理坐标和高程构成ZY3控制点数据在多视SAR影像中的坐标集(B_i,L_i,H_i,i_MSAR,j_MSAR)。其中，(B_i,L_i,H_i)为ZY3的地理坐标和高程；(i_MSAR,j_MSAR)为ZY3在多视SAR影像中的像点坐标。

优选的，ZY3控制点数据在多视SAR影像中的坐标集(B_i,L_i,H_i,i_MSAR,j_MSAR)，如表5所示。

控制点名称	纬度/度	经度/度	i<sub>MSAR</sub>	j<sub>MSAR</sub>
					F03	39.80890	116.09411	1285.46	5329.784
F18	39.74416	116.01044	1872.99	6389.175
					F05	39.79243	116.16860	632.348	5414.823
C24	39.88796	116.01093	2134.48	4395.387
					C22	39.89176	116.02806	1993.94	4309.637
W15	40.12384	116.22403	831.600	718.400
					C02	40.12662	116.22181	856.207	684.1301
W17	40.06218	116.28779	180.838	1448.448

表5 ZY3数据在多视SAR影像中的坐标

步骤S7，将所述资源三号卫星影像控制点数据在多视SAR影像中的像点坐标集转换为所述资源三号卫星影像控制点数据在原始SAR影像中的坐标集；

根据步骤S2设置的SAR影像的方位向(行方向)多视数为nAzML，距离向(列方向)多视数为nRgML，需要将资源三号卫星影像控制点数据在多视SAR影像中的坐标集(B_i,L_i,H_i,i_MSAR,j_MSAR)转换为资源三号卫星影像控制点数据在原始SAR影像(即，真实SAR影像)坐标集，即多视SAR影像的像素坐标分别乘以方位向多视数为nAzML，距离向多视数为nRgML，最后构成ZY3控制点数据在原始SAR影像中的坐标集(B_i,L_i,H_i,i_SAR,j_SAR)。计算公式为：

其中，(i_SAR,j_SAR)为真实SAR像素坐标，(i_MSAR,j_MSAR)为多视SAR像素坐标，nAzML为方位向多视数，nRgML为距离向多视数。

优选的，多视数nAzML＝4，nRgML＝4；计算得到ZY3控制点数据在真实SAR影像中的坐标集(B_i,L_i,H_i,i_SAR,j_SAR)，如表6所示。

控制点名称	纬度/度	经度/度	i<sub>SAR</sub>	j<sub>SAR</sub>
					F03	39.80890	116.09411	5141.87	21319.14
F18	39.74416	116.01044	7491.97	25556.70
					F05	39.79243	116.16860	2529.39	21659.29
C24	39.88796	116.01093	8537.95	17581.55
					C22	39.89176	116.02806	7975.77	17238.55
W15	40.12384	116.22403	3326.40	2873.600
					C02	40.12662	116.22181	3424.83	2736.521
W17	40.06218	116.28779	723.354	5793.793

表6 ZY3数据在真实SAR影像中的坐标

步骤S8，利用资源三号卫星影像控制点数据在原始SAR影像中的坐标集，进行SAR影像正射纠正处理，计算得到距离向和方位向多项式的改正系数；

具体的，利用步骤S7获取的ZY3控制点数据集(B_i,L_i,H_i,i_SAR,j_SAR)和步骤S1提供轨道状态信息，进行SAR影像正射纠正处理，计算得到距离向多项式的改正系数e₀和e₁，方位向多项式的改正系数f₀和f₁，参与后续计算。

SAR影像正射纠正的目的是建立目标点大地坐标与SAR影像坐标之间精确的映射关系。控制点的大地坐标(B,L,H)经过步骤S3.4构建的RD模型间接定位后得到控制点影像坐标(x',y')，这是利用定位模型得到的初步定位结果。显然对于一个精确的定位模型来说，(x',y')应该和实际测量的控制点影像坐标(x,y)是相同的。但由于直接利用SAR轨道参数信息计算得到的定位模型参数通常会存在一定的误差，实际上(x',y')与(x,y)并不相等。它们之间就建立起了目标点的初步定位结果(x',y')与其实际影像坐标(x,y)之间的关系。它们之间的几何关系表示为：

其中，(x,y)是控制点在影像上的实测坐标，(x',y')是利用RD模型由地理坐标计算得到的影像坐标，e₀,e₁分别为距离向多项式的常数项和一次项，f₀,f₁分别为方位向多项式的常数项和一次项。

步骤S8具体包括以下四个子步骤：

步骤S8.1，构建多项式误差模型，分别对SAR影像的方位向和距离向进行多项式系数的计算。考虑到SAR成像时，距离向和方位向近似垂直，两个方向上的误差相对独立，因此在影像的方位向和距离向上分别对影像坐标进行一次线性改正，表达式为：

其中，F(x)为距离向多项式函数，F(y)为方位向多项式函数，(x,y)是控制点在影像上的实测坐标，(x',y')是利用RD模型由地理坐标计算得到的影像坐标，e₀,e₁分别为距离向多项式的常数项和一次项，f₀,f₁分别为方位向多项式的常数项和一次项。

步骤S8.2，根据步骤S8.1的式(10)可以对ZY3控制点列如下线性方程，根据最小二乘平差求解影像的一次线性改正参数。解算误差方程如下：

其中，v_x、v_y分别为函数F(x)、F(y)的残差；

为函数F(x)对e₀的偏导数，

为函数F(x)对e₁的偏导数；

为函数F(y)对f₀的偏导数，

为函数F(y)对f₁的偏导数；Δe₀为距离向多项式的常数项改正量，Δe₁为距离向多项式的一次项改正量；Δf₀为方位向多项式的常数项改正量，Δf₁为方位向多项式的一次项改正量；F_x为已知的距离向坐标，F_y为已知的方位向坐标。

步骤S8.3，根据式(11)所列出的误差模型，列出距离向x的误差方程式

V₁＝A₁Δx₁-l₁ (18)

其中，

Δx₁＝[Δe₀ Δe₁]^T

其中，V₁为残差项矩阵，其中v_xG(i)表示第G(i)个控制点在F(x)函数上的残差，m_G为ZY3控制点的个数；A₁为距离向的偏导数系数矩阵，其中

为函数F(x)对e₀的偏导数，

为函数F(x)对e₁的偏导数；l₁为常数项矩阵，其中F_x0G(i)表示第G(i)个控制点在F(x)函数上的模型计算值；Δx₁为待求解的参数改正量矩阵，Δe₀,Δe₁分别为距离向多项式的常数项改正量和一次项改正量。

优选的，距离向多项式的常数项为e₀＝0.0741和一次项为e₁＝0.999985。

步骤S8.4，根据式(11)所列出的误差模型，列出方位向y的误差方程式

V₂＝A₂Δx₂-l₂ (19)

其中，

Δx₂＝[Δf₀ Δf₁]^T

其中，V₂为残差项矩阵，其中v_yG(i)表示第G(i)个控制点在F(y)函数上的残差，m_G为ZY3控制点的个数；A₂为距离向的偏导数系数矩阵，其中

为函数F(y)对f₀的偏导数，

为函数F(y)对f₁的偏导数；l₂为常数项矩阵，其中F_y0G(i)表示第G(i)个控制点在F(y)函数上的模型计算值；Δx₂为待求解的参数改正量矩阵，Δf₀,Δf₁分别为方位向多项式的常数项改正量和一次项改正量。

优选的，方位向多项式的常数项为f₀＝-0.2013和一次项为f₁＝1.000015。

步骤S9，根据所述距离向和方位向多项式的改正系数，对原始SAR影像进行SAR正射纠正，得到SAR正射纠正后的DOM影像。

具体的，根据步骤S8的公式(15)，结合步骤S3.4构建的RD模型，对原始SAR影像的每个像素进行校正，输出数字正射影像图(Digital Orthophoto Map,DOM)影像和定位精度报告。

步骤S9具体包括以下四个子步骤：

步骤S9.1，根据原始SAR影像的四个角点像素坐标、原始SAR成像参数和输入的DOM数据的纬度方向分辨率、经度方向分辨率，确定待校正的DOM数据的宽和高，以及DOM起始点地理坐标。步骤S9.1包括以下子步骤：

步骤S9.1.1，根据原始SAR影像的四个角点像素坐标和原始SAR成像参数，利用公式(2)和步骤S3.1.1，得到原始SAR影像的四个角点像素坐标对应的四个地理坐标(计算过程和步骤S3.1.1类似，这里不再赘述)。

优选的，原始SAR影像的四个角点像素坐标对应的地理坐标分别为左上角(40.160378,116.333731)、右上角(40.216899,115.925287)、左下角(39.658992,116.212917)、右下角(39.720849,115.815612)。

步骤S9.1.2，根据四个地理坐标计算得到纬度方向的最大最小值和经度方向的最大最小值，进而得到纬度方向的纬度差值dB_DOM和经度方向的经度差值dL_DOM和DOM起始点地理坐标(B_DOM,L_DOM)。

优选的，纬度方向的纬度差值dB_DOM＝0.557907，经度方向的经度差值dL_DOM＝0.518119，DOM起始点地理坐标为纬度方向最大值和经度方向最小值得组合，即(B_DOM,L_DOM)＝(40.216899,115.815612)。

步骤S9.1.3，根据步骤步骤S9.1.2获取的纬度方向的纬度差值dB_DOM和经度方向的经度差值dL_DOM，计算得到待校正的DOM数据的宽nWidth_DOM和高nHeight_DOM，公式如下：

nWidth_DOM＝dB_DOM/B_DOM,res

nHeight_DOM＝dL_DOM/L_DOM.res

其中，nWidth_DOM为待校正的DOM数据的宽，nHeight_DOM为待校正的DOM数据的高；dB_DOM为纬度方向的纬度差值，dL_DOM为经度方向的经度差值；B_DOM,res为DOM数据的纬度方向分辨率，L_DOM.res为DOM数据的经度方向分辨率。

优选的，输入的DOM数据的纬度方向分辨率和经度方向分辨率均为0.000052度；待校正的DOM数据的宽nWidth_DOM＝10002和高nHeight_DOM＝10770。

步骤S9.2，根据步骤S9.1的地理参数(包括DOM起始点地理坐标、输入的DOM数据的纬度方向分辨率、经度方向分辨率)，利用步骤S3.4构建的RD模型，逐像素对待校正的DOM影像进行计算，得到待校正的DOM像素在原始SAR影像中的像素坐标(x',y')。

步骤S9.3，对步骤S9.2得到的像素坐标(x',y')进行改正，得到改正后的原始SAR像素坐标(x,y)，获取原始SAR像素坐标(x,y)对应的像元值Pix_DOM(x,y)，将所述像元值Pix_DOM(x,y)保存到步骤S9.2所述的待校正的DOM像素坐标pix(x_DOM,y_DOM)中，遍历整个待校正的DOM像素，获得DOM正射纠正影像。

具体的，根据步骤S8得到的方位向和距离向多项式系数和步骤S9.2得到的像素坐标(x',y')，利用步骤S8的公式(15)进行改正，得到改正后的原始SAR像素坐标(x,y)，获取原始SAR像素坐标(x,y)对应的像元值Pix_DOM(x,y)，将所述像元值Pix_DOM(x,y)保存到步骤S9.2所述的待校正的DOM像素坐标pix(x_DOM,y_DOM)中，遍历整个待校正的DOM像素，获得DOM正射纠正影像。

根据步骤S8得到的方位向和距离向多项式系数和步骤S9.2得到的像素坐标(x',y')，利用步骤S8的公式(15)进行改正，得到改正后的SAR像素坐标(x,y)，如公式下式所示(公式15)

x＝e₀+e₁·x'

y＝f₀+f₁·y'

就是说，已知e₀,e₁分别为距离向多项式的常数项和一次项，f₀,f₁分别为方位向多项式的常数项和一次项，将通过步骤S9.2得到的像素坐标(x',y')，带入上述公式计算即可得到改正后的SAR像素坐标(x,y)。

优选的，步骤S9还可以包括：

步骤S9.4，对步骤S9.3得到的DOM正射纠正影像进行定位精度评价，输出校正后的精度报告。

优选的，根据表1提供的数据，对步骤S9.3输出的DOM正射纠正数据进行精度评价，统计得到纬度方向中误差为3.95m，经度方向中误差为3.01m，总的平面定位中误差为4.97m。对正射纠正影像DOM的定位精度进行统计，如表7所示。

表7 正射纠正影像DOM的定位精度统计

相对于现有技术，本发明通过外部DEM模拟SAR影像得到SAR模拟影像，并与多视SAR影像进行匹配，确定了高精度资源三号控制点库在SAR影像上的坐标集，比人工刺点精度要高；常规方法使用DEM参与SAR正射纠正精度较低，而现在通过模拟匹配的方法，使用高精度资源三号控制点库参与正射纠正，既解决了选点问题又达到了较好的纠正精度。通过匹配的方法，确定了高精度资源三号控制点库在SAR影像上的坐标集。

以上内容仅为本发明的较佳实施例，对于本领域的普通技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，本说明书内容不应理解为对本发明的限制。

Claims (10)

1.一种基于卫星控制点库和DEM的SAR影像正射纠正方法，其特征在于，所述方法具体包括以下步骤：

步骤S1，采集原始数据，所述原始数据包括资源三号卫星影像控制点数据、外部参考DEM数据、原始SAR影像以及原始SAR影像成像参数；

步骤S2，对所述原始SAR影像进行多视处理，得到多视SAR影像；

步骤S3，基于所述外部参考DEM数据和多视SAR影像，计算外接矩形数字高程模型RefDEM数据、DEM模拟SAR影像、DEM模拟SAR影像和多视SAR影像雷达坐标空间一一映射的坐标关系；

步骤S4，根据步骤S3获取的DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系和步骤S3获取的DEM模拟SAR影像，计算得到与多视SAR影像雷达坐标空间一致的SAR模拟影像；

步骤S5，对所述SAR模拟影像和所述多视SAR影像进行灰度匹配，计算得到新的DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系；

步骤S6，基于所述新的DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系，计算资源三号卫星影像控制点数据在多视SAR影像中的像点坐标集；

步骤S7，将所述资源三号卫星影像控制点数据在多视SAR影像中的像点坐标集转换为所述资源三号卫星影像控制点数据在原始SAR影像中的坐标集；

步骤S8，利用资源三号卫星影像控制点数据在原始SAR影像中的坐标集，进行SAR影像正射纠正处理，计算得到距离向和方位向多项式的改正系数；

步骤S9，根据所述距离向和方位向多项式的改正系数，对所述原始SAR影像进行SAR正射纠正，得到SAR正射纠正后的DOM影像。

2.根据权利要求1所述的基于卫星控制点库和DEM的SAR影像正射纠正方法，其特征在于，步骤2具体包括：根据下式对所述原始SAR影像进行多视处理，得到多视SAR影像：

其中，Pix(i,j)表示多视SAR影像中的像元(i,j)对应的像元值；nAzML表示原始SAR影像的方位向多视数，nRgML表示原始SAR影像的距离向多视数；Pix(n,m)表示原始SAR影像中的像元(n,m)对应的像元值，n表示原始SAR影像中的像元(n,m)的方位向坐标，取值范围为n＝1,K,nAzML，m表示原始SAR影像中的像元(n,m)的距离向坐标，取值范围为m＝1,K,nRgML。

3.根据权利要求1所述的基于卫星控制点库和DEM的SAR影像正射纠正方法，其特征在于，所述步骤S3包括以下子步骤：

步骤S3.1，基于所述外部参考DEM数据和多视SAR影像，计算外接矩形数字高程模型RefDEM数据，其中，所述外接矩形数字高程模型RefDEM数据包括所述外接矩形数字高程模型RefDEM在外部参考DEM数据中的宽、高、起始点地理坐标、纬度方向分辨率、经度方向分辨率和高程值；

步骤S3.2，根据步骤S3.1获取的外接矩形数字高程模型RefDEM的宽nWidth_DEM和高nHeight_DEM，创建两个大小均为(nHeight_DEM,nWidth_DEM)的实数空矩阵，用于分别存储外接矩形数字高程模型RefDEM数据和DEM模拟SAR影像数据；创建一个大小为(nHeight_DEM,nWidth_DEM)的复数空矩阵，用于存储DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系数据；

步骤S3.3，基于所述外部参考DEM数据和多视SAR影像，根据SAR严密成像模型，计算外接矩形数字高程模型RefDEM数据与多视SAR影像雷达坐标空间一一映射的坐标关系；

步骤S3.4，根据SAR影像模拟公式，根据外接矩形数字高程模型RefDEM数据模拟SAR像元(i,j)处的影像值Pix_Sim(i,j)，最终所有模拟SAR像元值Pix_Sim(i,j)的集合构成DEM模拟SAR影像。

4.根据权利要求3所述的基于卫星控制点库和DEM的SAR影像正射纠正方法，其特征在于，所述步骤S3.1可以包括：

步骤S3.1.1，利用外部参考DEM数据和多视SAR影像的轨道状态、成像参数信息，根据SAR严密成像模型计算得到多视SAR影像的四个角点坐标在外部参考DEM数据中的像素坐标；

步骤S3.1.2，根据步骤S3.1.1获取的多视SAR影像的四个角点在外部参考DEM数据中的像素坐标(i_DEM,j_DEM)，计算外接矩形数字高程模型RefDEM的四个角点坐标在外部参考DEM数据中的坐标，以及外接矩形数字高程模型RefDEM在外部参考DEM数据中的宽nWidth_DEM和高nHeight_DEM；

步骤S3.1.3，根据步骤S3.1.2得到的外接矩形数字高程模型RefDEM的左上角C_lt(i,j)角点坐标，结合外部参考DEM数据提供的起点地理坐标(B_ext,0,L_ext,0)、纬度方向的分辨率B_ext,res、经度方向的分辨率L_ext,res，根据下式计算外接矩形数字高程模型RefDEM的起始点地理坐标(B₁,L₁)；

B₁＝B_ext,0+B_ext,res*i

L₁＝L_ext,0+L_ext,res*j

其中，(B₁,L₁)为外接矩形数字高程模型RefDEM的起始点地理坐标；(B_ext,0,L_ext,0)为外部参考DEM数据提供的起点地理坐标，B_ext,res为外部参考DEM数据的纬度方向的分辨率，L_ext,res为外部参考DEM数据的经度方向的分辨率；i为外接矩形数字高程模型RefDEM在外部参考DEM数据中的左上角的行方向坐标，j为外接矩形数字高程模型RefDEM在外部参考DEM数据中的左上角的列方向坐标；

步骤S3.1.4，根据步骤S3.1.3获取的外接矩形数字高程模型RefDEM的起始点地理坐标(B₁,L₁)，计算外接矩形数字高程模型RefDEM的高程值。

5.根据权利要求4所述的基于卫星控制点库和DEM的SAR影像正射纠正方法，其特征在于，所述步骤S3.1可以包括以下子步骤：

步骤S3.1.1.1，根据多视SAR影像的角点的方位向坐标i_AZ,计算卫星方位向成像时间t_T，并根据卫星方位向成像时间t_T计算卫星位置矢量(X_s,Y_s,Z_s)和速度矢量

步骤S3.1.1.2，对多视SAR影像的中心点经纬度

进行坐标转换，将其从大地坐标系转换到地心空间直角坐标系，作为待求点的地心坐标初值(X_T ⁰,Y_T ⁰,Z_T ⁰)；

步骤S3.1.1.3，基于步骤S3.1.1.1求得的卫星位置矢量(X_s,Y_s,Z_s)和速度矢量

和步骤S3.1.1.2求得的待求点的地心坐标初值(X_T ⁰,Y_T ⁰,Z_T ⁰)，根据SAR严密成像模型，解算所述待求点的改正数；

步骤S3.1.1.4，根据下式，利用待求点的改正数更新待求点的地心坐标：

X_T ¹＝X_T ⁰+ΔX_T，Y_T ¹＝Y_T ⁰+ΔY_T，Z_T ¹＝Z_T ⁰+ΔZ_T；

其中，

是更新后的待求点的地心坐标，

是待求点的地心坐标初值，ΔX_T，ΔY_T，ΔZ_T是待求点的改正坐标；

步骤S3.1.1.5，利用更新后的待求点的地心坐标(X_T ¹,Y_T ¹,Z_T ¹)进行迭代解算，重复S3.1.1.3-S3.1.1.4，直至改正数(ΔX_T,ΔY_T,ΔZ_T)的绝对值都小于设定阈值ε迭代结束，确定待求点的地心坐标(X_T,Y_T,Z_T)；

步骤S3.1.1.6，对通过步骤S3.1.1.5获取的待求点的地心坐标(X_T,Y_T,Z_T)进行坐标转换，将其从地心空间直角坐标系转换到大地坐标系，得到并输出待求点的大地坐标(B_T,L_T,H_T)；

步骤S3.1.1.7，根据下述公式，计算得到多视SAR影像的角点坐标在外部参考DEM数据中的像素坐标(i_DEM,j_DEM)：

i_DEM＝(B_ext,0-B_i)/B_ext,res

j_DEM＝(L_ext,0-L_i)/L_ext,res

其中，(i_DEM,j_DEM)为多视SAR影像的角点坐标在外部参考DEM数据中的像素坐标；(B_ext,0,L_ext,0)为外部参考DEM数据的起点地理坐标；B_ext,res为纬度方向分辨率，L_ext,res为经度方向分辨率；B_i为多视SAR影像的角点坐标对应的大地坐标纬度值，L_i为多视SAR影像角点坐标对应的大地坐标经度值。

6.根据权利要求5所述的基于卫星控制点库和DEM的SAR影像正射纠正方法，其特征在于，步骤S3.1.1.3包括：

所述SAR严密成像模型为：

其中，Vx_S·(X_S-X_P)+Vy_S·(Y_S-Y_P)+Vz_S·(Z_S-Z_P)＝-f_d·λ·(R₀+ΔR·x)/2为多普勒方程，

为距离方程，(X_S,Y_S,Z_S)为卫星位置矢量，(X_P,Y_P,Z_P)为角点的地心空间直角坐标，

为卫星速度矢量，x为角点的距离向坐标，R₀为近地点斜距，ΔR为距离向分辨率，f_d为多普勒中心频率，λ为雷达波长；

根据上述公式(2)，构建观测方程，计算多视SAR影像的角点的误差方程系数，形成误差方程系数矩阵，再根据误差方程系数矩阵建立法方程：B^TBX-B^TL＝0，根据下式解算待求点的改正数：

X＝(B^TB)^-1B^TL

其中，X＝(ΔX_T,ΔY_T,ΔZ_T)，ΔX_T，ΔY_T，ΔZ_T是待求点的改正数坐标；B为误差方程系数矩阵；L为常数项矩阵；

其中，在步骤S3.1.4中，根据下述公式计算外接矩形数字高程模型RefDEM的高程值：

B_RefDEM＝B₁+B_ref,res*m,m＝0,1,L,nHeight_DEM-1

L_RefDEM＝L₁+L_ref,res*n,n＝0,1,L,nWidth_DEM-1

i_DEM＝(B_ext,0-B_RefDEM)/B_ext,res

j_DEM＝(L_ext,0-L_RefDEM)/L_ext,res

H_DEM＝Pix_DEM(i_DEM,j_DEM)

其中，(B_RefDEM,L_RefDEM)为外接矩形数字高程模型RefDEM的某个点的经纬度坐标，B_RefDEM为外接矩形数字高程模型RefDEM的某个点的纬度坐标，L_RefDEM为外接矩形数字高程模型RefDEM的某个点的经度坐标；(B₁,L₁)为外接矩形数字高程模型RefDEM的起始点地理坐标，B_ref,res为外接矩形数字高程模型RefDEM的纬度方向的分辨率，L_ref,res为外接矩形数字高程模型RefDEM的经度方向的分辨率，nWidth_DEM和nHeight_DEM分别为外接矩形数字高程模型RefDEM的宽和高，m为高度方向的坐标索引，其中m＝0,1,L,nHeight_DEM-1；n为宽度方向的坐标索引，其中n＝0,1,L,nWidth_DEM-1；(B_ext,0,L_ext,0)为外部参考DEM数据的起点地理坐标，B_ext,res为外部参考DEM(ExtDEM)在纬度方向的分辨率，L_ext,res为外部参考DEM在经度方向的分辨率；(i_DEM,j_DEM)为外接矩形数字高程模型RefDEM的某个点的经纬度坐标在外部参考DEM数据中的行列坐标；Pix_DEM(i_DEM,j_DEM)表示行列坐标为(i_DEM,j_DEM)的点在外部参考DEM数据中的高程值，H_DEM表示外接矩形数字高程模型RefDEM的高程值。

7.根据权利要求3所述的影像正射纠正方法，其特征在于，步骤S3.3包括以下子步骤：

步骤S3.3.1，基于公式(2)中的多普勒方程迭代计算外接矩形数字高程模型RefDEM数据中的每个像素点的初始方位向坐标y；

S3.3.2，根据公式(2)的距离方程计算外接矩形数字高程模型RefDEM数据中的每个像素点的初始距离向坐标x；

S3.3.3，组合步骤S3.3.1和步骤S3.3.2得到的外接矩形数字高程模型RefDEM数据中的每个像素点的初始方位向坐标y和初始距离向坐标x，得到外接矩形数字高程模型RefDEM数据中的每个像素点在多视SAR影像雷达坐标空间中的像点坐标(x,y)，所述像点坐标(x,y)为外接矩形数字高程模型RefDEM数据与多视SAR影像雷达坐标空间一一映射的坐标关系。

8.根据权利要求1所述的影像正射纠正方法，其特征在于，

步骤S4具体包括以下三个步骤：

步骤S4.1，获得DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系的每一个像元pix(i,j)的复数像元值(i_Value,j_Value)，如满足如下条件式，则保留该像元pix(i,j)，参与下一步计算：

其中，i_value为某一个像元pix(i,j)对应的复数像元值的高度方向坐标，j_value为某一个像元pix(i,j)对应的复数像元值的宽度方向坐标，nMLHeight为多视SAR影像的高，nMLWidth为多视SAR影像的宽；

步骤S4.2，根据步骤S4.1获取的像元pix(i,j)和对应的复数像元值(i_Value,j_Value)，在DEM模拟SAR影像中取(i,j)像元值填充到复数像元值(i_Value,j_Value)对应的SAR模拟影像坐标中；

步骤S4.3，对DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系的每一个像元均采用步骤S4.1和步骤S4.2，得到SAR模拟影像；

步骤S5具体包括以下三个子步骤：

步骤S5.1，采用匹配方法，在SAR模拟影像上规则选取一定数量的点，对搜索窗口数为nWin范围内的像素进行匹配，得到窗口内每个像元的匹配测度ρ，取最大的匹配测度ρ，以及对应的方位向偏移量记为off_Az，距离向偏移量记为off_Rg；匹配方法采用常用的NCC匹配方法，其公式如下：

其中，r_i表示多视SAR影像窗口内的像素值；m_i表示SAR模拟影像窗口内的像素值；

表示多视SAR影像窗口内的像素均值；

表示模拟SAR影像窗口内的像素均值；σ_r表示多视SAR影像窗口内的像素标准差；σ_m表示模拟SAR影像窗口内的像素标准差；n表示匹配窗口内的像素个数，n＝nWin²，nWin为搜索窗口大小；

步骤S5.2，对步骤S5.1获取的方位向的偏移量off_Az、距离向的偏移量off_Rg，进行多项式拟合，得到方位向和距离向的偏移量拟合系数；多项式拟合公式为：

y(x_Az)＝a₀+a₁x_Az

y(x_Rg)＝b₀+b₁x_Rg

其中，y(x_Az)、y(x_Rg)为方位向和距离向的偏移量矩阵，矩阵大小为k个元素；x_Az、x_Rg为对应的元素索引，取值为x_Az＝1,2,K,k，x_Rg＝1,2,K,k；a₀、a₁为方位向的多项式常数项和一次项系数；b₀、b₁为距离向的多项式常数项和一次项系数；

步骤S5.3，根据步骤S5.2获取的方位向和距离向的偏移量拟合系数，更新坐标映射数据1得到新的坐标映射数据2；更新的公式如下：

对坐标映射数据1的每一个像元(i,j)，得到其复数像元值(i_Value,1,j_Value,1)，根据方位向和距离向的偏移量，有

i_Value,2＝i_Value,1+a₀+a₁·i_Value,1

j_Value,2＝j_Value,1+b₀+b₁·j_Value,1

其中，i_Value,2、j_Value,2为更新后的坐标映射数据2的方位向和距离向像元值；i_Value,1、j_Value,1为坐标映射数据1的方位向和距离向像元值；a₀、a₁为方位向的多项式常数项和一次项系数；b₀、b₁为距离向的多项式常数项和一次项系数。

9.根据权利要求1所述的影像正射纠正方法，其特征在于，步骤S6具体包括以下子步骤：

步骤S6.1，计算ZY3数据在外接矩形数字高程模型RefDEM中的坐标(i_DEM,j_DEM)，计算公式如下所示：

i_DEM＝(B_i-B_ref,1)/B_ref,res

j_DEM＝(L_i-L_ref,1)/L_ref,res

其中，(B_ref,1,L_ref,1)为外接矩形数字高程模型RefDEM的起始点地理坐标；(B_i,L_i)为资源三号卫星影像控制点数据ZY3的经纬度地理坐标；B_ref,res为纬度方向外接矩形数字高程模型RefDEM分辨率，L_ref,res为经度方向外接矩形数字高程模型RefDEM分辨率；(i_DEM,j_DEM)为ZY3数据在外接矩形RefDEM中的坐标；

步骤S6.2，计算(i_DEM,j_DEM)在SAR模拟影像中的像素坐标(i_MSAR,j_MSAR)；

步骤S6.3，将所有ZY3数据重复步骤S6.1和S6.2，获得ZY3在SAR模拟影像中的像点坐标(i_MSAR,j_MSAR)，并基于所述像点坐标(i_MSAR,j_MSAR)和ZY3的地理坐标和高程构成ZY3控制点数据在多视SAR影像中的坐标集(B_i,L_i,H_i,i_MSAR,j_MSAR)；其中，(B_i,L_i,H_i)为ZY3的地理坐标和高程；(i_MSAR,j_MSAR)为ZY3在多视SAR影像中的像点坐标。

10.根据权利要求9所述的影像正射纠正方法，其特征在于，步骤S9具体包括以下子步骤：

步骤S9.1，根据原始SAR影像的四个角点像素坐标、原始SAR成像参数和输入的DOM数据的纬度方向分辨率、经度方向分辨率，确定待校正的DOM数据的宽和高，以及DOM起始点地理坐标；

步骤S9.2，根据步骤S9.1的地理参数，利用步骤S3.4构建的RD模型，逐像素对待校正的DOM影像进行计算，得到待校正的DOM像素在原始SAR影像中的像素坐标(x',y')；

步骤S9.3，对步骤S9.2得到的像素坐标(x',y')进行改正，得到改正后的原始SAR像素坐标(x,y)，获取原始SAR像素坐标(x,y)对应的像元值Pix_DOM(x,y)，将所述像元值Pix_DOM(x,y)保存到步骤S9.2所述的待校正的DOM像素坐标pix(x_DOM,y_DOM)中，遍历整个待校正的DOM像素，获得DOM正射纠正影像。

Images