CN105466400A

CN105466400A - 一种利用rpc探测多源卫星影像同名像点成像交会角方法

Info

Publication number: CN105466400A
Application number: CN201510940824.1A
Authority: CN
Inventors: 谢国靖; 陈柱; 杨磊
Original assignee: Guizhou Tdt Technology Co Ltd
Current assignee: Guizhou Tdt Technology Co Ltd
Priority date: 2015-12-16
Filing date: 2015-12-16
Publication date: 2016-04-06
Anticipated expiration: 2035-12-16
Also published as: CN105466400B

Abstract

一种利用RPC探测多源卫星影像同名像点成像交会角方法，是一种基于有理多项式(RPC)参数和全球SRTM数据、适合同源或异源卫星遥感影像同名像点成像交会角的探测方法，在全球SRTM数据的支持下，通过高程的改变，利用RPC模型来计算地面点位置矢量的变化，进而计算地面目标的入射角，以及任意两景影像同名像点地面目标的成像交会角。本发明在不知晓传感器位置的情况下，利用RPC参数可以获得重叠影像同名点的成像交会角，为RPC模型定位成像交会角探测提供一种有效的方法，能够在影像几何处理中预防因弱交会成像而导致的地面点坐标求解病态问题提供基础信息。

Description

一种利用RPC探测多源卫星影像同名像点成像交会角方法

技术领域

本发明涉及卫星影像摄影测量方法，特别涉及到利用有理多项式模型进行空中三角测量或三维重建时探测多源卫星影像同名像点成像交会角探测方法。

背景技术

基高比或成像交会角反映了影像像点对应三维地面点交会求解的稳定性。基高比是摄影基线长度与摄站高度之比。摄影基线是多站摄影时，遥感器连续两次曝光瞬间镜头中心间的距离，一般由曝光时间间隔来控制；摄站高度系镜头中心沿铅垂线到地面的距离，即相对航高。基高比决定相邻影像在航向上的重叠度。基高比愈大，航向重叠愈小，反之重叠度愈大。而影像航向重叠愈大，立体观测和高程量测精度愈低。重叠度大，相邻两幅影像组成立体像对时构成的立体模型，其垂直比例尺夸大不明显，亦即立体效应差，影响观测。成像交会角是地面点与立体遥感影像成像时相应传感器位置间的夹角，交会角的大小决定了遥感影像目标对地三维定位时的精度与稳健性。基高比和成像交会角是判断经典光学遥感影像进行定位质量的重要参数，满足特定基高比或成像交会角要求也是经典影像进行摄影测量制图时需要满足的必须条件。

RPC(有理多项式参数)模型目前正成为星载遥感影像定位的主流方法，但RPC参数中不包含传感器的位置，使得传感器至地面目标的光线或电磁波入射向量难以使用传统的方法计算。因此，利用RPC参数进行成像交会角的探测，对RPC进行影像几何处理，无疑有着重要的作用。

卫星影像在成像的过程中，受到透视投影、摄影轴倾斜、大气折光、地球曲率及地形起伏等诸多因素影响，致使影像中各像点产生不同程度的几何变形而失真，无法直接与真实地图套合以进行后续应用。因此要精确地使用遥感卫星影像，首要的工作当先解决影像之传感器的几何模型问题。通常建立传感器的模型的方法主要可分成两大类，除严密模型外，通用模型在近几年受到前所未有的关注，其主要进展是有理多项式模型逐渐替代传统的严密模型，用于光学与雷达卫星遥感影像的几何处理中。

常用的多项式模型比较，RPC模型实际上是各种传感器几何模型的一种抽象的表达方式，它适用于各类传感器包括最新的航空和航天传感器，是多项式模型更精确的形式。它以像点坐标(R,C)表达为以相应地面点空间坐标(L,B,H)为自变量的多项式比值：

\{\begin{matrix} R_{n} = \frac{N_{u m R} (L_{n}, B_{n}, H_{n})}{D_{e n R} (L_{n}, B_{n}, H_{n})} \\ C_{n} = \frac{N_{u m C} (L_{n}, B_{n}, H_{n})}{D_{e n C} (L_{n}, B_{n}, H_{n})} \end{matrix}

其中：

N_umR(L_n,B_n,H_n)＝a₁+a₂B_n+a₃L_n+a₄H_n+a₅B_nL_n+a₆B_nH_n+a₇L_nH_n+a₈B_n ²+a₉L_n ²+a₁₀H_n ²+a₁₁B_nL_nH_n+a₁₂B_n ³+a₁₃B_nL_n ²+a₁₄B_nH_n ²+a₁₅B_n ²L_n+a₁₆L_n ³+a₁₇L_nH_n ²+a₁₈B_n ²H_n+a₁₉L_n ²H_n+a₂₀H_n ³

D_enR(L_n,B_n,H_n)＝b₁+b₂B_n+b₃L_n+b₄H_n+b₅B_nL_n+b₆B_nH_n+b₇L_nH_n+b₈B_n ²+b₉L_n ²+b₁₀H_n ²+b₁₁B_nL_nH_n+a₁₂B_n ³+a₁₃B_nL_n ²+b₁₄B_nH_n ²+b₁₅B_n ²L_n+b₁₆L_n ³+b₁₇L_nH_n ²+b₁₈B_n ²H_n+b₁₉L_n ²H_n+b₂₀H_n ³

N_umC(L_n,B_n,H_n)＝c₁+c₂B_n+c₃L_n+c₄H_n+c₅B_nL_n+c₆B_nH_n+c₇L_nH_n+c₈B_n ²+c₉L_n ²+c₁₀H_n ²+c₁₁B_nL_nZ_n+c₁₂B_n ³+c₁₃B_nL_n ²+c₁₄B_nH_n ²+c₁₅B_n ²L_n+c₁₆L_n ³+c₁₇L_nH_n ²+c₁₈B_n ²H_n+c₁₉L_n ²H_n+c₂₀H_n ³

D_enC(L_n,B_n,H_n)＝d₁+d₂B_n+d₃L_n+d₄H_n+d₅B_nL_n+d₆B_nH_n+d₇L_nH_n+d₈B_n ²+d₉L_n ²+d₁₀H_n ²+d₁₁B_nL_nH_n+d₁₂B_n ³+d₁₃B_nL_n ²+d₁₄B_nH_n ²+d₁₅B_n ²L_n+d₁₆L_n ³+d₁₇L_nH_n ²+d₁₈B_n ²H_n+d₁₉L_n ²H_n+d₂₀H_n ³

式中(L_n,B_n,H_n)为地面点(L,B,H)的归一化物方坐标系坐标，N_umR(L_n,B_n,H_n),D_enR(L_n,B_n,H_n),N_umC(L_n,B_n,H_n),D_enC(L_n,B_n,H_n)为三元三次多项式,a_i，b_i，c_i，d_i(i＝1….20)为影像附带的RPC参数，b₁＝1，d₁＝1；(R_n,C_n)为像点影像坐标(R,C)的归一化影像坐标，根据下式计算：

\{\begin{matrix} R_{n} = (R - R_{o f f}) / R_{s c a l e} \\ C_{n} = (C - C_{o f f}) / C_{s c a l e} \end{matrix}

\{\begin{matrix} L_{n} = (L - L_{o f f}) / L_{s c a l e} \\ B_{n} = (B - B_{o f f}) / B_{s c a l e} \\ H_{n} = (H - H_{o f f}) / H_{s c a l e} \end{matrix}

其中R_off,C_off,L_off,B_off,H_off是像点的影像坐标和物方大地坐标的偏移参数，R_scale,C_scale,L_scale,B_scale,H_scale为像点和地面点坐标的比例系数。

与多项式模型一样，RPC模型参数的计算本质上是通过采样点来拟合。这里的采样点是计算RPC时已知的像点和地面点坐标的集合：

{R_i,C_i,L_i,B_i,H_i}(i＝1,2,3…为第i个已知点)

根据已知点的来源，RPC模型参数求解有与地形无关和与地形相关两种求解方式。在已有成像几何模型的情况下，采用与地形无关的求解方式，通过已有模型(通常是严密模型)产生虚拟的控制点；在严密模型未知的情况下，则采用与地形相关的求解方式，通过测量给定一定数目的控制点。

RPC模型在数据处理方式上虽然有与传感器和坐标系无关的特点。目前，已有大量文献对不同传感器影像RPC模型的精度进行了测试与应用，在一些高质量的高分影像如IKONOS和QuickBird中，通过平差定位精度可达到子像素级。

由于RPC模型属于传感器无关模型，它非常适合不同卫星影像的联合处理，在传感器摄影时摄站位置未知的情况下，探测地面目标的成像交会角和不同影像的弱交会条件，有重要的作用，但目前未出现有效的方法，本发明正是为解决这一难题而提出。

发明内容

本发明提供了一种利用RPC探测多源卫星影像同名像点成像交会角方法，是一种基于RPC参数和全球SRTM(航天飞机雷达地形测绘)数据的多源影像成像交会角的探测方法，在全球SRTM数据的支持下，通过地面点高程的变动，来计算地面点位置矢量的改变，进而计算任意影像像点对应的地面入射角，根据入射角计算任意两景影像同名像点对应的地面目标的成像交会角，并根据影像的性质判断多源影像的同名点是否属于弱交会成像。该方法通过以下步骤实现：

(1)利用卫星影像附带的RPC参数，高程采用全球SRTM数据，利用RPC模型计算各影像像点p(R,C)对应地面点P₀的大地坐标经度、纬度和高程(L₀,B₀,H₀)，R、C为像点p的影像行、列坐标；

(2)依次取新的高程值H₁＝H₀+1和H₂＝H₀+1000，利用RPC模型重新计算得到p点对应地面点P₁和P₂的大地坐标经度、纬度和高程(L₁,B₁H₁)、(L₂,B₂,H₂)；

(3)将P₀、P₁和P₂的大地坐标转换为地心直角坐标系三维坐标，分别得到P₀(X₀，Y₀，Z₀)、P₁(X₁，Y₁，Z₁)、P₂(X₂，Y₂，Z₂)，并计算地心直角坐标系中的向量和向量

(4)若p为光学卫星影像上的像点，对应地面点的入射向量若p为合成孔径雷达影像上的像点，入射向量为

(5)计算任意重叠影像1和2上的同名像点对应的成像交会角β＝arccos[α₁·α₂/(|α₁|·|α₂|)]，其中α₁、α₂分别按步骤(4)计算得到的影像1和2中同名像点对应地面点的入射向量；

(6)根据影像性质和成像交会角大小，判定同名点是否属于弱交会成像。

所述步骤(1)中，影像坐标为(R,C)的像点p对应地面点大地坐标P₀(L₀，B₀，H₀)的计算根据影像RPC参数和SRTM数据，利用RPC模型通过迭代的方法计算得到，包括以下步骤：

(1)给定影像点对应的地面点P₀所在区域的近似高程h₀作为初值，即令H₀＝h₀；

(2)将像点坐标(R,C)、高程H₀和RPC参数代入RPC模型，得到述方程组并求解地面点的归一化大地坐标经度和纬度(L_n,B_n)：

\{\begin{matrix} R_{n} = \frac{N_{u m R} (L_{n}, B_{n}, H_{n})}{D_{e n R} (L_{n}, B_{n}, H_{n})} \\ C_{n} = \frac{N_{u m C} (L_{n}, B_{n}, H_{n})}{D_{e n C} (L_{n}, B_{n}, H_{n})} \end{matrix}

其中，N_umR(L_n,B_n,H_n),D_enR(L_n,B_n,H_n),N_umC(L_n,B_n,H_n),D_enC(L_n,B_n,H_n)代表三元三次多项式，且有：

式中a_i，b_i，c_i，d_i(i＝1….20)为影像附带的RPC参数；R_n为影像行坐标R的归一化值，C_n为影像列坐标C的归一化值，H_n为地面点高程H₀的归一化值，它们根据下式计算：

\{\begin{matrix} R_{n} = (R - R_{o f f}) / R_{s c a l e} \\ C_{n} = (C - C_{o f f}) / C_{s c a l e} \\ H_{n} = (H_{0} - H_{o f f}) / H_{s c a l e} \end{matrix}

其中，R_off、C_off、H_off为像点影像行向坐标、列向坐标和地面点高程的平移参数，R_scale、C_scale、H_scale为像点影像行向坐标、列向坐标和地面点高程的缩放系数，它们均取自影像附带的RPC参数值；

(3)在得到像点的归一化坐标(L_n,B_n)后，根据下式求解地面点的大地坐标(L₀,B₀)：

\{\begin{matrix} L_{0} = L_{n} L_{s c a l e} + L_{o f f} \\ B_{0} = B_{n} B_{s c a l e} + B_{o f f} \end{matrix}

上式中，L_off、B_off为地面点大地坐标经度和纬度的平移参数，L_scale、B_scale为地面点大地坐标经度和纬度的缩放系数，它们均取自影像附带的RPC参数；

(4)将大地坐标(L₀，B₀，H₀)转换到横轴墨卡托地图投影坐标系中，得到坐标(x₀，y₀)，并从SRTM中内插出该投影坐标对应的高程值h，给定阈值10，若|h-H₀|<10，则(L₀，B₀，H₀)即为P₀点的大地坐标值，终止迭代，否则，取H₀＝h，重回步骤(2)。

所述步骤(3)中，P₀、P₁和P₂的大地坐标转换到地心直角坐标系中，坐标转换时所采用的椭球参数为WGS-84椭球参数，转换按下述公式进行：

e = \sqrt{\frac{a^{2} - b^{2}}{a^{2}}}

N = \frac{a}{\sqrt{1 - e^{2} \sin^{2} B}}

[\begin{matrix} X \\ Y \\ Z \end{matrix}] = [\begin{matrix} (N + H) \cos B \cos L \\ (N + H) \cos B \sin L \\ (N (1 - e^{2}) + H) \sin B \end{matrix}]

上式中，(X,Y,Z)为地面点在地心直角坐标系中的坐标，(L,B,H)为大地坐标的经度、纬度和高程，N为地球椭球的卯酉圈曲率半径，a,b,e分别为地球椭球的长半轴、短半轴和第一偏心率，其中a＝6378137，b＝6356752.3142；

所述步骤(6)中，弱交会判断方法为首先给定阈值σ，并根据影像性质进行判定：

(1)当同名点所在影像对同为光学影像，或同为合成孔径雷达影像时：若交会角|β|<σ时，认为同名点属于弱交会成像，否则认为不属于弱交会成像；

(2)当同名点所在影像对一景为光学影像，另一景为合成孔径雷达影像时：若交会角满足|β-90|≤σ时，认为同名点属于弱交会成像，否则认为不属于弱交会成像。

RPC模型定位成像交会角的探测方法，其特征是在RPC参数和SRTM数据的支持下，通过改变地面点的高程值，计算地面点的变化向量，并利用变化向量计算光学影像和雷达影像的对应地面点入射向量，以及任意同源或异源影像同名像点对应目标的成像交会角。

附图说明

图1是本发明提供的利用RPC探测多源卫星影像同名像点成像交会角方法的示意图。

具体实施方式

所述利用RPC探测多源卫星影像同名像点成像交会角方法是在全球SRTM数据的支持下，通过地面点高程的变动，计算地面点位置矢量的改变，进而计算任意影像像点对应的地面入射角，以及任意两景影像同名像点对应的地面目标的成像交会角。如图1所示，该方法通过以下步骤实现：

步骤1，利用卫星遥感影像附带的RPC参数，高程采用全球SRTM数据，利用RPC模型计算各像点p(R,C)对应地面点P₀的大地坐标经度、纬度和高程(L₀，B₀，H₀)；大地坐标P₀(L₀，B₀，H₀)的计算根据RPC参数和SRTM数据，利用RPC模型通过迭代的方法计算得到，包括以下步骤：

1)在SRTM中，给定P₀所在区域中附近的近似高程h₀，并令H₀＝h₀；

2)利用像点影像坐标(R,C)和RPC参数构建有理多项式模型，并按下式所构建的方程组求解物方坐标经纬度的归一化坐标(L_n,B_n)：

\{\begin{matrix} R_{n} = \frac{N_{u m R} (L_{n}, B_{n}, H_{n})}{D_{e n R} (L_{n}, B_{n}, H_{n})} \\ C_{n} = \frac{N_{u m C} (L_{n}, B_{n}, H_{n})}{D_{e n C} (L_{n}, B_{n}, H_{n})} \end{matrix}

其中N_umR(L_n,B_n,H_n),D_enR(L_n,B_n,H_n),N_umC(L_n,B_n,H_n),D_enC(L_n,B_n,H_n)为三元三次多项式,：

式中a_i，b_i，c_i，d_i(i＝1….20)为影像附带的RPC参数；(R_n,C_n)为像点影像坐标(R,C)的归一化坐标，H_n为地面点高程H的归一化坐标，它们根据下式计算：

\{\begin{matrix} R_{n} = (R - R_{o f f}) / R_{s c a l e} \\ C_{n} = (C - C_{o f f}) / C_{s c a l e} \\ H_{n} = (H_{0} - H_{o f f}) / H_{s c a l e} \end{matrix}

上式中，(R_off,C_off,H_off)为像点影像坐标和地面点高程的平移参数，(R_scale,C_scale,H_scale)为像点坐标和地面点高程的缩放系数，它们均取自影像附带的RPC参数；

3)在得到(L_n,B_n)后，根据下面公式求解地面点的大地坐标(L₀,B₀)：

\{\begin{matrix} L_{0} = L_{n} L_{s c a l e} + L_{o f f} \\ B_{0} = B_{n} B_{s c a l e} + B_{o f f} \end{matrix}

上式中，(L_off,B_off)为地面点大地坐标经度和纬度的平移参数，(L_scale,B_scale)为地面点大地坐标经度和纬度的缩放系数，它们均取自影像附带的RPC参数；

4)将大地坐标(L₀，B₀，H₀)转换到UTM投影坐标系中，得到坐标(x₀，y₀)，并从SRTM中内插出该平面坐标对应的高程值h，给定阈值10，若|h-H₀|<10，坐标(L₀，B₀，H₀)即为P₀点的坐标，终止迭代，否则，取H₀＝h，重回步骤2)。

步骤2，分别取新的高程值H₁＝H₀+1和H₂＝H₀+1000，利用RPC重新计算像点p对应地面点的大地坐标P₁(B₁,L₁,H₁)和P₂(B₂,L₂,H₂)；

步骤3，将P₀、P₁和P₂的大地坐标转换到地心直角坐标中，坐标转换时所采用的椭球参数为WGS-84椭球参数，转换按下述公式进行：

e = \sqrt{\frac{a^{2} - b^{2}}{a^{2}}}

N = \frac{a}{\sqrt{1 - e^{2} \sin^{2} B}}

[\begin{matrix} X \\ Y \\ Z \end{matrix}] = [\begin{matrix} (N + H) \cos B \cos L \\ (N + H) \cos B \sin L \\ (N (1 - e^{2}) + H) \sin B \end{matrix}]

上式中，(X,Y,Z)为地面点在地心直角坐标系中的坐标，(L,B,H)为大地经纬度和大地高，N为椭球的卯酉圈曲率半径，a,b,e分别为地球椭球的长半轴、短半轴和第一偏心率，其中a＝6378137，b＝6356752.3142；

利用上式分别计算P₀、P₁和P₂的地心直角坐标P₀(X₀，Y₀，Z₀)、P₁(X₁，Y₁，Z₁)、P₂(X₂，Y₂，Z₂)，并计算空间直角坐标系向量和向量

步骤4，计算像点p对应地面点入射向量，对于光学卫星影像，入射向量对于合成孔径雷达影像，入射向量

步骤5，计算任意两个影像1和影像2上的同名像点对应的成像交会角：β＝arccos[α₁·α₂/(|α₁|·|α₂|)]，其中α₁、α₂分别为影像1和影像2中同名像点对应地面点的入射向量；

步骤6，根据影像性质和成像交会角大小，判定同名点是否属于弱交会成像。弱交会成像判断方法为首先给定阈值σ，并根据影像性质进行判定：

通过图1可以看出本发明主要包括以下步骤：(1)利用卫星遥感影像附带的RPC参数，高程采用全球SRTM数据，计算同名影像像点坐标对应物方坐标P₀(L₀,B₀,H₀)；

(2)分别取新的高程值H₁＝H₀+1和H₂＝H₀+1000，利用RPC参数重新计算地面点P₁和P₂的大地坐标；

(3)将P₀、P₁和P₂的大地坐标转换为地心直角坐标系中的坐标，并计算地心直角坐标系中空间向量和向量

(4)根据向量和向量分别针对光学传感器和雷达传感器，计算地面入射向量；

(5)根据任意两个重叠影像上的同名点对应的入射向量，计算重叠影像同名像点的成像交会角；

(6)根据立体影像的性质和成像交会角的大小，判定同名点是否属于弱交会成像。弱交会成像判断方法为首先给定阈值σ，并根据影像性质进行判定：

1)若立体影像同为光学影像，或同为合成孔径雷达影像时：在成像交会角满足|β|<σ时,认为是弱交会成像，否则认为不属于弱交会成像；

2)若立体影像一景为光学影像，另一景为合成孔径雷达影像时：若成像交会角满足|β-90|≤σ,则认为同名点属于弱交会成像，否则认为不属于弱交会成像。

Claims

1.一种利用RPC探测多源卫星影像同名像点成像交会角方法，包括以下步骤：

2.根据权利要求1所述利用RPC探测多源卫星影像同名像点成像交会角方法，其特征在于：所述步骤(1)中，影像坐标为(R,C)的像点p对应地面点大地坐标P₀(L₀，B₀，H₀)的计算根据影像RPC参数和SRTM数据，利用RPC模型通过迭代的方法计算得到，包括以下步骤：

(2)将像点坐标(R,C)、高程H₀和RPC参数代入有理多项式RPC模型，得到下述方程组并求解地面点的归一化大地坐标经度和纬度(L_n,B_n)：

式中a_i，b_i，c_i，d_i(i＝1….20)为取自影像附带的RPC参数数据集；R_n为影像行坐标R的归一化值，C_n为影像列坐标C的归一化值，H_n为地面点高程H₀的归一化值，它们根据下式计算：

其中，R_off、C_off、H_off为像点影像行向坐标、列向坐标和地面点高程的平移参数，R_scale、C_scale、H_scale为像点影像行向坐标、列向坐标和地面点高程的缩放系数，它们均取自影像附带的RPC参数数据集；

3.根据权利要求1所述的利用RPC探测多源卫星影像同名像点成像交会角方法，其特征在于：所述步骤(3)中，P₀、P₁和P₂的大地坐标转换到地心直角坐标系中，坐标转换时所采用的椭球参数为WGS-84椭球参数，转换按下述公式进行：

上式中，(X,Y,Z)为地面点在地心直角坐标系中的坐标，(L,B,H)为大地坐标的经度、纬度和高程，N为地球椭球的卯酉圈曲率半径，a,b,e分别为地球椭球的长半轴、短半轴和第一偏心率，其中a＝6378137，b＝6356752.3142。

4.根据权利要求1所述的利用RPC探测多源卫星影像同名像点成像交会角方法，其特征在于：所述步骤(6)中，弱交会成像判断方法为首先给定阈值σ，并根据影像性质进行判定：