CN102654576B - 基于sar图像和dem数据的图像配准方法 - Google Patents

Abstract

本发明公开了一种基于SAR图像和DEM数据的图像配准方法，主要解决现有技术对SAR图像中地面目标点的定位精度不高和在实际工程中实现难的问题。其实现过程为：根据距离-多普勒RD原理及地球模型方程，用最速下降法求解出地面目标点在球心坐标下的经度L，纬度B和高程值h′，并将其值转换到大地坐标系；查找数字高程模型DEM数据库中与B、L对应的高程值h″，将|h′-h″|与设定的阈值ε比较，若|h′-h″|>ε，则令h=h″，返回重新解方程组，若|h′-h″|≤ε，停止迭代，则h′就是地面目标点的高程值，所对应的纬度B和经度L为地面目标点的经纬度；对经纬度进行线性修正消除系统误差，得到最终的地面目标点定位结果。本发明具有运算量小、定位精度高的优点，可用于实际工程中对SAR图像中地面目标点的定位。

Description

基于SAR图像和DEM数据的图像配准方法

技术领域

本发明属于图像处理技术领域，涉及图像配准，特别是一种地面目标定位和配准方法，可用于SAR图像中的地面目标定位和SAR影像的几何校正处理过程的预处理环节。

背景技术

合成孔径雷达SAR技术以其高分辨率、全天候、全天时的工作能力以及对目标的识别特性，在军事、地形测绘、海洋监测等方面发挥着日益重要的作用。图像配准是将来自同一区域、在相同时间或不同时间、不同传感器或不同视点所获取的遥感数字图像在空间上进行对准。利用SAR图像和数字高程模型DEM数据库的图像配准，可以对图像的像素进行精确的地理定位，SAR图像和DEM数据的图像配准过程就是对SAR图像中目标进行定位的过程。

在对SAR图像中目标进行定位时，目前有3种定位方法：基于多项式模型的定位方法、基于共线方程模型的定位方法和基于距离-多普勒RD原理的定位方法。

一．基于多项式模型的定位方法。这种方法是通过构建多项式模型对目标进行定位，多项式模型用一组多项式对输入空间和输出空间的变换进行简单的函数逼近，是一种基于平面的变换。通过在地形图上和卫星遥感影像上选取参考点并读取其地理坐标和影像坐标，用最小二乘法建立两者的多项式变换关系，然后将图像坐标系中的象素以此关系转换到地图坐标系中从而实现对地定位。

多项式模型定位方法不依赖于卫星传感器成像系统本身，它要求选择参考点，以定出一系列多项式系数，因此与卫星轨道、传感器姿态等参数全然无关。该方法操作简单且原理清晰，但在实际应用上具有明显的局限性，特别是在宽阔的洋面上，或是在一望无际的沙漠里，或是在有大片浓云覆盖的地区，很难找到用于图像定位的可靠的参考点，因此也就很难谈得上有多少定位精确度。

二．基于共线方程模型的定位方法。这种方法是通过构建共线方程模型进行目标定位，共线方程模型的基本思想是将侧视SAR图像看成是一个中心投影遥感器获取的等效图像，然后借用光学遥感的中心投影方式进行校正，从而把SAR的距离投影和侧视几何成像转化为中心投影的透视关系。该方法是建立在影像坐标和地面坐标严格变换关系的基础上的，是对成像空间几何形态的直接描述，因而，从理论上讲比多项式方法更严密些。

对于诸如Landsat、SPOT等卫星的TM光学成像机理而言，共线方程模型无疑在理论上很好地反映了中心投影的透视几何关系。但是对于SAR成像而言，虽然也可以将距离投影转换为中心投影，仿照摄影测量中惯用的共线方程来表达，但共线方程只是一个虚构的中心投影构像模型，并不能真实地反映出SAR距离投影的几何形态，即位于传感器所发出之探测波的同一波前球面上的所有地物点将成像于同一点。因此SAR不适用共线方程的概念。

三．基于距离-多普勒RD原理的定位方法。这种方法是1982年首先由Curlander等人提出的，此方法利用地球物理模型方程、SAR多普勒方程、SAR斜距方程对图像像素进行定位，并于1989年完成了SAR图像自动校正地理编码后处理系统。该系统可以对SAR图像进行精确定位，精度可以达到±200m。

距离-多普勒RD原理主要是根据合成孔径雷达SAR的成像原理，这样点目标的空间定位由两方面来决定：一是根据目标回波时间的长短来确定点目标到雷达的距离，二是根据目标回波的多普勒特性。如图1所示，在SAR的照射区域内，分布着等时延的同心圆束和等多普勒频移的双曲线束，图中S是卫星位置，D是星下点，P是地面目标点，O是地球中心。点目标必定在等时延的同心圆束和一个等多普勒曲面的交线上，如果预先知道点目标所在的约束平面，如地表面，则像点的空间位置可以唯一确定。根据距离向上回波信号的时延和方位向上的多普勒频移这两个信息来确定点目标位置的方法，就是距离-多普勒RD定位方法。

根据RD原理，列出RD方程组如下：

1)多普勒频率方程：

$f_{\mathrm{DC}}=-\frac{2}{\mathrm{\λR}}V\·R=-\frac{2}{\mathrm{\λR}}(V_s-V_t)\·(R_s-R_t),$

式中，f_DC为多普勒频率，λ为雷达波长，R是卫星到地面的距离，V表示速度矢量，R表示位置矢量，v_s、v_t分别指卫星平台和地面目标的速度矢量，且V_t=ω_e×R_t，ω_e是地球自转角速度。

2)斜距方程：

R=|R_s-R_t|=[(R_s-R_t)·(R_s-R_t)]^1/2，

其中，R_s和R_t分别为卫星平台和地面目标的位置矢量，它们可通过卫星下行数据影像头文件中的星历参数并通过函数预测得到，设R_s=(x_s,y_s,z_s)^T，R_t=(x_t,y_t,z_t)^T。

3)地球模型方程：

$\frac{(x_t^2+y_t^2)}{{(R_e+h)}^2}+\frac{z_t^2}{R_p^2}=1,$

式中，R_e为平均赤道半径，R_p为极轴半径，且R_p=[(1-f)(R_e+h)]²，平坦度因子f=1/298.255，(x_t,y_t,z_t)为地面目标点的位置矢量，h是SAR图像中目标所在区域的平均高程值。

这种方法虽然能解决基于多项式模型的定位方法在实际工程应用中局限性大且定位精度低和基于共线方程模型的定位方法不能很好地用在SAR图像定位上的缺陷。但对这种定位方法的求解，大多采用最小二乘迭代算法，求解过程复杂，且求解过程中目标点的高度h未知，需要用地面目标位置矢量(x_t,y_t,z_t)来表示再通过计算求出高度h。同时要对计算出的地面目标位置矢量进行坐标转换，这样在求解过程中，由于已知的数据如地球赤道半径、地球极半径等为近似数据，故得到的计算结果误差就会被放大使定位精度降低，而且无法验证所得数据的合理性。

发明内容

本发明的目的在于针对上述已有方法中的不足，提出了一种基于SAR图像和DEM数据的图像配准方法，以降低定位时求解方程的复杂性，提高定位精度，满足实际工程需要。

实现本发明目的的技术思路是：通过将SAR图像头文件中给出的所需参数代入距离-多普勒RD原理定位方法中的3个联立方程，用最速下降法及将SAR图像和DEM数据相结合的方法，求解方程得到目标点的经纬度值和高度值，对所得数据进行线性修正得到更准确、高精度的地面点目标位置。其具体步骤包括如下：

（1）查找SAR图像的所带头文件，得到地面目标点所在区域的平均高程值，作为高程值h的初始值；

（2）将初始值h分别代入由距离-多普勒RD原理确定的斜距方程、多普勒频率方程及地球模型方程中，将这3个方程联立，用最速下降法求解出地面目标点的位置矢量R_t=(x_t,y_t,z_t)，其中，x_t、y_t、z_t分别是球心坐标系中地面目标点的X轴坐标、Y轴坐标、Z轴坐标；

（3）将计算得到的位置矢量R_t=(x_t,y_t,z_t)转换为大地坐标系中对应的纬度B、经度L和高程值h′；

（4）根据纬度B和经度L查找数字高程模型DEM数据库中对应的高程值h″；

（5）将数字高程模型DEM数据库中对应的高程值h″与计算得到的高程值h′做差并求绝对值，再与所设定的阈值ε=0.3作比较，如果该绝对值大于阈值ε，则令h=h″，返回步骤（2），如果该绝对值小于等于阈值ε，则停止迭代，计算得到的高程值h′就是地面目标点的高程值，所对应的纬度B和经度L就是地面目标点的经纬度；

（6）对所得地面目标点的纬度B和经度L进行线性修正，得到修正后的纬度B′和经度L′，以消除系统误差：

$B^{\′}=B-0.001*(\frac{r}{2}-|\frac{r}{2}-i\left|\right)/\frac{r}{2}+0.002*(j-\frac{c}{2})/\frac{c}{2},$

$L^{\′}=L+0.001*(i-1)/r+0.001*(\frac{c}{2}-j)/\frac{c}{2}+0.0065*(\frac{r}{2}-|i-1\left|\right)/\frac{r}{2},$

其中，r是SAR图像的总行数，c是SAR图像的总列数，i是地面目标点所在行号，j是地面目标点所在列号；

（7）根据修正后的经度和纬度，得到地面目标点在大地坐标系中的位置。

本发明与现有技术相比具有以下优点：

1）本发明由于在计算中采用最速下降法，而没有采用通常所用的最小二乘迭代法，大大减少了运算量，并且提高了计算结果的精度。

2）本发明由于将SAR图像和DEM数据结合起来对地面目标进行定位，不仅对目标的定位结果更加精确，而且有利于在实际工程上的使用。

3）本发明由于对所得地面目标点的位置矢量进行了线性修正，提高了定位精度。

实验结果表明，本发明对地面目标点的定位可以满足实际工程中要求的。

附图说明

图1是本发明中的距离-多普勒RD定位原理图；

图2是本发明的实现流程图。

具体实施方式

参照图2，本发明的具体实施步骤如下：

步骤1.从SAR图像所带头文件提供的参数中，查获SAR图像中地面目标点所在区域的平均高程值h。

步骤2.联立由距离-多普勒RD原理确定的斜距方程、多普勒频率方程及地球模型方程，将这3个方程作为一个非线性方程组，用最速下降法求解得到地面目标点在球心坐标中的位置。

2a）根据距离-多普勒RD原理确定的斜距方程、多普勒频率方程和地球模型方程：

斜距方程：R′=|R_s-R_t|=[(R_s-R_t)·(R_s-R_t)]^1/2，

多普勒频率方程： $f_{\mathrm{DC}}=-\frac{2}{\mathrm{\λ}{R}^{\′}}V\·R=-\frac{2}{\mathrm{\λ}{R}^{\′}}(V_s-V_t)\·(R_s-R_t),$

地球模型方程： $\frac{(x_t^2+y_t^2)}{{(R_e+h)}^2}+\frac{z_t^2}{R_p^2}=1,$

其中，R′是斜距，R_s、R_t分别是卫星平台和地面目标点的位置矢量，λ是雷达波长，V是速度矢量，R是位置矢量，v_s、v_t分别是卫星平台和地面目标点的速度矢量，v_t=ω_e×R_t，ω_e是地球自转角速度，x_t、y_t、z_t分别是球心坐标系中地面目标点的X轴坐标、Y轴坐标、Z轴坐标，R_e是平均赤道半径，R_p是极轴半径且R_p=[(1-f)(R_e+h)]²，平坦度因子f=1/298.255，h是SAR图像中地面目标点所在区域的平均高程值；

2b）将这3个方程联立一个非线性方程组为：

$\left\{\begin{array}{c}{f}_1=R^{\′}={[(R_s-R_t)\·(R_s-R_t)]}^{1/2}\\ f_2=f_{\mathrm{DC}}=-\frac{2}{\mathrm{\λ}{R}^{\′}}(V_s-V_t)\·(R_s-R_t)\\ f_3=(x_t^2+y_t^2)/{(R_e+h)}^2+z_t^2/R_p^2-1\end{array}\right.;$

2c）将步骤1中查到的平均高程值h作为高程的初始值，代入非线性方程组中；

2d）给定一组值(x₀,y₀,z₀)，作为步骤2b）中非线性方程组解的初始值，代入非线性方程组得到f₁、f₂和f₃的值，将所得值代入目标函数F(x₀,y₀,z₀)，计算目标函数值：

$F=F(x_0,y_0,z_0)=\underset{i=1}{\overset{3}{\mathrm{\Σ}}}{f}_i^2,$

其中，F是目标函数值；

2e）若目标函数值F<ε，则(x₀,y₀,z₀)^T就是非线性方程组的一组实根，停止计算，否则继续执行下面的计算步骤，其中，ε是设定的阈值，其值为0.3，T表示转置；

2f）计算目标函数F(x₀,y₀,z₀)在点(x₀,y₀,z₀)处的偏导数：

$\frac{\&PartialD;F}{\&PartialD;x_0}=2\underset{j=1}{\overset{3}{\mathrm{\&Sigma;}}}{f}_j^2\&CenterDot;\frac{\&PartialD;f_j}{\&PartialD;x_0},$

$\frac{\&PartialD;F}{\&PartialD;y_0}=2\underset{j=1}{\overset{3}{\mathrm{\&Sigma;}}}{f}_j^2\&CenterDot;\frac{\&PartialD;f_j}{\&PartialD;y_0},$

$\frac{\&PartialD;F}{\&PartialD;z_0}=2\underset{j=1}{\overset{3}{\mathrm{\&Sigma;}}}{f}_j^2\&CenterDot;\frac{\&PartialD;f_j}{\&PartialD;z_0},$

其中，表示函数A对变量a求偏导；

2g）计算中间变量D的值：

$D={\left(\frac{\&PartialD;F}{\&PartialD;x_0}\right)}^2+{\left(\frac{\&PartialD;F}{\&PartialD;y_0}\right)}^2+{\left(\frac{\&PartialD;F}{\&PartialD;z_0}\right)}^2;$

2h）计算地面目标点在(x_t,y_t,z_t)处的值：

$x_t=x_0-\mathrm{\&lambda;}\frac{\&PartialD;F}{\&PartialD;x_0},$

$y_t=y_0-\mathrm{\&lambda;}\frac{\&PartialD;F}{\&PartialD;y_0},$

$z_t=z_0-\mathrm{\&lambda;}\frac{\&PartialD;F}{\&PartialD;z_0},$

其中，λ=F/D；

2i）重复步骤2d）~2h），直到满足x_t、y_t和z_t的精度值0.0201，得到SAR图像中地面目标点在球心坐标系中的位置坐标(x_t,y_t,z_t)。

步骤3.将计算得到的地面目标点位置坐标(x_t,y_t,z_t)转换到大地坐标系下。

3a）将大地坐标系中纬度B，经度L和高程h′转换到球心坐标系下地面目标点的位置坐标(x_t,y_t,z_t)：

$\left\{\begin{array}{c}{x}_t=(N+h^{\&prime;})\cos B\cos L\\ y_t=(N+h^{\&prime;})\cos B\sin L\\ z_t=(N(1-e^2)+h^{\&prime;})\sin B\end{array}\right.,$

其中，卯酉曲率半径

e为地球的第一偏心率，

a是地球赤道半径，b是地球极半径；

3b）将球心坐标系中地面目标点的位置坐标(x_t,y_t,z_t)代入步骤3a）中的大地坐标系中纬度B，经度L和高程h′公式，得到其在大地坐标系中纬度B、经度L和高程值h′的计算公式：

$B=\arctan (\frac{z_t+{\left(e^{\&prime;}\right)}^{2}b{\sin }^3\mathrm{\&theta;}}{P-e^{2}a{\cos }^3\mathrm{\&theta;}},$

$L=\arctan \left(\frac{y_t}{x_t}\right),$

$h^{\&prime;}=\frac{P}{\cos B}-N,$

其中，P是球心坐标系中地面目标点到Z轴的距离且

θ是中间变量且

e′为地球的第二偏心率，

3c）将步骤2i）得到的地面目标点位置坐标(x_t,y_t,z_t)代入步骤3b）中的纬度B、经度L和高程值h′计算公式中，得到其对应的纬度B、经度L和高程值h′。

步骤4.根据大地坐标系中的纬度B和经度L，查找数字高程模型DEM数据库中和此经纬度值所对应的高程值h″。

步骤5.比较数字高程模型DEM数据库中对应的高程值h″与计算得到的高程值h′差的绝对值与阈值ε的大小，以确定计算得到的高程值h′就是地面目标点的高程值。

5a）取阈值ε=0.3，如果|h′-h″>ε，则令h=h″，得到大地坐标系中的坐标(B，L，h″)，将其转换到球心坐标系，作为方程的初始值(x₀,y₀,z₀)，返回步骤2d）重新执行；

5b）如果|h′-h″|≤ε，则停止迭代，计算得到的高程值h′就是地面目标点的高程值，所对应的纬度B和经度L就是地面目标点的经纬度。

步骤6.对所得地面目标点的纬度B和经度L进行线性修正，以消除系统误差。

6a）将经纬度代入以下修正公式：

$B^{\&prime;}=B-0.001*(\frac{r}{2}-|\frac{r}{2}-i\left|\right)/\frac{r}{2}+0.002*(j-\frac{c}{2})/\frac{c}{2},$

$L^{\&prime;}=L+0.001*(i-1)/r+0.001*(\frac{c}{2}-j)/\frac{c}{2}+0.0065*(\frac{r}{2}-|i-1\left|\right)/\frac{r}{2},$

其中，r是SAR图像的总行数，c是SAR图像的总列数，i是地面目标点所在行号，j是地面目标点所在列号；

6b）根据以上的修正公式，得到修正后地面目标点的纬度B′和经度L′。

步骤7.由修正后的经度、纬度及计算得到的高程值，最终得到地面目标点在大地坐标系中的坐标(B′,L′,h′)。

本发明的效果可以通过如下仿真实验进一步说明。

1）仿真条件：

硬件平台为：Intel Core2Duo CPU E65502.33GHZ、2GB RAM

软件平台为：VC++6.0

2）仿真内容：

应用本发明分别对SAR图像头文件中给定的5个标准参考点和在SAR图像中取的20个地面目标点进行实验，在球心坐标系将计算值与Goole Earth中对应的各点值进行比较，来评价本发明的精度和实际应用性。

仿真1，应用本发明对SAR图像提供头文件中给定的5个标准参考点进行仿真实验，在球心坐标系下将计算值与GooleEarth中对应的各点值进行比较，实验结果如表1所示。

表15个标准参考点的实验结果

由表1可见，本发明中对5个标准参考点的定位精度在50m范围之内，查找每个点的时间为0.015s，说明本发明能满足实际工程中的定位精度需要和实时性要求。

仿真2，应用本发明在SAR图像中取20个点进行仿真实验，在球心坐标系下将计算值与Goole Earth中对应的各点值进行比较，实验结果如表2所示。

表2  20个地面目标点的实验结果

由表2可见，本发明中20个地面目标点的定位误差均控制在50m范围之内，查找每一个点的计算时间为0.015s，说明本发明能满足实际工程中的定位精度需要和实时性要求。

Claims (2)

1.一种基于SAR图像和DEM数据的图像配准方法，包括如下步骤：

（1）查找SAR图像的所带头文件，得到地面目标点所在区域的平均高程值，作为高程值h的初始值；

（2）将初始值h分别代入由距离-多普勒RD原理确定的斜距方程、多普勒频率方程及地球模型方程中，将这3个方程联立，用最速下降法求解出地面目标点的位置矢量R_t＝(x_t,y_t,z_t)，其中，x_t、y_t、z_t分别是球心坐标系中地面目标点的X轴坐标、Y轴坐标、Z轴坐标；

（3）将计算得到的位置矢量R_t＝(x_t,y_t,z_t)转换为大地坐标系中对应的纬度B、经度L和高程值h'；

（4）根据纬度B和经度L查找数字高程模型DEM数据库中对应的高程值h′′；

（5）将数字高程模型DEM数据库中对应的高程值h′′与计算得到的高程值h'做差并求绝对值，再与所设定的阈值ε＝0.3作比较，如果该绝对值大于阈值ε，则令h＝h''，返回步骤（2），如果该绝对值小于等于阈值ε，则停止迭代，计算得到的高程值h'就是地面目标点的高程值，所对应的纬度B和经度L就是地面目标点的经纬度；

（6）对所得地面目标点的纬度B和经度L进行线性修正，得到修正后的纬度B′和经度L′，以消除系统误差：

$B^{\&prime;}=B-0.001*(\frac{r}{2}-|\frac{r}{2}-i\left|\right)/\frac{r}{2}+0.002*(j-\frac{c}{2})/\frac{c}{2},$

$L^{\&prime;}=L+0.001*(i-1)/r+0.001*(\frac{c}{2}-j)/\frac{c}{2}+0.0065*(\frac{r}{2}-|i-1\left|\right)/\frac{r}{2},$

其中，r是SAR图像的总行数，c是SAR图像的总列数，i是地面目标点所在行号，j是地面目标点所在列号；

（7）根据修正后的经度和纬度，得到地面目标点在大地坐标系中的位置；

所述步骤（2）中用最速下降法解方程组，按如下步骤进行：

2a）确定非线性方程组为：

f₁＝R′＝|R_s-R_t|＝[(R_s-R_t)·(R_s-R_t)]^1/2，

$f_2=\frac{2}{{\mathrm{\&lambda;R}}^{\&prime;}}V\&CenterDot;R=-\frac{2}{{\mathrm{\&lambda;R}}^{\&prime;}}(V_s-V_t)\&CenterDot;(R_s-R_t),$

$f_3=\frac{(x_t^2+y_t^2)}{{(R_e+h)}^2}+\frac{z_t^2}{R_p^2}-1,$

其中，f₁、f₂分别是由RD原理确定的斜距方程和多普勒频率方程，f₃是地球模型方程，R′是斜距，R_s、R_t分别是卫星平台和地面目标点的位置矢量，λ是雷达波长，V是速度矢量，R是位置矢量，V_s、V_t分别是卫星平台和地面目标点的速度矢量，V_t＝ω_e×R_t，ω_e是地球自转角速度，R_e是平均赤道半径，R_p是极轴半径且R_p＝[(1-f)(R_e+h)]²，平坦度因子f＝1/298.255，h是SAR图像中地面目标点所在区域的平均高程值；

2b）给定一组值(x₀,y₀,z₀)，作为步骤2a）中非线性方程组解的初始值，分别得到f₁、f₂和f₃的值，将所得值代入目标函数F(x₀,y₀,z₀)，计算目标函数值：

$F=F(x_0,y_0,z_0)=\underset{i=1}{\overset{3}{\mathrm{\&Sigma;}}}{f_i}^2,$

其中，F是目标函数值；

2c）若目标函数值F＜ε，则(x₀,y₀,z₀)^T为非线性方程组的一组实根，停止计算，否则继续执行下面的计算步骤，其中，ε是设定的阈值，其值为0.3；

2d）计算目标函数F(x₀,y₀,z₀)在点(x₀,y₀,z₀)处的偏导数：

$\frac{\&PartialD;F}{{\&PartialD;x}_0}=2\underset{j=1}{\overset{3}{\mathrm{\&Sigma;}}}{f}_j^2\&CenterDot;\frac{{\&PartialD;f}_j}{{\&PartialD;x}_0},$

$\frac{\&PartialD;F}{{\&PartialD;y}_0}=2\underset{j=1}{\overset{3}{\mathrm{\&Sigma;}}}{f}_j^2\&CenterDot;\frac{{\&PartialD;f}_j}{{\&PartialD;y}_0},$

$\frac{\&PartialD;F}{{\&PartialD;z}_0}=2\underset{j=1}{\overset{3}{\mathrm{\&Sigma;}}}{f}_j^2\&CenterDot;\frac{{\&PartialD;f}_j}{{\&PartialD;z}_0},$

其中，

表示函数A对变量a求偏导；

2e）计算中间变量D的值：

$D={\left(\frac{\&PartialD;F}{{\&PartialD;x}_0}\right)}^2+{\left(\frac{\&PartialD;F}{{\&PartialD;y}_0}\right)}^2+{\left(\frac{\&PartialD;F}{{\&PartialD;z}_0}\right)}^2;$

2f）计算地面目标点在位置(x_t,y_t,z_t)处的值：

$x_t=x_0-\mathrm{\&lambda;}\frac{\&PartialD;F}{{\&PartialD;x}_0},$

$y_t=y_0-\mathrm{\&lambda;}\frac{\&PartialD;F}{{\&PartialD;y}_0},$

$z_t=z_0-\mathrm{\&lambda;}\frac{\&PartialD;F}{{\&PartialD;z}_0},$

其中，λ＝F/D，x_t、y_t、z_t分别是球心坐标系中地面目标点的X轴坐标、Y轴坐标、Z轴坐标；

2g）重复步骤2b）～2f），直到满足x_t、y_t和z_t的精度为0.0201，得到的(x_t,y_t,z_t)就是地面目标点的球心坐标。

2.根据权利要求1所述的基于SAR图像和DEM数据的图像配准方法，其中步骤（3）所述的坐标系转换，按如下步骤进行：

3a）将大地坐标系中的BLh′转换到球心坐标系中的x_ty_tz_t：

$\{\begin{array}{c}{x}_t=(N+h^{\&prime;})\cos B\cos L\\ y_t=(N+h^{\&prime;})\cos B\sin L\\ z_t=(N(1-e^2)+h^{\&prime;})\sin B\end{array},$

其中，卯酉曲率半径

e为地球的第一偏心率且

a是地球赤道半径，b是地球极半径，B是纬度，L是经度，h'是高程值；

3b）将球心坐标系中的(x_t,y_t,z_t)代入步骤3a）中的公式，得到在大地坐标系中纬度B、经度L和高程值h'的计算公式：

$B=\arctan (\frac{z_t+{\left(e^{\&prime;}\right)}^2{b\sin }^3\mathrm{\&theta;}}{P-e^2{a\cos }^3\mathrm{\&theta;}},$

$L=\arctan \left(\frac{y_t}{x_t}\right),$

$h^{\&prime;}=\frac{P}{\cos B}-N,$

其中，P是球心坐标系中地面目标点到Z轴的距离且

θ是中间变量且

e′为地球的第二偏心率且

x_t、y_t、z_t分别是球心坐标系中地面目标点的X轴坐标、Y轴坐标、Z轴坐标；

3c）将步骤2g）得到的(x_t,y_t,z_t)代入步骤3b）中的B、L、h′公式中，得到地面目标点的纬度B和经度L及高程值h'。

Images