CN103235304B - 一种基于误差等效rd模型修正的sar几何校正方法 - Google Patents

Abstract

本发明提供了一种基于误差等效RD模型修正的SAR几何校正方法，仅需很少的控制点，且对控制点的均匀分布特性要求不高；而且采用常规的数值求解方法即可得到很好的结果，无需迭代。第一步：建立基于误差等效RD模型修正的距离-多普勒方程，第二步：对第一步中得到的方程进行泰勒级数一级展开得到误差方程；第三步：采用最小二乘的方法求解第二步得到的误差方程；第四步：求得误差方程的解后，采用修正后的距离多普勒方程对每个像素进行几何定位，即实现了SAR的几何校正。

Description

一种基于误差等效RD模型修正的SAR几何校正方法

技术领域

本发明是一种基于误差等效RD模型修正的SAR几何校正方法，涉及合成孔径雷达（SAR）信号处理领域，具体涉及合成孔径雷达图像几何校正处理领域。

背景技术

图像几何校正是星载SAR数据处理的重要环节，只有通过几何校正，将SAR图像的每个像素赋予地理位置信息，SAR图像才能真正成为地理空间信息的载体，支撑后续的应用。SAR图像几何校正的精度对SAR图像应用效益的发挥有着至关重要的影响。SAR几何校正通常可分为无控校正（又称***级几何校正）和有控校正两种。其中***级几何校正是指基于轨道或航迹测量结果、SAR***参数以及成像处理参数对SAR图像进行几何校正，然而限于轨道测量精度、SAR***时间精度等，***级几何校正后的SAR图像定位精度有时无法满足应用的需求，尤其在我国目前***研制水平还不是很高的情况下，这种现象比较突出。为此，通常需要利用地面控制点，也即地理位置精确标定的点，来进行有控校正，提高SAR图像几何定位精度。

有控校正的方法有基于多项式模型的方法，基于共线模型的方法[Konecny G,Schuhr W.Reliability of radar image data.Proceedings of16_th ISPRS Congress,Kyoto,1998,1～10.]，以及基于SAR成像原理的距离-多普勒（RD）模型方法等。其中基于RD模型的方法符合成像几何故从原理上讲具有最佳的精度，因此成为目前SAR处理中的首选方法。现有技术中，基于RD模型的有控校正方法普遍通过精化轨道模型来完成几何精校正[尤红建，付琨，《合成孔径雷达图像精准处理》，北京:科学出版社，2011,ISBN987-7-03-031168-9]，即将定位误差的原因等效归结为由轨道位置、速度存在误差而引起，通常用SAR天线的初始位置（[X₀,Y₀,Z₀]^T），速度（[V_x,V_y,V_z]^T）和加速度（[a_x,a_y,a_z]^T）来对卫星轨道或机载航迹进行建模，并基于控制点建立距离-多普勒方程，方程如下

距离方程： $R_{\mathrm{near}}+N_i\mathrm{\δr}=\sqrt{{(X_{\mathrm{si}}-x_i)}^2+{(Y_{\mathrm{si}}-y_i)}^2+{(Z_{\mathrm{si}}-z_i)}^2}$

多普勒方程： $f_d=-\frac{2}{\mathrm{\λ}}\frac{V_{\mathrm{sxi}}(X_{\mathrm{si}}-x_i)+V_{\mathrm{syi}}(Y_{\mathrm{si}}-y_i)+V_{\mathrm{szi}}(Z_{\mathrm{si}}-z_i)}{R_{\mathrm{near}}+N_i\mathrm{\δr}}$

其中R_near表示图像距离向起始像素对应的斜距，即近距；δr为图像的距离向像素间隔；λ为波长；i＝1,...,K为控制点号，K为控制点总数；N_i为第_i个控制点的距离向像素号；[x_i,y_i,z_i]^T为第i个控制点的三维位置坐标；f_d为成像处理所采用的多普勒中心频率；其余参数如下，

$\left\{\begin{array}{c}{X}_{\mathrm{si}}=X_0+V_x(t_i-t_0)+\frac{1}{2}{a}_x{(t_i-t_0)}^2\\ Y_{\mathrm{si}}=Y_0+V_y(t_i-t_0)+\frac{1}{2}{a}_y{(t_i-t_0)}^2\\ Z_{\mathrm{si}}=Z_0+V_z(t_i-t_0)+\frac{1}{2}{a}_z{(t_i-t_0)}^2\\ V_{\mathrm{sxi}}=V_x=a_x(t_i-t_0)\\ V_{\mathrm{syi}}=V_y+a_y(t_i-t_0)\\ V_{\mathrm{szi}}=V_z+a_z(t_i-t_0)\end{array}\right.$

[X₀,Y₀,Z₀]^T为参考时刻t₀时的SAR天线位置，t_i为控制点i对应的方位时间。

基于上述方程，以现有轨道或航迹的测量值为初值，通过方程迭代求解来获得[X₀,Y₀,Z₀,V_x,V_y,V_z,a_x,a_y,a_z这12个参数的修正值，然后利用修正后的轨道对SAR图像进行重新定位，以提高SAR图像的定位精度。虽然上述技术方案从定位模型上符合SAR成像机理，而实际上将定位误差都归因于轨道位置速度误差其实也只是一种等效。实际中，由于***时间误差、大气传输延迟误差等的存在，距离-多普勒方程中的近距R_near、方位时间t_i等也可能存在误差。

因此说，现有技术方案的轨道修正模型也是对于SAR***定位误差源的一种等效模型。这种误差等效模型存在如下几个缺点：

1）所需的控制点较多且要求分布比较均匀。由于一个控制点可以建立2个方程，求解12个参数从理论上需要6个控制点。然而，实际中上述方程往往呈现病态，尤其是在星载SAR中，由于卫星离成像区域比较远，达几百甚至上千公里，而地面上一景图像则范围相对比较有限，通常在几十公里左右，因此各个控制点的误差方程存在很强的线性相关性，方程呈现明显的病态。故而，上述方程求解对控制点的要求非常高，一方面要求控制点较多且控制点自身误差较小；另一方面，要求控制点比较均匀的分布于场景中。工程实践表明，通常需要十个以上均匀分布的控制点以达到较为理想的定位结果。

2）病态方程的求解效率较低，且可能存在迭代不收敛的情况。如前所述，当SAR离场景较远（如星载SAR）的情况下，现有技术方案的方程往往呈现很严重的病态，为了解决病态方程的求解问题，人们提出了诸如岭估计、主成分估计、谱修正估计[王新洲，刘丁酉.最小二乘估计中法方程的迭代解法.湖北民族大学（自然科学版），2002,20(3):1～4]等估计方法，其中，谱修正估计方法是能较好处理星载SAR病态方程求解的一种迭代方法。然而，一方面，迭代会导致求解的效率下降，另一方面，即使是谱修正等比较好的求解方法，还是可能出现迭代不收敛的情况。

总之，现有技术方案的稳健性、鲁棒性还有待提高，对控制点数及均匀分布的依赖性也需进一步降低。

发明内容

针对现有技术方案的缺点，本发明提供了一种基于误差等效RD模型修正的SAR几何校正方法，仅需很少的控制点，且对控制点的均匀分布特性要求不高；而且采用常规的数值求解方法即可得到很好的结果，无需迭代。

本发明的技术方案如下：

该基于误差等效RD模型修正的SAR几何校正方法，包括以下步骤：

第一步：建立基于误差等效RD模型修正的距离-多普勒方程，如下：

距离方程： $F_{1i}(\mathrm{\ΔR},p)=R_{\mathrm{near}}+\mathrm{\ΔR}+p{N}_i+N_i\mathrm{\δr}-\sqrt{{(X_{\mathrm{si}}-x_i)}^2+{(Y_{\mathrm{si}}-y_i)}^2+{(Z_{\mathrm{si}}-z_i)}^2}=0;$

多普勒方程： $F_{2i}(\mathrm{\ΔR},p,\mathrm{\Δ}{f}_d,q)=f_d+\mathrm{\Δ}{f}_d+q{M}_i+\frac{2}{\mathrm{\λ}}+\frac{V_{\mathrm{sxi}}(X_{\mathrm{si}}-x_i)+V_{\mathrm{syi}}(Y_{\mathrm{si}}-y_i)+V_{\mathrm{szi}}(Z_{\mathrm{si}}-z_i)}{R_{\mathrm{near}}+\mathrm{\ΔR}+p{N}_i+N_i\mathrm{\δ}r}=0;$ 其中R_near表示图像距离向起始像素对应的斜距，即近距；δr为图像的距离向像素间隔；λ为波长；i＝1,...,K为控制点号，K为控制点总数；N_i为第i个控制点的距离向像素号；M_i是第i个控制点的方位向像素号；[x_i,y_i,z_i]^T为第i个控制点的三维位置坐标；f_d为成像处理所采用的多普勒中心频率；[X_si,Y_si,Z_si]^T为SAR卫星与第i个控制点相对位移引起的多普勒频率为f_d时，卫星的位置矢量；[V_sxi,V_syi,V_szi]^T为SAR卫星与第i个控制点相对位移引起的多普勒频率为f_d时，卫星与目标的相对速度矢量；ΔR,p,Δf_d,q为四个待定的误差等效模型参数；此处称ΔR+pN_i为距离修正量，ΔR和p为距离模型修正参数；称Δf_d+qM_i为多普勒修正量，Δf_d和q为多普勒模型修正参数；

第二步：对第一步中得到的方程进行泰勒级数一级展开得到误差方程；

第三步：采用最小二乘的方法求解第二步得到的误差方程；

第四步：求得误差方程的解后，采用修正后的距离多普勒方程对每个像素进行几何定位，即实现了SAR的几何校正。

其中第（1）步所述的距离修正量采用其中，L_p为距离修正量随像素距离位置变化的阶数，L_g为距离修正量随像素方位位置变化的阶数，为N_i的l次方,为M_i的l次方，p_l和g_l为待定参数，该式子表明距离修正量为距离向和方位向位置的多项式函数；或/和多普勒中心频率修正量采用为多普勒修正量随像素距离位置变化的阶数，L_h为多普勒修正量随像素方位位置变化的阶数，为N_i的l次方,为M_i的l次方，q_l和h_l为待定参数，该式子表明中心频率修正模型为距离向和方位向位置的多项式函数。

其中第（1）步所述的距离修正量采用中心频率修正量采用其中x_i为与距离像素号N_i呈线性关系的变量，y_i为与方位向像素号M_i呈线性关系的变量。

本发明的有益效果：

1、本发明仅需很少的控制点（采用两个参数时，理论上只需要1个控制点），且对控制点的均匀分布特性要求不高；同时采用常规的数值求解方法即可得到很好的结果，无需迭代，提高了运算效率。

2、本发明建立了误差等效修正模型的RD方程，引入距离修正量和多普勒中心频率修正量两类等效的模型修正参数，以提高几何定位精度。

附图说明

图1一种基于误差等效RD模型修正的SAR几何校正方法技术流程图。

具体实施方式

下面结合附图对本发明作进一步介绍。

如图1所示，本发明的基于误差等效RD模型修正的SAR几何校正方法，包括以下步骤：

第一步：建立基于误差等效RD模型修正的距离-多普勒方程，如下：

距离方程： $F_{1i}(\mathrm{\ΔR},p)=R_{\mathrm{near}}+\mathrm{\ΔR}+p{N}_i+N_i\mathrm{\δr}-\sqrt{{(X_{\mathrm{si}}-x_i)}^2+{(Y_{\mathrm{si}}-y_i)}^2+{(Z_{\mathrm{si}}-z_i)}^2}=0;$

多普勒方程： $F_{2i}(\mathrm{\ΔR},p,\mathrm{\Δ}{f}_d,q)=f_d+\mathrm{\Δ}{f}_d+q{M}_i+\frac{2}{\mathrm{\λ}}+\frac{V_{\mathrm{sxi}}(X_{\mathrm{si}}-x_i)+V_{\mathrm{syi}}(Y_{\mathrm{si}}-y_i)+V_{\mathrm{szi}}(Z_{\mathrm{si}}-z_i)}{R_{\mathrm{near}}+\mathrm{\ΔR}+p{N}_i+N_i\mathrm{\δ}r}=0;$ 其中R_near表示图像距离向起始像素对应的斜距，即近距；δr为图像的距离向像素间隔；λ为波长；i＝1,...,K为控制点号，K为控制点总数；N_i为第i个控制点的距离向像素号；M_i是第i个控制点的方位向像素号；[x_i,y_i,z_i]^T为第i个控制点的三维位置坐标；f_d为成像处理所采用的多普勒中心频率；[X_si,Y_si,Z_si]^T为SAR卫星与第i个控制点相对位移引起的多普勒频率为f_d时，卫星的位置矢量；[V_sxi,V_syi,V_szi]^T为SAR卫星与第i个控制点相对位移引起的多普勒频率为f_d时，卫星与目标的相对速度矢量；ΔR,p,Δf_d,q为四个待定的误差等效模型参数。此处称ΔR+pN_i为距离修正量，ΔR和p为距离模型修正参数；称Δf_d+qM_i为多普勒修正量，Δf_d和q为多普勒模型修正参数；

第二步：对上述方程进行泰勒级数一级展开，得到如下误差方程：

$\left[\begin{array}{cccc}\frac{\∂F_{11}}{\∂\mathrm{\ΔR}}& \frac{\∂F_{11}}{\∂p}& 0& 0\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ \frac{\∂F_{1K}}{\∂\mathrm{\ΔR}}& \frac{\∂F_{1K}}{\∂p}& 0& 0\\ \frac{\∂F_{21}}{\∂\mathrm{\ΔR}}& \frac{\∂F_{21}}{\∂p}& \frac{\∂F_{21}}{\∂\mathrm{\Δ}{f}_d}& \frac{\∂F_{21}}{\∂q}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ \frac{\∂F_{2K}}{\∂\mathrm{\ΔR}}& \frac{\∂F_{2K}}{\∂p}& \frac{\∂F_{2K}}{\∂\mathrm{\Δ}{f}_d}& \frac{\∂F_{2K}}{\∂q}\end{array}\right]\left[\begin{array}{c}\mathrm{\ΔR}\\ p\\ \mathrm{\Δ}{f}_d\\ q\end{array}\right]=\left[\begin{array}{c}-F_{11}\left(\mathrm{0,0}\right)\\ .\\ .\\ .\\ -F_{1K}\left(\mathrm{0,0}\right)\\ -F_{21}\left(\mathrm{0,0,0,0}\right)\\ .\\ .\\ .\\ -F_{2K}\left(\mathrm{0,0,0,0}\right)\end{array}\right]$

记上述方程为AX＝Y，其中

$\left\{\begin{array}{c}\frac{\∂F_{1i}}{\∂\mathrm{\ΔR}}=1\\ \frac{\∂F_{1i}}{\∂p}=N_i\\ \frac{\∂F_{2i}}{\∂\mathrm{\ΔR}}=\frac{2}{\mathrm{\λ}}\frac{(-1)[V_{\mathrm{sxi}}(X_{\mathrm{si}}-x_i)+V_{\mathrm{syi}}(Y_{\mathrm{si}}-y_i)+V_{\mathrm{szi}}(Z_{\mathrm{si}}-z_i)]}{{(R_{\mathrm{near}}+\mathrm{\ΔR}+p{N}_i+N_i\mathrm{\δr})}^2}{|}_{\mathrm{\ΔR}=0,p=0}\\ \frac{\∂F_{2i}}{\∂p}=\frac{2}{\mathrm{\λ}}\frac{(-N_i)[V_{\mathrm{sxi}}(X_{\mathrm{si}}-x_i)+V_{\mathrm{syi}}(Y_{\mathrm{si}}-y_i)+V_{\mathrm{szi}}(Z_{\mathrm{si}}-z_i)]}{{(R_{\mathrm{near}}+\mathrm{\ΔR}+p{N}_i+N_i\mathrm{\δr})}^2}{|}_{\mathrm{\ΔR}=0,p=0}\\ \frac{\∂F_{2i}}{\∂\mathrm{\Δ}{f}_d}=1\\ \frac{\∂F_{2i}}{\∂\mathrm{\Δq}}=M_i\end{array}\right.$

则对于K个点而言，有2K个方程，要求解4个位置数，理论上最少只需2个控制点。

第三步：采用最小二乘的方法求解上述方程，有X＝(A^TA)_-1A^TY，其中上标“^T”表示转置，“^-1”表示求逆。A^TA是满秩的，且通常条件数较小。

第四步：求得上述参数后，按照下面的方程对每个像素进行几何定位，定位求解方法与传统几何校正定位求解方法相同：

$\left\{\begin{array}{c}{R}_{\mathrm{near}}+\mathrm{\ΔR}+p{N}_j+N_j\mathrm{\δr}=\sqrt{{(X_{\mathrm{sj}}-x_t)}^2+{(Y_{\mathrm{sj}}-y_t)}^2+{(Z_{\mathrm{sj}}-z_t)}^2}\\ f_d+\mathrm{\Δ}{f}_d+q{M}_j=-\frac{2}{\mathrm{\λ}}\frac{V_{\mathrm{sxj}}(X_{\mathrm{sj}}-x_t)+V_{\mathrm{syj}}(Y_{\mathrm{sj}}-y_t)+V_{\mathrm{szj}}(Z_{\mathrm{sj}}-z_t)}{R_{\mathrm{near}}+\mathrm{\ΔR}+p{N}_j=N_j\mathrm{\δr}}\\ \frac{x_t^2+y_t^2}{{(R_e+h_t)}^2}+\frac{z_t^2}{R_{\mathrm{pt}}^2}=1,R_{\mathrm{Pt}}=(1-f)(R_e+h_t)\end{array}\right.$

其中[N_j,M_j]为距离向、方位向像素号，(x_t,y_t,z_t)^T为待求的目标位置，h_t为目标处高程；R_e为地球模型赤道半径，f为地球模型扁率因子，在WGS84地球模型下，R_e为6378137m，f约为0.003352。

需要说明的是，实际处理中可以根据需要采用2个误差参数（ΔR，Δf_d），或采用如下更多的参数：距离修正模型：其中L_p为距离修正模型的阶数。中心频率修正模型：其中L_q为中心频率修正模型的阶数。模型阶数的确定通常根据无控条件下的定位误差分布情况而定，如果定位误差比较一致，即畸变很小，则采用2个误差参数进行修正即可；如定位误差呈现沿距离向或/和方位向的变化，则根据变化情况采用合适的高阶修正模型。

如距离修正模型参数采用即距离修正模型为距离向和方位向位置的多项式函数，或/和多普勒中心频率修正模型参数采用即中心频率修正模型为距离向和方位向位置的多项式函数，也属于本发明的范畴之内。

如距离修正模型参数采用中心频率修正模型采用其中x_i为与距离像素号N_i呈线性关系的变量（如斜距等），y_i为与方位向像素号M_i呈线性关系的变量（如方位向距离等），则也属于本发明的范畴之内。

下面通过实际数据的处理实例来验证本发明的优点。我们在北京地区用差分GPS实测了22个点，下面给出利用这些点作控制点（GCP）或检查点的北京地区2幅星载SAR图像几何校正后的定位误差。表1为现有技术方案的校正结果，可见，采用现有技术方案，误差方程求解矩阵的条件数非常大，因此，当GCP点分布不均匀时，缺少控制点的区域定位误差明显增大，最终导致图像总体的位置中误差非常大。表2给出了本发明的几何校正结果，首先由于这两幅图像在几何校正前，其位置误差的方差就比较小，因此，选用距离修正量和多普勒修正量的一阶形式，即只采用ΔR和Δf_d这两个待定参数的模型修正方法进行修正，从校正结果可见，误差方程的求解矩阵条件数接近1，因此，方程求解是非病态的，可见，该方法是稳健的；并且其校正精度与现有技术方案控制点较多且分布均匀时的结果相当，甚至更佳，验证了本发明方法的优越性。

表1现有技术方案几何校正结果

表2本发明方法几何校正结果（采用2个参数修正模型）

Claims (3)

1.基于误差等效RD模型修正的SAR几何校正方法，其特征在于，包括以下步骤：

第一步：建立基于误差等效RD模型修正的距离-多普勒方程，如下：

距离方程： $F_{1i}(\mathrm{\ΔR},p)=R_{\mathrm{near}}+\mathrm{\ΔR}+p{N}_i+N_i\mathrm{\δr}-\sqrt{{(X_{\mathrm{si}}-x_i)}^2+{(Y_{\mathrm{si}}-y_i)}^2+{(Z_{\mathrm{si}}-z_i)}^2}=0;$

多普勒方程： $F_{2i}(\mathrm{\ΔR},p,\mathrm{\Δ}{f}_d,q)=f_d+\mathrm{\Δ}{f}_d+q{M}_i+\frac{2}{\mathrm{\λ}}\frac{V_{\mathrm{sxi}}(X_{\mathrm{si}}-x_i)+V_{\mathrm{syi}}(Y_{\mathrm{si}}-y_i)+V_{\mathrm{szi}}(Z_{\mathrm{si}}-z_i)}{R_{\mathrm{near}}+\mathrm{\ΔR}+p{N}_i+N_i\mathrm{\δr}}=0;$

其中R_near表示图像距离向起始像素对应的斜距，即近距；δr为图像的距离向像素间隔；λ为波长；i＝1,...,K为控制点号，K为控制点总数；N_i为第i个控制点的距离向像素号；M_i是第i个控制点的方位向像素号；[x_i,y_i,z_i]^T为第i个控制点的三维位置坐标；f_d为成像处理所采用的多普勒中心频率；[X_si,Y_si,Z_si]^T为SAR卫星与第i个控制点相对位移引起的多普勒频率为f_d时，卫星的位置矢量；[V_sxi,V_syi,V_szi]^T为SAR卫星与第i个控制点相对位移引起的多普勒频率为f_d时，卫星与目标的相对速度矢量；ΔR,p,Δf_d,q为四个待定的误差等效模型参数；此处称ΔR+pN_i为距离修正量，ΔR和p为距离模型修正参数；称Δf_d+qM_i为多普勒中心频率修正量，Δf_d和q为多普勒模型修正参数；

第二步：对第一步中得到的方程进行泰勒级数一级展开得到误差方程；

第三步：采用最小二乘的方法求解第二步得到的误差方程；

第四步：求得误差方程的解后，采用修正后的距离多普勒方程对每个像素进行几何定位，即实现了SAR的几何校正。

2.如权利要求1所述的基于误差等效RD模型修正的SAR几何校正方法，其特征在于，其中第一步所述的距离修正量采用其中，L_p为距离修正量随像素距离位置变化的阶数，L_g为距离修正量随像素方位位置变化的阶数，p_l和g_l为待定参数，该式子表明距离修正量为距离向和方位向位置的多项式函数；或/和多普勒中心频率修正量采用L_q为多普勒中心频率修正量随像素距离位置变化的阶数，L_h为多普勒中心频率修正量随像素方位位置变化的阶数，q_l和h_l为待定参数，该式子表明多普勒模型为距离向和方位向位置的多项式函数，所述为N_i的l次方,为M_i的l次方。

3.如权利要求1所述的基于误差等效RD模型修正的SAR几何校正方法，其特征在于，其中第一步所述的距离修正量采用多普勒中心频率修正量采用其中x_i为与距离向像素号N_i呈线性关系的变量，y_i为与方位向像素号M_i呈线性关系的变量。