CN101339244A - 一种机载sar图像自动目标定位方法 - Google Patents

Abstract

一种机载SAR图像自动目标定位方法，利用载机的飞行信息(位置、速度等)和SAR图像信息(分辨率、中心距离、中心多谱勒频率等)，将目标在SAR图像中的坐标转换为目标在地球上的实际地理坐标。本发明在距离－多谱勒方程的基础上，结合地球模型方程进行目标定位，既可实现无控制点的绝对定位，又考虑了地球模型对定位精度的影响；另一方面，利用泰勒级数展开对非线性方程组线性化，通过迭代的方法求解方程，提高了定位速度；同时，根据雷达工作方式适当选择迭代计算的初值点，控制了迭代中收敛的方向，避免了目标混淆。作为一种无控制点的定位方法，本发明具有精度高、速度快的优点，可进行实时处理，适合各种机载SAR系统。

Description

一种机载SAR图像自动目标定位方法

技术领域

本发明涉及一种机载SAR图像自动目标定位方法，利用载机的飞行信息和SAR图像信息将目标在SAR图像中的坐标转换为目标在地球上的实际地理坐标，适合于在机载SAR系统中对目标实行精确定位。

背景技术

合成孔径雷达(Synthetic Aperture Radar，SAR)本身是一种利用微波感知的主动式传感器，它不受天气、光照等条件限制，可对感兴趣目标进行全天候、全天时的侦察；此外，由于微波具有一定的穿透探测能力，可以探测到一定深度地表下的或是其他经过伪装或掩盖的目标，这使得SAR在目标探测方面具有巨大的应用潜力。

通常情况下，目标检测的结果得到的是目标在SAR图像中的位置。而在实际的应用中，不管是在军事方面，如常规境外侦察定位，飞行器导航、火控、目标跟踪等；或是民用方面，如自然灾害(森林火灾)的评估和预测、船只监测和海难救援等，仅得到目标在图像中的位置是不够的，这时更加重要的是得到目标在地球上的位置。因此，将目标在图像中的坐标转换为目标在地球上的地理坐标的SAR图像自动目标定位技术具有重要的研究意义。

传统的SAR图像目标定位方法一般是进行有参考点定位，即在SAR成像范围内找出一些位置已经确知的参考点(人造的或自然的)，然后再去确定图像中其它目标点相对于参考点的位置；这类方法依赖于较高精度的地面控制点，并且需要通过人工判读的方法确定控制点在影像上的位置，因此很难实现实时处理。目前，关于机载SAR目标定位的专门性文章并不多，在常见的目标定位算法中，多数是忽略地球模型对定位精度的影响，直接利用机载SAR系统成像原理的几何关系，计算目标点与机下点位置的相对偏离关系，再根据机下点的地理位置和载机的航向来获得目标的地理坐标；这类方法虽然能实现快速定位，但是其定位精度不高。而一些星载SAR自动目标定位方法，虽然考虑了地球模型对定位精度的影响，但是却忽略了迭代初值点对目标定位的影响；迭代初值点的设定影响着迭代算法的收敛方向，在初值点设置与SAR工作模式矛盾的情况下，会产生目标混淆的情况，使最终的定位结果产生很大的偏差。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提供一种以距离-多谱勒原理为基础的SAR图像自动目标定位方法，该方法是一种无控制点的绝对定位方法，既考虑了地球模型的影响，提高了定位精度；又通过泰勒级数展开将非线性方程组线性化，提高定位速度；同时根据SAR工作模式选择迭代初值点，避免了目标混淆。

本发明的技术解决方案为：一种SAR图像自动目标定位方法，其特征在于利用载机的飞行信息(位置、速度等)和SAR图像信息(分辨率、中心距离、中心多谱勒频率等)，在距离-多谱勒方程的基础上，结合地球模型方程，将目标在SAR图像中的坐标转换为目标在地球上的实际地理坐标，其具体步骤如下：

(1)确定载机位置参数，将载机位置的地理坐标(经度、纬度、高度)转换为空间直角坐标；

(2)确定载机飞行速度，将载机在雷达惯性坐标系中的速度矢量(北向、天向、东向)转换为空间直角坐标系中的速度矢量；

(3)根据SAR图像距离向和方位向的分辨率，以及图像原点的距离和多谱勒信息，求出图像坐标为(x，y)的目标点的距离R和多谱勒频率f_D；

(4)根据雷达波束照射方向，计算照射区域内距离为R的等距离圆上零多谱勒频率点的直角坐标，并以此坐标作为迭代计算的初值；

(5)联立距离-多谱勒方程和地球模型方程，利用泰勒级数展开将此非线性方程组线性化，进行迭代计算求得目标点的直角坐标；

(6)将目标点的直角坐标转化为目标在地球上的经纬度坐标，输出最终定位结果。

本发明的原理是：SAR传感器所发出之探测波的同一波前球面上的所有地物点将成像于一点，因此，在雷达照射区域内，分布着等时延的同心圆束和等多普勒频移的双曲线束。同一回波时延的点目标具有不同的多普勒频移，而具有相同多普勒频移的点目标又有不同的延时。根据距离向上回波信号的时延和方位向上的多普勒频移这两个信息，即已知目标点的距离-多谱勒方程，再结合地球模型方程，就可以将点目标区别确定。

本发明与现有技术相比的优点在于：本发明利用载机的飞行信息(位置、速度等)和SAR图像信息(分辨率、中心距离、中心多谱勒频率等)，将目标在SAR图像中的坐标转换为目标在地球上的实际地理坐标，在距离-多谱勒方程的基础上，结合地球模型方程进行目标定位，既可实现无控制点的绝对定位，又考虑了地球模型对定位精度的影响；较之于依赖高精度地面控制点的传统SAR目标定位方法，本发明以SAR成像的距离-多谱勒原理为基础，可以实现无控制点的绝对定位；较之于一般的以SAR成像几何关系为基础的机载SAR目标定位方法，本发明考虑了地球模型的影响，有着较高的定位精度；同时，本发明利用泰勒级数展开对非线性方程组线性化，通过迭代的方法求解方程，提高了算法的计算和定位速度，适合实时处理；最后，较之于其他SAR目标定位的迭代算法，本发明根据SAR的工作模式适当选择迭代计算的初值点，控制了迭代中收敛的方向，避免了目标混淆的发生。

总之，本发明作为一种无控制点的定位方法，具有精度高、速度快的优点，可进行实时处理，适合各种机载SAR系统。

附图说明

图1为本发明的SAR图像自动目标定位方法流程图；

图2为空间直角坐标系示意图；

图3为WGS84地理坐标系示意图；

图4为目标混淆示意图；

图5为初值计算平面图。

具体实施方式

首先简单对本发明中用到的几种不同的坐标系作简单的介绍：

(1)空间直角坐标系：以地球球心O为坐标原点，X轴是起始子午面与赤道面交线，指向格林威治子午线上半分支与赤道的交点；Z轴沿地球的自转轴指向正北极；Y轴也在赤道平面内，使该坐标系构成右手直角坐标系O-XYZ。在该坐标系中，P点的位置用(x，y，z)表示，如图2所示。

(2)WGS84地理坐标系：用一个参考椭球体模型化地球，椭球中心与地球质心重合，椭球的短轴与地球自转轴重合，椭球的起始子午面与格林尼治起始子午面重合。如图3所示，假设P为椭球面上的点，过P点的子午面与起始子午面夹角叫P点的经度L，由起始子午面起算，向东为东经，向西为西经；而过P点的椭球法线与赤道平面的夹角叫P点的纬度B，由赤道面起算，向北为北纬，向南为南纬。对于不在椭球面上的点，还需要另一参数——大地高H，定义为该点沿椭球法线方向至椭球面的距离。

(3)雷达惯性坐标系：以当地水平面为XOZ平面，X轴正方向指正北，Y轴过原点垂直于当地水平面，向上为正，而Z轴正方向则指向正东，XYZ三轴构成右手系。由于各地水平面不平行，所以在不同位置的雷达惯性坐标系之间也是存在一定夹角的。因此雷达惯性坐标系也称为局部坐标系，其坐标轴的指向受坐标原点所处的经纬度位置的影响。

(4)图像平面坐标系：SAR图像既无框标，也无像主点，故图像平面坐标系要人为选取。图像平面坐标通常以截取的影像中某一近距离点为坐标原点，以影像行方向(距离向)为X轴，平台运动方向(方位向)为Y轴，构成X-Y图像平面坐标系。

如图1所示，本发明方法基于距离-多谱勒原理和地球模型方程进行目标定位，将目标的图像坐标转换成地理坐标。方程组采用的是大地直角坐标系，因此首先需要通过步骤(1)和(2)将载机的位置和速度矢量归一化到直角坐标系中；而步骤(3)通过目标的图像坐标获得目标的距离和多谱勒信息，为列出距离-多谱勒方程提供条件；步骤(4)根据几何关系确定迭代初值后，利用步骤(5)的迭代算法求解出目标的大地直角坐标；最终由步骤(6)转换为经纬度坐标得到最后的结果，其具体实施步骤如下：

(1)确定载机位置参数，假设载机所在的地理坐标为(B_s，L_s，H_s)，B_s，L_s，H_s分别为载机所在位置的纬度、经度和高度，将此坐标转化为空间直角坐标(X_s，Y_s，Z_s)，如式(9)：

$\left[\begin{array}{c}{X}_s\\ Y_s\\ Z_s\end{array}\right]=\left[\begin{array}{c}(N+H)\cos B_s\cos L_s\\ (N+H)\cos B_s\sin L_s\\ (N+H)\sin B_s\end{array}\right]---\left(9\right)$

其中， $N=\sqrt{{\left(R_e\cos B_s\right)}^2+{\left(R_p\sin B_s\right)}^2}$ 为近似椭球卯酉圈曲率半径，R_e和R_p分别为地球的赤道半径和极半径的长度。

(2)确定载机飞行速度，所设载机在雷达惯性坐标系中的速度矢量为(V_N，V_U，V_E)，V_N，V_U，V_E分别为北向、天向、东向的速度矢量，将其转换为空间直角坐标系中的速度矢量(V_X，V_Y，V_Z)，如式(10)：

$\left(\begin{array}{c}{V}_X\\ V_Y\\ V_Z\end{array}\right)=\left(\begin{array}{ccc}-\sin B_s\cos L_s& \cos B_s\sin L_s& -\sin L_s\\ -\sin B_s\sin L_s& \cos B_s\sin L_s& \cos L_s\\ \cos B_s& \sin B_s& 0\end{array}\right)\left(\begin{array}{c}{V}_N\\ V_U\\ V_E\end{array}\right)---\left(10\right)$

其中，B_s和L_s分别为载机所在位置的纬度和经度。

(3)假设目标的图像坐标为(x，y)，根据SAR图像信息，求出目标点的距离R和多谱勒频率f_D，如式11：

$\left\{\begin{array}{c}R=R_0+x\·\mathrm{\Δx}\\ f_D=f_{D0}+y\·\mathrm{\Δy}(\frac{2v}{\mathrm{\λR}}-\frac{\mathrm{\λ}{f}_{D0}^2}{2\mathrm{vR}})\end{array}\right.---\left(11\right)$

其中，Δx和Δy分别为图像距离向和方位向的采样间隔，λ为雷达工作波长、v为载机飞行速度、R₀和f_D0分别为图像原点的距离和多谱勒频率。

(4)根据雷达波束照射方向，计算照射区域内距离为R的等距离圆上零多谱勒频率点的直角坐标，并以此坐标作为迭代计算的初值。

假设某时刻载机S的地理坐标为(B_s，L_s，H_s)，空间直角坐标为(X_s，Y_s，Z_s)；则其垂直地面的垂足G的地理坐标为(B_s，L_s，0)，空间直角坐标为(X_G，Y_G，Z_G)。此时，地球表面上与载机距离为R的等距线可看成一个以G为圆心的圆，如图4所示。在距离确定的情况下，目标的多谱勒频率可由目标的方位角唯一确定。而在图4中，假设载机的速度矢量v为(V_N，V_U，V_E)，T₁和T₂都在以点G为圆心的交线圆上，∠vST₁和∠vST₂大小都等于θ，此时T₁和T₂具有相同的多谱勒频率，从而产生了目标混淆。同时，从图4也可以看出，T₁和T₂分布在载机飞行航迹的两侧，因此可以根据雷达波束照射方向来选择适当的初值，使迭代过程收敛到与初值点同侧的目标点。当雷达照射方向向左时，选取点H₁为初值点；反之，当雷达照射方向向右时，则选取点H₂为初值点，H₁和H₂都是多谱勒频率为零的点，与飞行方向垂直。具体求解过程如下：

首先计算交线圆的半径。如图5所示，三角形SOH三条边的长度都是已知的：

$\left\{\begin{array}{c}a=\left|\right|\mathrm{OS}\left|\right|=\sqrt{X_s^2+Y_s^2+Z_s^2}\\ b=\left|\right|\mathrm{OH}\left|\right|=\left|\right|\mathrm{OG}\left|\right|=\sqrt{X_G^2+Y_G^2+Z_G^2}\\ c=\left|\right|\mathrm{SH}\left|\right|=R\end{array}\right.---\left(12\right)$

根据式(13)的余弦定理可以求出夹角θ的值：

$\mathrm{\θ}=\arccos \left[\frac{(a^2+b^2-c^2)}{2\mathrm{ab}}\right]---\left(13\right)$

因此可求得交线圆的半径：

r＝||GH||＝b·sinθ    (14)

在确定了交线圆的半径之后，需要确定初值点的方向，该方向与载机飞机方向垂直。载机速度v为(V_N，V_U，V_E)，若雷达波束照射方向向左，则初值点方向为(V_E，0，-V_N)；若雷达波束照射方向向右，则初值点方向为(-V_E，0，V_N)。将初值点的方向矢量转换为直角坐标(X_ver，Y_ver，Z_ver)，则初值点坐标(X_H，Y_H，Z_H)为：

$\left\{\begin{array}{c}{X}_H=X_G+X_{\mathrm{ver}}\·(r/v)\\ Y_H=Y_G+Y_{\mathrm{ver}}\·(r/v)\\ Z_H=Z_G+Z_{\mathrm{ver}}\·(r/v)\end{array}\right.---\left(15\right)$

(5)联立距离-多谱勒方程和地球模型方程，利用泰勒级数展开将此非线性方程组线性化，进行迭代计算求得目标点的直角坐标。

在SAR作用区域内，分布着等时延的同心圆束和等多普勒频移的双曲线束。同一回波时延的点目标具有不同的多普勒频移，而具有相同多普勒频移的点目标又有不同的延时。结合地球模型，根据距离向上回波信号的时延和方位向上的多普勒频移这两个信息，即可将点目标区别确定。

假设在直角坐标系中，载机位置的坐标为(X_s，Y_s，Z_s)，速度为(V_x，V_y，V_z)，则雷达与目标点的距离为：

(X_t-X_S)²+(Y_t-Y_S)²+(Z_t-Z_S)²＝R²    (16)

用矢量形式表示则有：

$R^2=({\stackrel{\→}{R}}_S-{\stackrel{\→}{R}}_t)\·({\stackrel{\→}{R}}_S-{\stackrel{\→}{R}}_t)---\left(17\right)$

对上式两边求导，有：

$R\stackrel{\·}{R}={\stackrel{\→}{V}}_S\·({\stackrel{\→}{R}}_S-{\stackrel{\→}{R}}_t)---\left(18\right)$

其中

为载机在空间直角坐标系中的速度矢量。考虑到回波的双程性，回波信号中由回波延迟带来的相位为：

$\mathrm{\φ}=-\frac{4\mathrm{\π}}{\mathrm{\λ}}R---\left(19\right)$

对应多普勒频率为：

$f_D=\frac{\stackrel{\·}{\mathrm{\φ}}}{2\mathrm{\π}}=-\frac{2}{\mathrm{\λ}}\stackrel{\·}{R}---\left(20\right)$

将式(18)代入式(20)，就得到了雷达成像的多普勒条件方程：

$f_D=-\frac{2}{\mathrm{\λR}}[{\stackrel{\→}{V}}_S\·({\stackrel{\→}{R}}_S-{\stackrel{\→}{R}}_t)]---\left(21\right)$

以坐标形式表示为：

$(X_t-X_S)V_x+(Y_t-Y_S)V_y+(Z_t-Z_S)V_z=-\frac{\mathrm{\λR}{f}_D}{2}---\left(22\right)$

同时，近似地球模型方程为

$\frac{X_t^2+Y_t^2}{{(R_e+h)}^2}+\frac{Z_t^2}{R_p^2}=1---\left(23\right)$

其中，R_e和R_p分别为地球的赤道半径和极半径的长度，h是相对假设模型的局部地域高度。

从以上分析可以看出，式(16)表示与载机距离为R的等距面，式(22)表示多谱勒频率为f_D的等多谱勒面，而式(23)表示目标为地球表面上的点，由于目标点就是等距面、等多谱勒面和地球模型三者的交点，因此在已知载机位置、速度和雷达工作参数的情况下，联立式(16)、式(22)和式(23)，求解如式(24)的方程组就能唯一地确定出目标点(X，Y，Z)的位置，从而实现SAR图像目标定位。

$\left\{\begin{array}{c}{F}_0={(X-X_S)}^2+{(Y-Y_S)}^2+{(Z-Z_S)}^2-R^2=0\\ F_1=(X-X_S)V_x+(Y-Y_S)V_y+(Z-Z_S)V_z+\mathrm{\λR}{f}_D/2=0\\ F_2=\frac{X^2+Y^2}{{(R_e+h)}^2}+\frac{Z^2}{R_p^2}-1=0\end{array}\right.---\left(24\right)$

注意到上式是非线性方程组，不能直接计算。为此，可推导出其一次项近似式，即进行线性化。一般线性化的方法是将非线性函数按泰勒级数展开，取至一次项，即得原函数的一次项近似公式。通常采用迭代算法求解，假定第k次迭代时的目标位置为(X_k，Y_k，Z_k)，式(24)线性化后的误差方程为：

ε＝AΔ+L    (25)

式中

$A=\left[\begin{array}{ccc}{A}_{11}& A_{12}& A_{13}\\ A_{21}& A_{22}& A_{23}\\ A_{31}& A_{32}& A_{33}\end{array}\right]---\left(26\right)$

$A_{11}=\frac{d{F}_0}{\mathrm{dX}}=2(X_k-X_S);$ $A_{12}=\frac{d{F}_0}{\mathrm{dY}}=2(Y_k-Y_S);$ $A_{13}=\frac{d{F}_0}{\mathrm{dZ}}=2(Z_k-Z_S);$

$A_{21}=\frac{d{F}_1}{\mathrm{dX}}=V_x;$ $A_{22}=\frac{d{F}_1}{\mathrm{dY}}=V_y;$ $A_{23}=\frac{d{F}_1}{\mathrm{dZ}}=V_z;$

$A_{31}=\frac{d{F}_2}{\mathrm{dX}}=\frac{2X_k}{{(R_e+h)}^2};$ $A_{32}=\frac{d{F}_2}{\mathrm{dY}}=\frac{2Y_k}{{(R_e+h)}^2};$ $A_{33}=\frac{d{F}_2}{\mathrm{dZ}}=\frac{2Z_k}{R_p^2};$

而Δ＝[dX dY dZ]^T为地面点坐标修正值，L＝[dF₀ dF₁ dF₂]^T为常数项矢量，其中：

dF₀＝(X_k-X_S)²+(Y_k-Y_S)²+(Z_k-Z_S)²-R²

dF₁＝(X_k-X_S)V_x+(Y_k-Y_S)V_y+(Z_k-Z_S)V_z+λRf_D/2    (27)

${\mathrm{dF}}_2=\frac{{X_k}^2+{Y_k}^2}{{(R_e+h)}^2}+\frac{{Z_k}^2}{R_p^2}-1$

从而可由最小二乘法求使得ε^Tε最小的坐标修正值

Δ＝-(A^TA)^-1A^TL    (28)

将修正值加上X，Y，Z得：X_k+1＝X_k+dX，Y_k+1＝Y_k+dY，Z_k+1＝Z_k+dZ，这样循环迭代，直到方程解的误差小于给定限差为止，最后得到目标点的大地直角坐标。

(6)将目标点的直角坐标(X_t，Y_t，Z_t)转化为经纬度坐标(B_t，L_t)，输出最终定位结果，如式(29)：

$\left[\begin{array}{c}{L}_t\\ B_t\end{array}\right]=\left[\begin{array}{c}\arctan \left(\frac{Y_t}{X_t}\right)\\ \arctan \left(\frac{Z_t}{\sqrt{X_t^2+Y_t^2}}\right)\end{array}\right]---\left(29\right)$

本发明说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。

尽管为说明目的公开了本发明的最佳实施例和附图，但是本领域的技术人员可以理解：在不脱离本发明及所附的权利要求的精神和范围内，各种替换、变化和修改都是可能的。因此，本发明不应局限于最佳实施例和附图所公开的内容。

Claims (4)

1、一种机载SAR图像自动目标定位方法，其特征在于包括步骤如下：

(1)确定载机位置参数，将载机位置的地理坐标，包括经度、纬度、高度转换为空间直角坐标；

(2)确定载机飞行速度，将载机在雷达惯性坐标系中的速度矢量，包括北向、天向、东向转换为空间直角坐标系中的速度矢量；

(3)根据SAR图像距离向和方位向的分辨率，以及图像原点的距离和多谱勒信息，求出图像坐标为(x，y)的目标点距离R和多谱勒频率f_D；

(4)根据雷达波束照射方向，计算照射区域内目标点距离为R的等距离圆上零多谱勒频率点的直角坐标，并以此坐标作为迭代计算的初值；

(5)联立距离-多谱勒方程和地球模型方程，利用泰勒级数展开将此非线性方程组线性化，进行迭代计算求得目标点的直角坐标；

(6)将目标点的直角坐标转化为目标在地球上的经纬度坐标，输出最终定位结果。

2、根据权利要求1所述的SAR图像自动目标定位方法，其特征在于：所述的步骤(3)中，假设图像距离向和方位向的采样间隔Δx和Δy，利用雷达工作波长λ、载机飞行速度v、图像原点的距离R₀和多谱勒频率f_D0，则图像坐标为(x，y)的目标点的距离R和多谱勒频率f_D如下：

$\left\{\begin{array}{cc}R& =R_0+x\·\mathrm{\Δx}\\ f_D& =f_{D0}+y\·\mathrm{\Δy}(\frac{2v}{\mathrm{\λR}}-\frac{\mathrm{\λ}{f}_{D0}^2}{2\mathrm{vR}})\end{array}\right.---\left(1\right).$

3、根据权利要求1所述的SAR图像自动目标定位方法，其特征在于：所述的步骤(4)中，若载机和机下点的直角坐标分别为(x_s，y_s，z_s)和(x_G，y_G，z_G)，假设 $a=\sqrt{x_s^2+y_s^2+z_s^2},$ $b=\sqrt{x_G^2+y_G^2+z_G^2},$ c＝R，则以机下点为圆心的距离为R的等距离圆的半径r为：

$r=b\·\sin \left\{\arccos \right[\frac{(a^2+b^2-c^2)}{2\mathrm{ab}}\left]\right\}---\left(2\right)$

假设在雷达惯性坐标系中，载机的速度v＝(N，U，E)，若雷达波束照射方向向左，则初值点方向为(E，0，-N)；若雷达波束照射方向向右，则初值点方向为(-E，0，N)，将初值点方向转换为直角坐标(x_ver，y_ver，z_ver)，则初值点坐标为：

$\left\{\begin{array}{c}{x}_H=x_G+x_{\mathrm{ver}}\·(r/v)\\ y_H=y_G+y_{\mathrm{ver}}\·(r/v)\\ z_H=z_G+z_{\mathrm{ver}}\·(r/v)\end{array}\right.---\left(3\right).$

4、根据权利要求1所述的SAR图像自动目标定位方法，其特征在于：所述的步骤(5)中，假设在直角坐标系中，载机位置的坐标为(X_s，Y_s，Z_s)，速度为(V_x，V_y，V_z)，R_e和R_p为地球的赤道半径和极半径，h为目标局部高度，则目标点坐标(X，Y，Z)满足以下的由距离-多谱勒方程和地球模型方程联立的方程组：

$\left\{\begin{array}{c}{F}_0={(X-X_S)}^2+{(Y-Y_S)}^2+{(Z-Z_S)}^2-R^2=0\\ F_1=(X-X_S)V_x+(Y-Y_S)V_y+(Z-Z_S)V_z+\mathrm{\λR}{f}_D/2=0\\ F_2=\frac{X^2+Y^2}{{(R_e+h)}^2}+\frac{Z^2}{R_p^2}-1=0\end{array}\right.---\left(4\right)$

利用泰勒级数展开将方程组线性化，线性化后的误差方程为

ε＝AΔ+L    (5)

式中

$A=\left[\begin{array}{ccc}{A}_{11}& A_{12}& A_{13}\\ A_{21}& A_{22}& A_{23}\\ A_{31}& A_{32}& A_{33}\end{array}\right]---\left(6\right)$

$A_{11}=\frac{d{F}_0}{\mathrm{dX}}=2(X_k-X_S);$ $A_{12}=\frac{d{F}_0}{\mathrm{dY}}=2(Y_k-Y_S);$ $A_{13}=\frac{d{F}_0}{\mathrm{dZ}}=2(Z_k-Z_S);$

$A_{21}=\frac{d{F}_1}{\mathrm{dX}}=V_x;$ $A_{22}=\frac{d{F}_1}{\mathrm{dY}}=V_y;$ $A_{23}=\frac{d{F}_1}{\mathrm{dZ}}=V_z;$

$A_{31}=\frac{d{F}_2}{\mathrm{dX}}=\frac{2X_k}{{(R_e+h)}^2};$ $A_{32}=\frac{d{F}_2}{\mathrm{dY}}=\frac{2Y_k}{{(R_e+h)}^2};$ $A_{33}=\frac{d{F}_2}{\mathrm{dZ}}=\frac{2Z_k}{R_p^2};$

而Δ＝[dX dY dZ]^T为地面点坐标修正值，L＝[dF₀ dF₁ dF₂]^T为常数项矢量，其中：

dF₀＝(X_k-X_S)²+(Y_k-Y_S)²+(Z_k-Z_S)²-R²

dF₁＝(X_k-X_S)V_x+(Y_k-Y_S)V_y+(Z_k-Z_S)V_z+λRf_D/2    (7)

$d{F}_2=\frac{{X_k}^2+{Y_k}^2}{{(R_e+h)}^2}+\frac{{Z_k}^2}{R_p^2}-1$

从而可由最小二乘法求使得ε^Tε最小的坐标修正值

Δ＝-(A^TA)^-1A^TL    (8)

将修正值加上X，Y，Z得：X_k+1＝X_k+dX，Y_k+1＝Y_k+dY，Z_k+1＝Z_k+dZ，这样循环迭代，直到方程解的误差小于给定限差为止，最后得到地面点坐标。

Images