CN102607533A - 线阵ccd光学与sar影像联合区域网平差定位方法 - Google Patents

Abstract

本发明涉及线阵CCD光学与SAR影像联合区域网平差定位方法，本发明是针对线阵CCD光学与SAR遥感影像联合区域网平差方法，实现方案是首先进行区域网概算，设定每张影像的外定向参数和加密点坐标的初值，然后根据联合区域网平差模型建立误差方程式和相应的权矩阵，接着按照最小二乘原理对误差方程进行法化答解，得到每张影像的外定向参数，最后利用多片前方交会计算加密点的地面坐标。本发明遥感影像联合区域网平差定位方法的目的是将不同传感器或平台获取的同一区域的遥感影像组合起来进行平差计算，以获得所有影像的外定向参数与加密点的地面坐标。

Description

线阵CCD光学与SAR影像联合区域网平差定位方法

技术领域

本发明属于摄影测量与遥感技术，图像处理领域，特别是涉及摄影测量领域中的一种针对线阵CCD光学与SAR遥感影像联合区域网平差定位方法。

背景技术

在立体摄影测量技术中完成影像定向需要在每幅影像上选择尽可能多且分布均匀的控制点，这是一项非常繁琐的工作。摄影测量中的区域网平差是解决这一问题的重要手段，其通过对同一区域的多幅影像进行联合处理，利用少量的控制点即可实现所有影像的定向并同时获得更加密集的控制点。区域网平差是摄影测量中的一项重要工作，是生成大范围DEM的必要工序。

在“复合式”像对立体定位方法中，SAR影像的定向精度往往不高，其主要原因就是控制点的量测很困难且精度不高，而采用不同类型遥感影像进行联合区域网平差的方法，则可以尽可能地减少SAR影像上控制点的数量，而利用精度较高的光学影像上的控制点来计算SAR影像的定向参数，同时获得较高精度的地面点坐标。

早在模拟摄影测量阶段，人们已经在全能型立体测图仪上对航空像片进行空中三角测量。在二十世纪60年代，由于计算机的诞生而发展了解析空中三角测量。二十世纪60-70年代，解析空中三角测量被应用于雷达航空影像，成功实现了对云层遮盖地区的测图，但精度较低。传统区域网平差的研究对象主要是同一种传感器获取的影像，而对于多源遥感影像联合区域网平差的研究则很少得到关注。Young-Ran Lee(Young-Ran Lee，2003)采用改进的共线条件方程，对Landsat-7、KOMPSAT-1、SPOT-1、IKONOS等光学卫星影像进行了联合区域网平差试验，在5个控制点、10个连接点和10个检查点的情况下，获得平面2.76m，高程3.1m的精度。Thierry Toutin(ThierryToutin，2004)尝试了将多种光学和SAR卫星遥感影像联立进行光束法区域网平差，其使用了Landsat-7ETM+、SPOT-4HRV、ASTER、RADARSAT-1、ERS-1五颗卫星共49幅影像进行了实验，取得了25-30m的定位精度。他在2006年(Thierry Toutin，2006)又对SPOT-5HRS、SPOT-5HRG、IKONOS、Quick Bird等影像进行了多星多传感器遥感影像的联合区域网平差试验，获得平面2.54m、高程2.1m的精度。但上述研究主要是实验验证了多源遥感影像联合区域网平差思想的可行性，并未提出相关的方法、模型。国内目前尚未见相关技术的研究成果。

发明内容

本发明遥感影像联合区域网平差定位方法的目的是将不同传感器或平台获取的同一区域的遥感影像组合起来进行平差计算，以获得所有影像的外定向参数与加密点的地面坐标。

线阵CCD光学与SAR影像联合区域网平差定位方法，步骤如下：1)对于同一区域两幅及两幅以上线阵CCD光学和/或SAR遥感图像，任意选择两幅或两幅以上像片；

2)每幅像片使用其各自的控制点，根据每幅像片对应的构像模型，通过空间后方交会确定每幅像片的外方位元素初值，然后由选择的所有像片通过前方交会确定所有加密点地面坐标的初值；

3)将步骤2)得到的各加密点地面坐标、各像片的外方位元素的初值带入各像片分别对应的构像模型，得到加密点在各像片上对应像点的坐标；

4)根据加密点和控制点所在像片类型的构像模型，建立加密点和控制点的误差方程式并进行解算，得到各像片的外方位元素和加密点的地面坐标的近似值；

5)重复3)——4)，直到连续两次得到的各像片的外方位元素改正数和加密点地面坐标改正数小于设定的限差；

6)采用多片前方交会，计算未参加区域网平差各点的地面坐标，实现定位。

所述线阵CCD光学影像采用线中心投影构像模型；所述SAR影像采用F.Leberl构像模型。步骤4)包括根据所述两构像模型的联合区域网平差模型建立误差方程式和相应的权矩阵，接着按照最小二乘原理对误差方程进行法化答解。步骤4)中，如果是线阵CCD光学影像，增加伪观测值误差方程式。

本发明是针对线阵CCD光学与SAR遥感影像联合区域网平差方法，实现方案是首先进行区域网概算，设定每张影像的外定向参数和加密点坐标的初值，然后根据联合区域网平差模型建立误差方程式和相应的权矩阵，接着按照最小二乘原理对误差方程进行法化答解，得到每张影像的外定向参数，最后利用多片前方交会计算加密点的地面坐标。

SAR影像的定向精度往往不高，其主要原因就是控制点的量测很困难且精度不高。而采用不同类型遥感影像进行联合区域网平差的方法，则可以尽可能地减少SAR影像上控制点的数量，而利用精度较高的光学影像上的控制点来计算SAR影像的定向参数，同时获得较高精度的地面点坐标。本发明遥感影像联合区域网平差定位方法能够有效减少摄影测量作业中所需地面控制点的数量，并能够同时将影像的外定向参数和加密点坐标一起解算出来，可以有效提高影像间的相对精度，并得到更多分布均匀的地面控制点。相比传统的区域网平差方法具有更广泛的适用性；以线中心投影模型和F.Leberl模型为基础构建了联合区域网平差模型，相比采用统一近似构像模型的方法具有更高的精度和可靠性；能够有效减少摄影测量作业中所需地面控制点的数量，并能够同时将影像的外定向参数和加密点坐标一起解算出来；能够有效提高影像间的相对定向精度，并得到更多分布均匀的地面控制点。

附图说明

图1线阵CCD的物像关系；

图2SAR传感器距离条件示意图；

图3SAR传感器多普勒条件示意图；

图4线阵CCD光学与SAR遥感影像联合区域网平差方法的流程图；

图5四幅卫星影像对应区域及控制点分布示意图；

图6-图16为(1)到(11)所有方案的影像覆盖范围以及控制点、加密点的分布示意图。(图5、图6中‘△’表示SPOT-5影像上的控制点，表示SPOT-4影像上的控制点，‘*’表示Radarsat-1影像上的控制点，‘+’表示ERS-2影像上的控制点，‘o’表示加密点)

具体实施方式

本发明是一种基于多种传感器构像模型的“复合式”像对立体定位方法，其核心内容是联合区域网平差模型，它是在线阵CCD光学影像和SAR影像平差模型的基础上建立起来的。

1、两种平差模型

(1)线阵CCD光学影像的平差模型

线阵CCD影像通常采用线中心投影模型，是因为它是由线阵列传感器沿飞行方向推扫而成的，每一扫描行图像与被摄物体之间具有严格的中心投影关系，每一扫描行都具有各自的外方位元素。但在平台飞行过程中，其姿态变化可认为是相当平稳的，所以可以假设每景影像的像平面坐标原点为中央扫描行的中点，每一扫描行的外方位元素随y值变化，则外方位元素可描述为：

式中为第i扫描行的外方位元素，y为该扫描行沿飞行方向的像平面坐标，

为中心扫描行的外方位元素，

为外方位元素的一阶变率。

如图1所示，i扫描行上的像点与相应地面点间的中心投影关系为：

$\left.\begin{array}{c}{x}_i=-f\frac{a_1(X-X_{\mathrm{si}})+b_1(Y-Y_{\mathrm{si}})+c_1(Z-Z_{\mathrm{si}})}{a_3(X-X_{\mathrm{si}})+b_3(Y-Y_{\mathrm{si}})+c_3(Z-Z_{\mathrm{si}})}\\ 0=-f\frac{a_2(X-X_{\mathrm{si}})+b_2(Y-Y_{\mathrm{si}})+c_2(Z-Z_{\mathrm{si}})}{a_3(X-X_{\mathrm{si}})+b_3(Y-Y_{\mathrm{si}})+c_3(Z-Z_{\mathrm{si}})}\end{array}\right\}---\left(2\right)$

写成矩阵形式为：

$\left[\begin{array}{c}{x}_i\\ 0\\ -f\end{array}\right]=\frac{1}{\mathrm{\λ}}{M}_i^T\left[\begin{array}{c}X-X_{\mathrm{Si}}\\ Y-Y_{\mathrm{Si}}\\ Z-Z_{\mathrm{Si}}\end{array}\right]=\frac{1}{\mathrm{\λ}}{\left[\begin{array}{ccc}{a}_1& a_2& a_3\\ b_1& b_2& b_3\\ c_3& c_2& c_3\end{array}\right]}^T\left[\begin{array}{c}X-X_{\mathrm{Si}}\\ Y-Y_{\mathrm{Si}}\\ Z-Z_{\mathrm{Si}}\end{array}\right]---\left(3\right)$

式中(x_i，0)为第i行上像点的像平面坐标，(X，Y，Z)为其对应地面点的地面坐标，(X_si，Y_si，Z_si)为第i行的摄站坐标，λ为比例尺系数，a_i，b_i，c_i(i＝1，2，3)是由第i行的外方位角元素

ω_i，κ_i所确定的旋转矩阵中的九个元素。

对(2)式线性化后可得：

式中为各外方位元素的改正数，dX、dY、dZ为地面点坐标改正数，C_ij(i＝1，2；j＝1，…，6)为各改正数的系数，l_x、l_y为常数项。

$\stackrel{\‾}{Z}=a_3(X-X_S)+b_3(Y-Y_S)+c_3(Z-Z_S)---\left(6\right)$

对于每张像片上的每个加密点，都可按照(4)式列出一组误差方程。而对于控制点，则需要根据控制点的类型采用不同形式的误差方程。如果是平高控制点，则上式中的dX＝dY＝dZ＝0；如果是平面控制点，则上式中的dX＝dY＝0；如果是高程控制点，则上式中的dZ＝0。

如果考虑控制点存在一定误差，则可以将控制点作为带权观测值，假设此时地面点坐标近似值取为该点已知值，则有

$\left[\begin{array}{c}{V}_X\\ V_Y\\ V_Z\end{array}\right]=\left[\begin{array}{c}\mathrm{dX}\\ \mathrm{dY}\\ \mathrm{dZ}\end{array}\right]-\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right]$ 权P_T(7)

由于线阵CCD光学影像的外方位元素间存在较强的相关性，为了保证法方程的正确答解，需要进行一些特殊的处理。

一种方法是增加一组“伪观测值”的误差方程式。引进如下的“伪观测值”

ω＝0；κ＝0；X_S0＝近似值；Y_S0＝0；Z_S0＝h

${\stackrel{\·}{Z}}_S=0$

把它们作为带权观测值，且取近似值为伪观测值，有如下误差方程式：

权PS(8)

基本误差方程式中包含两类未知数，即每景影像的外方位元素和加密点的地面坐标，其中后一类未知参数的数目往往比较巨大，并且随连接点的多寡而变化。此时可以对误差方程式进行改化，建立等效误差方程，先解算各景影像的外方位元素，再利用多片前方交会的方法，逐点计算加密点的地面坐标。在计算机满足要求的情况下，也可以同时答解出影像的外方位元素和加密点的地面坐标。

(2)SAR影像的平差模型

F.Leberl模型是基于SAR成像原理来表达雷达图像瞬间构像的数学模型。

①距离条件

如图2所示，D_S0为近距边扫描延迟，V_S为沿飞行方向的速度，R_S为天线中心S到地面点P的斜距，H为天线中心S到数据归化平面(基准面)的航高，y_S为地面点P在斜距显示图像上的距离向像坐标，y_g为地面点P在地距显示图像上的距离向像坐标，R₀为扫描延迟在数据归化平面上的投影，M_y为斜距显示图像的距离向像元分辨率，m_y为平距显示图像的距离向像元分辨率。

对于斜距显示图像有：

(X-X_S)²+(Y-Y_S)²+(Z-Z_S)²＝(y_sM_y+D_S0)²(9)

式中：(X，Y，Z)为地面点P的物方空间坐标，(X_S，Y_S，Z_S)为天线中心瞬时位置S的物方空间坐标(是飞行时间T的多项式函数)，可表示为：

$\left\{\begin{array}{c}{X}_S=X_{S0}+X_{V0}T+\frac{1}{2}{X}_{a0}{T}^2+......\\ Y_S=Y_{S0}+Y_{V0}T+\frac{1}{2}{Y}_{a0}+T^2+......\\ Z_S=Z_{S0}+Z_{V0}T+\frac{1}{2}{Z}_{a0}{T}^2+......\\ T=K_x\·x\end{array}\right.---\left(10\right)$

其中X_S0，Y_S0，Z_S0为对应于像坐标原点的雷达天线中心瞬时物方空间坐标，X_V0，Y_V0，Z_V0为飞行器对应于像坐标原点的速度矢量的分量(外方位元素的一阶变率)，X_a0，T_a0，Z_a0为飞行器对应于像坐标原点的加速度矢量的分量(外方位元素的二阶变率)，T为像坐标x相对于原点时刻的飞行时间，x为雷达图像的方位向像平面坐标，K_x为图像方位向每行所用的扫描时间。

设任意时刻天线中心的速度矢量为(X_V，Y_V，Z_V)，由(10)式可得：

$\left\{\begin{array}{c}{X}_V=\frac{{\∂X}_S}{\∂T}=X_{V0}+X_{a0}T+......\\ Y_V=\frac{{\∂Y}_S}{\∂T}=Y_{V0}+Y_{a0}T+......\\ Z_V=\frac{\∂Z_S}{\∂T}=Z_{V0}+Z_{a0}T+......\end{array}\right.---\left(11\right)$

同理，对于平距显示图像有：

(X-X_S)²+(Y-Y_S)²+(Z-Z_S)²＝(y_gm_y+R₀)²+H²(12)

②多普勒条件

如图3所示，τ为等多普勒圆锥面与垂直方向的夹角，D-XYZ为地面坐标系，其它变量定义同图2。

雷达图像的多普勒条件公式为

$f_{\mathrm{DC}}=\frac{2({\stackrel{\→}{V}}_S-{\stackrel{\→}{V}}_P)\·(\stackrel{\→}{S}-\stackrel{\→}{P})}{\mathrm{\λ}|\stackrel{\→}{S}-\stackrel{\→}{P}|}=\frac{2|{\stackrel{\→}{V}}_S-{\stackrel{\→}{V}}_P|}{\mathrm{\λ}}\·\frac{({\stackrel{\→}{V}}_S-{\stackrel{\→}{V}}_P)\·(\stackrel{\→}{S}-\stackrel{\→}{P})}{|{\stackrel{\→}{V}}_S-{\stackrel{\→}{V}}_P|\·|\stackrel{\→}{S}-\stackrel{\→}{P}|}=\frac{2|{\stackrel{\→}{V}}_S-{\stackrel{\→}{V}}_P|}{\mathrm{\λ}}\·\sin \mathrm{\τ}$

(13)

其中f_DC为多普勒中心频率，λ为雷达波长，

和

分别为S和P的速度矢量，和

分别为S和P的位置矢量。假设地面点是静止的，即

此时(13)式可以简化为

$X_V(X-X_S)+Y_V(Y-Y_S)+Z_V(Z-Z_S)=-\frac{\mathrm{\λ}{R}_S}{2}{f}_{\mathrm{DC}}---\left(14\right)$

其中R_S为地面点P的斜距。

F.Leberl采用的是正侧视角条件或零多普勒条件，即τ＝0，此时卫星飞行速度矢量与S指向P的位置矢量垂直，f_DC＝0，故称为零多普勒条件。此时(14)式变为

X_V(X-X_S)+Y_V(Y-Y_S)+Z_V(Z-Z_S)＝0(15)

(9)或(12)和(15)式共同构成了F.Leberl模型，构像方程中有X_S0、Y_S0、Z_S0、X_V0、Y_V0、Z_V0、X_a0、Y_a0、Z_a0、D_S0、K_x、M_y(或m_y)12个定向参数，通常D_S0、K_x、M_y(或m_y)能够由雷达的系统参数直接给定。

若用F₁和F₂分别表示(9)和(15)式，则有：

$\left\{\begin{array}{c}{F}_1={(y{M}_y+D_{S0})}^2-{(X-X_S)}^2-{(Y-Y_S)}^2-{(Z-Z_S)}^2=0\\ F_2=X_V(X-X_S)+Y_V(Y-Y_S)+Z_V(Z-Z_S)=0\end{array}\right.---\left(16\right)$

式(16)为非线性函数，不能直接用于平差计算，因此，需要对其进行线性化。其线性化后可得：

$\left\{\begin{array}{c}{A}_{11}d{X}_{S0}+A_{12}d{Y}_{S0}+A_{13}d{Z}_{S0}+A_{14}{\mathrm{dX}}_{V0}\\ +A_{15}d{Y}_{V0}+A_{16}{\mathrm{dZ}}_{V0}+A_{17}{\mathrm{dX}}_{a0}+A_{18}{\mathrm{dY}}_{a0}+A_{19}{\mathrm{dZ}}_{a0}\\ +B_{11}\mathrm{dX}+B_{12}\mathrm{dY}+B_{13}\mathrm{dZ}-F_1=C_{11}{v}_x+C_{12}{v}_y\\ A_{21}d{X}_{S0}+A_{22}d{Y}_{S0}+A_{23}d{Z}_{S0}+A_{24}d{X}_{V0}\\ +A_{25}d{Y}_{V0}+A_{26}d{Z}_{V0}+A_{27}{\mathrm{dX}}_{a0}+A_{28}{\mathrm{dY}}_{a0}+A_{29}d{Z}_{a0}\\ +B_{21}\mathrm{dX}+B_{22}\mathrm{dY}+B_{23}\mathrm{dZ}-F_2=C_{21}{v}_x+C_{22}{v}_y\end{array}\right.---\left(17\right)$

其中：

定向参数的改正数为

Δ＝[dX_S0，dY_S0，dZ_S0，dX_V0，dY_V0，dZ_V0，dX_a0，dY_a0，dZ_a0]^T

常数项为L＝-[F₁，F₂]^T

系数阵 $A=\left[\begin{array}{ccccccccc}{A}_{11}& A_{12}& A_{13}& A_{14}& A_{15}& A_{16}& A_{17}& A_{18}& A_{19}\\ A_{21}& A_{22}& A_{23}& A_{24}& A_{25}& A_{26}& A_{27}& A_{28}& A_{29}\end{array}\right]$

$B=\left[\begin{array}{ccc}{B}_{11}& B_{12}& B_{13}\\ B_{21}& B_{22}& B_{23}\end{array}\right]$

$C=\left[\begin{array}{cc}{C}_{11}& C_{12}\\ C_{21}& C_{22}\end{array}\right]$

且有：

$\left\{\begin{array}{c}{A}_{11}=\∂F_1/\∂X_{S0}=2(X-X_S)\\ A_{12}=\∂F_1/\∂Y_{S0}=2(Y-Y_S)\\ A_{13}=\∂F_1/\∂Y_{S0}=2(Z-Z_S)\\ A_{14}=\∂F_1/\∂X_{V0}=2(X-X_S)T\\ A_{15}=\∂F_1/\∂Y_{V0}=2(Y-Y_S)T\\ A_{16}=\∂F_1/\∂Z_{V0}=2(Z-Z_S)T\\ A_{17}=\∂F_1/\∂X_{a0}=(X-X_S)T^2\\ A_{18}=\∂F_1/\∂Y_{a0}=(Y-Y_S)T^2\\ A_{19}=\∂F_1/\∂Z_{a0}=(Z-Z_S)T^2\end{array}\right.$ $\left\{\begin{array}{c}{A}_{21}=\∂F_2/\∂X_{S0}=-X_V\\ A_{22}=\∂F_2/\∂Y_{S0}=-Y_V\\ A_{23}=\∂F_2/\∂Z_{S0}=-Z_V\\ A_{24}=\∂F_2/{\∂X}_{V0}=(X-X_S)-X_{V}T\\ A_{25}=\∂F_2/\∂Y_{V0}=(Y-Y_S)-Y_{V}T\\ A_{26}=\∂F_2/\∂Z_{V0}=(Z-Z_S)-Z_{V}T\\ A_{27}={\∂F}_2/\∂X_{a0}=(X-X_S)T-X_V{T}^2/2\\ A_{28}=\∂F_2/\∂Y_{a0}=(Y-Y_S)T-X_V{T}^2/2\\ A_{29}=\∂F_2/\∂Z_{a0}=(Z-Z_S)T-X_V{T}^2/2\end{array}\right.$

$\left\{\begin{array}{c}{B}_{11}=\∂F_1/\∂X=-2(X-X_S)\\ B_{12}=\∂F_1/\∂Y=-2(Y-Y_S)\\ B_{13}=\∂F_1/\∂Z=-2(Z-Z_S)\end{array}\right.$ $\left\{\begin{array}{c}{B}_{21}=\∂F_2/\∂X=X_v\\ B_{22}=\∂F_2/\∂Y=Y_V\\ B_{23}=\∂F_2/\∂Z=Z_V\end{array}\right.$

$\left\{\begin{array}{c}{C}_{11}=\∂F_1/\∂x=2K_x\·F_2\\ C_{12}=\∂F_1/\∂y=2M_y(y{M}_y+D_{S0})\end{array}\right.$ $\left\{\begin{array}{c}{C}_{21}=\∂F_2/\∂X=-K_x({X_V}^2+{Y_V}^2+{Z_V}^2)\\ C_{22}=\∂F_2/\∂Y=0\end{array}\right.$

其中dX_S0...dZ_a0为各外方位元素的改正数，dX、dY、dZ为地面点坐标改正数，A_ij为各外方位元素改正数的系数，B_ij为各地面点坐标改正数的系数，F₁、F₂为常数项，v_x、v_y为像坐标残差，C_ij为像坐标残差的系数。

①像点观测值的误差方程式

对于每张像片上的每个加密点和平高控制点，可按照(17)式列出一组误差方程式；

②控制点的误差方程式

当控制点作为带权观测值时，可按照(7)式列出一组误差方程式。

误差方程式中同样包含有两类未知数：一类是每张像片的定向参数，另一类是每个加密点的地面坐标。

答解未知数有多种途径：

整体答解：联立所有误差方程式进行法化，然后直接答解法方程组，一并解出像片的定向参数和加密点的地面坐标。

建立等效误差方程式：消去加密点的地面坐标，得到简化法方程式，平差答解各帧像片的定向参数，然后再用多片前方交会逐点计算加密点的地面坐标。

两类未知数交替趋近法：先利用单片后方交会计算各帧像片的定向参数，再通过多片前方交会计算加密点的地面坐标，两个步骤交替进行。

2、联合区域网平差模型

(1)误差方程式

根据(4)式，对于线阵CCD光学影像上的每个加密点或控制点，可列出一组误差方程式，即

其中

为各外方位元素的改正数，dX dY dZ为地面点坐标改正数

为各外方位元素改正数的系数，

为各地面点坐标改正数的系数，

为常数项，

为像坐标改正数。(5-3)式写成矩阵形式为

A^L·Δ^L+B^L·X^L-L^L＝V^L    权阵P^L    (19)

根据(17)式，对于SAR影像上的每个加密点或控制点，也可列出一组误差方程式，即

$\left.\begin{array}{c}{A}_{11}^{S}d{X}_{s0}+A_{12}^{S}d{Y}_{s0}+A_{13}^{S}d{Z}_{s0}+A_{14}^S{\mathrm{dX}}_{v0}+A_{15}^S{\mathrm{dY}}_{v0}+A_{16}^{S}d{Z}_{v0}\\ +A_{17}^{S}d{X}_{a0}+A_{18}^{S}d{Y}_{a0}+A_{19}^S{\mathrm{dZ}}_{a0}\\ +B_{11}^S\mathrm{dX}+B_{12}^S\mathrm{dY}+B_{13}^S\mathrm{dZ}-l_x^S=C_{11}{v}_x^S+C_{12}{v}_y^S\\ A_{21}^S{\mathrm{dX}}_{s0}+A_{22}^S{\mathrm{dY}}_{s0}+A_{23}^S{\mathrm{dZ}}_{s0}+A_{24}^S{\mathrm{dX}}_{v0}+A_{25}^{S}d{Y}_{v0}+A_{26}^S{\mathrm{dZ}}_{v0}\\ +A_{27}^S{\mathrm{dX}}_{a0}+A_{28}^S{\mathrm{dY}}_{a0}+A_{29}^S{\mathrm{dZ}}_{a0}\\ +B_{21}^S\mathrm{dX}+B_{22}^S\mathrm{dY}+B_{23}^S\mathrm{dZ}-l_y^S=C_{21}{v}_x^S+C_{22}{v}_y^S\end{array}\right\}---\left(20\right)$

其中dX_s0...dZ_a0为各外方位元素的改正数，C_ij为像坐标改正数的系数，其它各项与(18)式意义相同。(20)式写成矩阵形式为

A^S·Δ^S+B^S·X^S-L^S＝C·Vs^S＝V^S    权阵P^S    (21)

根据(19)和(21)式，对重叠范围内的所有加密点可列出下面一组误差方程

$\left.\begin{array}{ccc}{A}_P^L\·{\mathrm{\Δ}}^L& +B_P^L\·X_P-L_P^L=V_P^L& P_P^L\\ A_P^S\·{\mathrm{\Δ}}^S& +B_P^S\·X_P-L_P^S=V_P^S& P_P^S\end{array}\right\}---\left(22\right)$

写成矩阵形式为

$\left(\begin{array}{ccc}{A}_P^L& 0& B_P^L\\ 0& A_P^S& B_P^S\end{array}\right)\left(\begin{array}{c}{\mathrm{\Δ}}^L\\ {\mathrm{\Δ}}^S\\ X_P\end{array}\right)-\left(\begin{array}{c}{L}_P^L\\ L_P^S\end{array}\right)=\left(\begin{array}{c}{V}_P^L\\ V_P^S\end{array}\right)$ 权阵 $\left(\begin{array}{cc}{P}_P^L& 0\\ 0& P_P^S\end{array}\right)---\left(23\right)$

对所有平高控制点也可得到形如(23)式的误差方程

$\left(\begin{array}{ccc}{A}_G^L& 0& B_G^L\\ 0& A_G^S& B_G^S\end{array}\right)\left(\begin{array}{c}{\mathrm{\Δ}}^L\\ {\mathrm{\Δ}}^S\\ X_G\end{array}\right)-\left(\begin{array}{c}{L}_G^L\\ L_G^S\end{array}\right)=\left(\begin{array}{c}{V}_G^L\\ V_G^S\end{array}\right)$ 权阵 $\left(\begin{array}{cc}{P}_G^L& 0\\ 0& P_G^S\end{array}\right)---\left(24\right)$

上式中，如果认为控制点坐标准确，则X_G为零；若考虑到控制点有误差时，则X_G为非零矩阵，且要增加一组控制点观测误差方程

$E\·X_G-L_G^X=V_G^X$ 权阵 $P_G^X$ (25)

此外，对于线阵CCD光学影像，为保证在定向参数高度相关的情况下解的稳定性，还需要引进如下的“伪观测值”误差方程

$E\·{\mathrm{\Δ}}^L-L_{\mathrm{\Δ}}^L=V_{\mathrm{\Δ}}^L$ 权阵 $P_{\mathrm{\Δ}}^L$ (26)

将(23)-(26)式联立可构成如下的误差方程组

$\left(\begin{array}{cccc}{A}_P^L& & B_P^L& \\ & A_P^S& B_P^S& \\ A_G^L& & & B_G^L\\ & A_G^S& & B_G^S\\ & & & E\\ E& & & \end{array}\right)\left(\begin{array}{c}{\mathrm{\Δ}}^L\\ {\mathrm{\Δ}}^S\\ X_P\\ X_G\end{array}\right)-\left(\begin{array}{c}{L}_P^L\\ L_P^S\\ L_G^L\\ L_G^S\\ L_G^X\\ L_{\mathrm{\Δ}}^L\end{array}\right)=\left(\begin{array}{c}{V}_P^L\\ V_P^S\\ V_G^L\\ V_G^S\\ V_G^X\\ V_{\mathrm{\Δ}}^L\end{array}\right)$

权阵 $\left(\begin{array}{cccccc}{P}_P^L& & & & & \\ & P_P^S& & & & \\ & & P_G^L& & & \\ & & & P_G^S& & \\ & & & & P_G^X& \\ & & & & & P_{\mathrm{\Δ}}^L\end{array}\right)---\left(27\right)$

误差方程可总体表示为

A·Δ-L＝V    权P    (28)

经法化后可得法方程组为

$N\·\mathrm{\Δ}=\stackrel{\‾}{L}---\left(29\right)$

其中

$N=A^T\mathrm{PA}=\left(\begin{array}{cc}{N}_{11}& N_{12}\\ N_{21}& N_{22}\end{array}\right)$

$N_{11}=\left(\begin{array}{cc}{\left(A_P^L\right)}^T{P}_P^L{A}_P^L+{\left(A_G^L\right)}^T{P}_G^L{A}_G^L+P_{\mathrm{\Δ}}^L& 0\\ 0& {\left(A_P^S\right)}^T{P}_P^S{A}_P^S+{\left(A_G^S\right)}^T{P}_G^S{A}_G^S\end{array}\right)$

$N_{12}={\left(N_{21}\right)}^T=\left(\begin{array}{cc}{\left(A_P^L\right)}^T{P}_P^L{B}_P^L& {\left(A_G^L\right)}^T{P}_G^L{B}_G^L\\ {\left(A_P^S\right)}^T{P}_P^S{B}_P^S& {\left(A_G^S\right)}^T{P}_G^S{B}_G^S\end{array}\right)$

$N_{22}=\left(\begin{array}{cc}{\left(B_P^L\right)}^T{P}_P^L{B}_P^L+{\left(B_P^S\right)}^T{P}_P^S{B}_P^S& 0\\ 0& {\left(B_G^L\right)}^T{P}_G^L{B}_G^L+{\left(B_G^S\right)}^T{P}_G^S{B}_G^S+P_G^X\end{array}\right)$

$\stackrel{\‾}{L}=A^T\mathrm{PL}=\left(\begin{array}{c}{\left(A_P^L\right)}^T{P}_P^L{L}_P^L+{\left(A_G^L\right)}^T{P}_G^L{L}_G^L+P_{\mathrm{\Δ}}^L{L}_{\mathrm{\Δ}}^L\\ {\left(A_P^S\right)}^T{P}_P^S{L}_P^S+{\left(A_G^S\right)}^T{P}_G^S{L}_G^S\\ {\left(B_P^L\right)}^T{P}_P^L{L}_P^L+{\left(B_P^S\right)}^T{P}_P^S{L}_P^S\\ {\left(B_G^L\right)}^T{P}_G^L{L}_G^L+{\left(B_G^S\right)}^T{P}_G^S{L}_G^S+P_G^X{L}_G^X\end{array}\right)$

答解(29)式可得

$\mathrm{\Δ}=N^{-1}\·\stackrel{\‾}{L}$

(2)权的确定

由(27)式可知，误差方程组中包含有六类方程，其对应的权分别为：线阵CCD光学影像上加密点的权

和控制点的权

SAR影像上加密点的权

和控制点的权

控制点观测值的权

以及线阵CCD光学影像“伪观测值”的权

下面分别介绍各个权值的确定方法。

1) $P_P^L$

与像点的量测精度有关，而像点的量测精度主要由线阵CCD光学影像的地面分辨率来决定。如果参加平差影像的地面分辨率相同，则权值相同。否则的话，根据权的定义，影像j的权值P_j可按下式确定

$P_j=(\frac{r_{\max x}}{r_{\mathrm{jx}}}\·\frac{r_{\max y}}{r_{\mathrm{jy}}})\·P_{\max }---\left(30\right)$

其中P_max、r_maxx、r_maxy分别为区域网中所有线阵CCD光学影像中分辨率最高的影像的权值(设为1)及其x方向和y方向的分辨率，r_jx、r_jy分别为影像j的x方向和y方向分辨率。

2) $P_G^L$

$P_G^L=k^L\·P_P^L---\left(31\right)$

其中k^L为放大系数，通常大于1，以保证控制点发挥控制作用。

3) $P_P^S$

由(21)式可知，V^S并非像点观测值的直接改正数，而是间接改正数，因此其权值

可以由像点观测值的权值计算得到。按照权与协方差的关系，可以得到

$P_P^S={\left(Q_P^S\right)}^{-1}={\left(\frac{1}{{\mathrm{\σ}}_0^2}\mathrm{\Σ}{v}^S\right)}^{-1}---\left(32\right)$

其中σ₀为单位权中误差。

根据C·Vs^S＝V^S以及协方差的传播规律可得

$\mathrm{\Σ}{v}^S=C\·\mathrm{\Σ}{\mathrm{Vs}}^S\·C^T=C\·\left({\mathrm{\σ}}_0^{2}E\right)\·C^T={\mathrm{\σ}}_0^2\·{\mathrm{CC}}^T---\left(33\right)$

将(33)式代入(32)式得

$P_P^S={\left({\mathrm{CC}}^T\right)}^{-1}---\left(34\right)$

4) $P_G^S$

$P_G^S=k^S\·P_P^S---\left(35\right)$

其中k^S为放大系数，通常大于1，以保证控制点发挥控制作用。

5) $P_G^X$

该值可根据控制点与加密点的测量精度的比例来确定，通常在[0，1]区间内取值。

6) $P_{\mathrm{\Δ}}^L$

对于(18)式中具有非零近似值的外方位元素，

中的相应权值可取接近于0的小值或直接设为0。那些接近0的外方位元素，

中的相应权值可取接近于1的值。

(3)多片前方交会

如果线阵CCD光学影像和SAR影像的定向参数已知，则此时(18)和(20)式可表示为

$\left.\begin{array}{c}{B}_{11}^L\mathrm{dX}+B_{12}^L\mathrm{dY}+B_{13}^L\mathrm{dZ}-l_x^L=v_x^L\\ B_{21}^L\mathrm{dX}+B_{22}^L\mathrm{dY}+B_{23}^L\mathrm{dZ}-l_y^L=v_y^L\end{array}\right\}---\left(36\right)$

$\left.\begin{array}{c}{B}_{11}^S\mathrm{dX}+B_{12}^S\mathrm{dY}+B_{13}^S\mathrm{dZ}-l_x^S=C_{11}{v}_x^S+C_{12}{v}_y^s\\ B_{21}^S\mathrm{dX}+B_{22}^S\mathrm{dY}+B_{23}^S\mathrm{dZ}-l_y^S=C_{21}{v}_x^S+C_{22}{v}_y^S\end{array}\right\}---\left(37\right)$

即

B^L·X^L-L^L＝V^L    权阵P^L    (38)

B^S·X^S-L^S＝V^S    权阵P^S    (39)

(38)和(39)式是多片前方交会的基本误差方程。与传统的多片前方交会不同的是，传统的多片前方交会中误差方程只有一种形式，而这里的多片前方交会中每个像点需要根据其所在像片的类型来判断其误差方程形式，最后形成的误差方程组是由多种不同形式的误差方程混合组成的。

在给定地面点坐标初值后，通过迭代计算即可得到地面坐标X，Y，Z。地面点坐标的初值可以直接设为0，或是取地面控制点的平均值。

3、具体实施步骤

如图4为多源遥感影像联合区域网平差的流程图。

1.区域网概算

区域网概算的工作是设置未知数的初值。其中像片的外方位元素初值可以利用单张像片空间后方交会的方法计算，然后利用像片外方位元素初值可通过多片前方交会得到加密点地面坐标的初始值。

2.区域网平差

(1)根据各个像片的外方位元素和加密点地面坐标的初始值，利用(3)和(16)式计算各个像点的像点坐标；

(2)根据像点所在像片的类型按(18)和(20)式分别建立各加密点和控制点的误差方程式，并根据(26)式对线阵CCD光学影像增加“伪观测值”误差方程式；

(3)根据(30)-(35)式建立权矩阵，并对误差方程组进行法化；

(4)答解法方程组，求出所有未知数的近似值，并改正前一次的近似值；

(5)重复(1)-(4)步，直到各片外方位元素的改正数和各点地面坐标改正数小于规定的限差为止；

(6)采用多片前方交会，计算未参加区域网平差各点的地面坐标。

具体实施例一

附表说明

表1实验影像的参数统计；

表2参与平差影像的控制点数目以及外方位元素初值精度；

表3所有平差方案的加密点数目和加密点地面坐标初值精度；

表4实验一中所有平差方案的结果。

选择了四颗卫星在不同时间获取的北京地区的遥感影像，其具体参数见表1。图5显示了四幅影像的覆盖区域以及地面控制点的分布情况，其中三角形表示的是控制点，共50个，蓝色三角形表示的是光学影像上的控制点，共27个，红色三角形表示的是SAR影像上的控制点，共23个，绿色圆形表示的是检查点，共35个。这里将检查点作为平差时的加密点。

表1实验影像的参数统计

对四幅影像进行不同的组合可以产生11种平差方案(见图6)，包括两幅影像构成的6种平差方案、三幅影像构成的4种平差方案以及四幅影像构成的1种平差方案。每幅影像使用其各自的控制点，且外方位元素初值由各片的空间后方交会计算，平差中加密点地面坐标的初值由多片前方交会确定。各个影像的控制点数目以及定向精度列于表2，各个平差方案的加密点数目以及加密点地面坐标初值的精度列于表3。

在实验中，我们对11个平差方案分别进行了计算，并通过实验结果分析不同组合方案的平差性能及特点。11个平差方案的平差结果列于表4。从表4的结果来看，(2)、(5)、(6)、(8)四个方案的平差结果较好，其中影像外方位元素的精度与初值很接近，加密点地面坐标的平面精度为20-30米，高程精度为25-33米。

实验结果表明，采用本发明对多幅线阵CCD光学和SAR遥感影像进行大面积的区域网整体平差是正确可行的。

表2参与平差影像的控制点数目以及外方位元素初值精度

表3所有平差方案的加密点数目和加密点地面坐标初值精度

表4实验一中所有平差方案的结果

Claims (4)

1.线阵CCD光学与SAR影像联合区域网平差定位方法，其特征在于，步骤如下：

1)对于同一区域两幅及两幅以上线阵CCD光学和/或SAR遥感图像，任意选择两幅或两幅以上像片；

2)每幅像片使用其各自的控制点，根据每幅像片对应的构像模型，通过空间后方交会确定每幅像片的外方位元素初值，然后由选择的所有像片通过前方交会确定所有加密点地面坐标的初值；

3)将步骤2)得到的各加密点地面坐标、各像片的外方位元素的初值带入各像片分别对应的构像模型，得到加密点在各像片上对应像点的坐标；

4)根据加密点和控制点所在像片类型的构像模型，建立加密点和控制点的误差方程式并进行解算，得到各像片的外方位元素和加密点的地面坐标的近似值；

5)重复3)——4)，直到连续两次得到的各像片的外方位元素改正数和加密点地面坐标改正数小于设定的限差；

6)采用多片前方交会，计算未参加区域网平差各点的地面坐标，实现定位。

2.根据权利要求1所述的线阵CCD光学与SAR影像联合区域网平差定位方法，其特征在于，所述线阵CCD光学影像采用线中心投影构像模型；所述SAR影像采用F.Leberl构像模型。

3.根据权利要求2所述的线阵CCD光学与SAR影像联合区域网平差定位方法，其特征在于，步骤4)包括根据所述两构像模型的联合区域网平差模型建立误差方程式和相应的权矩阵，接着按照最小二乘原理对误差方程进行法化答解。

4.根据权利要求3所述的线阵CCD光学与SAR影像联合区域网平差定位方法，其特征在于，步骤4)中，如果是线阵CCD光学影像，增加伪观测值误差方程式。

Images