多角度多时段导航卫星双基地PS-InSAR三维形变反演方法

技术领域

本发明属于雷达信号处理的技术领域，具体涉及一种多角度多时段导航卫星双基地PS-InSAR三维形变反演方法。

背景技术

基于导航卫星的双基地重轨差分干涉SAR利用导航卫星作为发射平台，地面配置接收机，对需要探测形变的区域进行连续成像、差分干涉和形变反演处理可获得监测区域的一维形变。

然而，单颗导航卫星的重轨时间较长，同时仅能实现一维形变的反演，精度受限，实用性较差。实际上，导航卫星系统由多颗卫星组成，对于同一地区始终有多颗卫星的信号连续覆盖，而且随着时间的不同，卫星相对于一个特定场景的几何关系也会发生变化，即可提供不同视角的信息。因此，该系统在提高形变监测精度方面具有很大潜力。

因此，开发一种基于多角度多时段导航卫星照射的双基地PS-InSAR的三维连续形变监测方法，对于高精度形变测量领域具有重要意义。

发明内容

有鉴于此，本发明提供了一种多角度多时段导航卫星双基地PS-InSAR三维形变反演方法，能够利用多卫星、多时段的回波数据实现照射场景的三维形变测量，其精度高、成本低且时间连续。

实现本发明的技术方案如下：

步骤一：利用不同角度不同时段下的直达波信号对相应角度相应时段的回波信号进行同步和成像处理，得到同一场景的不同角度的双基地SAR图像序列；

步骤二：对所述双基地SAR图像序列，利用PS方法实现各角度下的PS点识别与PS点的一维形变量的反演；

步骤三：根据所述PS点的一维形变量的反演结果，利用Kriging插值方法对所述双基地SAR图像序列的每个像素的一维形变量进行反演；

步骤四：设置一个统一的时间轴，利用时间序列分析方法对所述双基地SAR图像序列的每个像素的一维形变量的反演结果进行时间插值，得到同一场景的不同角度下的每个像素在统一的时间轴下的一维形变量历史；

步骤五：利用加权最小二乘估计方法，对所述每个像素在统一的时间轴下的一维型变量历史进行处理，实现每个像素的三维形变量历史的反演。

有益效果：

本发明对比已有技术，充分利用了导航卫星数目多、导航信号覆盖范围广、照射时间长的特点，与监测区域附近的接收机配置构成双基地SAR系统，可对监测区域实现时间、空间连续的三维形变测量，能够获取更加全面的形变信息。同时作为被动测量方式，本发明成本低，配置灵活，可广泛应用在形变监测领域。

附图说明

图1为本发明系统构型示意图。

图2为本发明算法总流程图。

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。

本发明方法所利用的系统构型如图1所示，不同卫星从不同角度发射导航信号到场景，场景反射回波到回波天线，同时，一个全向直达波天线接收卫星的直达波信号。对于每颗卫星，采集其不同重轨时段的信号。以接收机直达波天线作为原点O，正东为X轴，正北为Y轴，竖直向上为Z轴建立空间直角坐标系O-XYZ。如图2所示为本发明方法的流程，该方法具体步骤如下：

步骤一：利用不同角度不同时段下的直达波信号对相应角度相应时段的回波信号进行同步和成像处理，得到同一场景的不同角度下的双基地SAR图像序列。设第i颗卫星第j组时段(一组时段对应同一个角度)的M_ij+1幅重轨图像序列为Q_i,j：

$Q_{i,j}={[q_{i,j,0},q_{i,j,1},···,q_{i,j,M_{\mathrm{ij}}}]}^T,i=\mathrm{0,2},···S-1,j=\mathrm{0,1},···K_{i-1}---\left(1\right)$

其中，S表示卫星总数，K_i表示第i颗卫星的时段总数。q_i,j,0表示第i颗卫星第j组时段的第一幅图像，q_i,j,1表示第i颗卫星第j组时段的第二幅图像，依此类推。

第i颗卫星第j组时段的卫星重轨时间序列T_i,j为：

$T_{i,j}={[t_{i,j,0},t_{i,j,1},···,t_{i,j,M_{\mathrm{ij}}}]}^T---\left(2\right)$

其中，t_i,j,0表示第i颗卫星第j组时段的第一个重轨时刻，t_i,j,1表示第i颗卫星第j组时段的第二个重轨时刻，其他依此类推。

由于不同时段，卫星相对于观测场景的观测角度不同，因此，所有卫星所有时段的观测角度总数L为：

$L=\underset{i=0}{\overset{S-1}{\Sigma }}{K}_i---\left(3\right)$

步骤二：对第i颗卫星第j组时段的M_ij+1幅重轨图像进行PS点识别与PS点的一维形变量反演，得到各个PS点沿着等效视线方向的一维形变反演结果ΔL_i,j(P_i,j；t_i,j,0)，设为(设q_i,j,0为主图像)：

${\mathrm{\Delta L}}_{i,j}(P_{i,j};t_{i,j,0})=\left(\begin{array}{c}{\mathrm{\Delta l}}_{i,j}(P_{i,j};t_{i,j,0},t_{i,j,1})\\ {\mathrm{\Delta l}}_{i,j}(P_{i,j};t_{i,j,0},t_{i,j,2})\\ .\\ .\\ .\\ {\mathrm{\Delta l}}_{i,j}(P_{i,j};t_{i,j,0},t_{i,j,M_{ij}})\end{array}\right)---\left(4\right)$

其中，P_i,j为第i颗卫星的第j组图像中的PS点集合。Δl_i,j(P_i,j；t_i,j,0,t_i,j,1)表示第i颗卫星的第j组图像中PS点沿着等效视线方向的t_i,j,1时刻相对于t_i,j,0时刻的形变反演结果。Δl_i,j(P_i,j；t_i,j,0,t_i,j,2)表示第i颗卫星的第j组图像中PS点沿着等效视线方向的t_i,j,2时刻相对于t_i,j,0时刻的形变反演结果。其他依此类推。

步骤三：由于不同卫星不同角度下探测到的PS点不同，无法直接计算某个像素的三维形变量，因此，根据不同角度下SAR图像中各个PS点的一维形变反演结果，利用Kriging插值方法，对每个角度下整个场景的双基地SAR图像序列的每个像素进行空间插值，得到空间插值后的每个角度下场景中每个像素的形变量矩阵：

${\mathrm{\Delta L}}_{i,j}(\tilde{P};t_{i,j,0})=\left(\begin{array}{c}{\mathrm{\Delta l}}_{i,j}(\tilde{P};t_{i,j,0},t_{i,j,1})\\ {\mathrm{\Delta l}}_{i,j}(\tilde{P};t_{i,j,0},t_{i,j,2})\\ .\\ .\\ .\\ {\mathrm{\Delta l}}_{i,j}(\tilde{P};t_{i,j,0},t_{i,j,M_{ij}})\end{array}\right)---\left(5\right)$

其中，表示场景中的像素点集合。表示空间插值后的Δl_i,j(P_i,j；t_i,j,0,t_i,j,1)，表示空间插值后的Δl_i,j(P_i,j；t_i,j,0,t_i,j,2)，其他依此类推。

步骤四：由于卫星的重轨时刻未必完全相同，为实现场景中每个像素三维形变量历史的提取，需要将所获取的不同角度下的形变量历史的时间轴进行统一。设统一时间轴为：

t＝[t₀,t₁,…,t_N-1]>

其中，t₀为时间轴的第1个时刻，t_k＝t₀+k·T(k＝1,2,…,N-1)为该时间轴的第k+1个时刻，其中，T表示统一时间轴的时间间隔，N表示统一时间轴的时刻总数。

那么，对进行基于时间序列分析方法的时间插值得到空间、时间插值后的形变量矩阵：

$\Delta {\tilde{L}}_{i,j}(\tilde{P};t)=\left(\begin{array}{c}\Delta {\tilde{l}}_{i,j}(\tilde{P};t_0)\\ \Delta {\tilde{l}}_{i,j}(\tilde{P};t_1)\\ .\\ .\\ .\\ \Delta {\tilde{l}}_{i,j}(\tilde{P};t_{N-1})\end{array}\right)---\left(7\right)$

其中，表示时间插值后的形变量在t₀时刻的形变量，其他依此类推。

步骤五：

设在第t_m(m＝0,1,…N-1)时刻，场景中的点的三维形变量为：

$\Delta D(\tilde{p};t_m)={\left(\begin{array}{ccc}\Delta x(\tilde{p};t_m)& \Delta y(\tilde{p};t_m)& \Delta z(\tilde{p};t_m)\end{array}\right)}^T---\left(8\right)$

其中，和分别为点在第t_m时刻沿着前述坐标系O-XYZ上的X、Y和Z轴方向的形变量。

则有：

$\Delta D{(\tilde{p};t_m)}^T{g}_{i,j}\left(\tilde{p}\right)=\Delta {\tilde{l}}_{i,j}(\tilde{p};t_m)+w_{i,j}(\tilde{p};t_m)---\left(9\right)$

其中，为噪声，标准差为为点在第i颗卫星第j组时段的形变量进行时间插值和空间插值后在时刻t_m的形变量，为等效视线方向的向量，可写为：

$g_{i,j}\left(\tilde{p}\right)=\frac{\left(\Psi \right(S_{ij})-\tilde{p})}{|\Psi \left(S_{ij}\right)-\tilde{p}|}+\frac{\left(\Psi \right(R)-\tilde{p})}{|\Psi \left(R\right)-\tilde{p}|}---\left(10\right)$

其中：

Ψ(S_ij)＝[x_ij,y_ij,z_ij]^T>

表示卫星i在第j组时段的孔径中心时刻的三维位置，其中x_ij、y_ij和z_ij分别表示卫星i在第j组时段的孔径中心时刻在前述坐标系O-XYZ上的X、Y和Z轴方向的坐标；

Ψ(R)＝[x_R,y_R,z_R]^T>R、y_R和z_R分别表示回波天线R在前述坐标系O-XYZ上的X、Y和Z轴方向的坐标；

$\tilde{p}={[x_{\tilde{p}},y_{\tilde{p}},z_{\tilde{p}}]}^T---\left(13\right)$

为点的三维位置。其中和分别表示点在前述坐标系O-XYZ上的X、Y和Z轴方向的坐标；

那么，点处三维形变量的加权最小二乘估计为：

$\Delta \hat{D}(\tilde{p};t_m)={\left[G^T\left(\tilde{p}\right)\Phi G\left(\tilde{p}\right)\right]}^{-1}{G}^T\left(\tilde{p}\right)\Phi \Delta \tilde{L}(\tilde{p};t_m)---\left(14\right)$

其中：

$\Delta \tilde{L}(\tilde{p};t_m)=\left(\begin{array}{c}\Delta {\tilde{l}}_{0,0}(\tilde{p};t_m)\\ \Delta {\tilde{l}}_{0,1}(\tilde{p};t_m)\\ .\\ .\\ .\\ \Delta {\tilde{l}}_{0,K_0-1}(\tilde{p};t_m)\\ .\\ .\\ .\\ \Delta {\tilde{l}}_{i,0}(\tilde{p};t_m)\\ \Delta {\tilde{l}}_{i,1}(\tilde{p};t_m)\\ .\\ .\\ .\\ \Delta {\tilde{l}}_{i,K_i-1}(\tilde{p};t_m)\\ .\\ .\\ .\\ \Delta {\tilde{l}}_{S-1,0}(\tilde{p};t_m)\\ \Delta {\tilde{l}}_{S-1,1}(\tilde{p};t_m)\\ .\\ .\\ .\\ \Delta {\tilde{l}}_{S-1,K_{S-1}-1}(\tilde{p};t_m)\end{array}\right)---\left(15\right)$

$G\left(\tilde{p}\right)=\left(\begin{array}{c}{g}_{0,0}\left(\tilde{p}\right)\\ g_{0,1}\left(\tilde{p}\right)\\ .\\ .\\ .\\ g_{0,K_0-1}\left(\tilde{p}\right)\\ .\\ .\\ .\\ g_{i,0}\left(\tilde{p}\right)\\ g_{i,1}\left(\tilde{p}\right)\\ .\\ .\\ .\\ g_{i,K_i-1}\left(\tilde{p}\right)\\ .\\ .\\ .\\ g_{S-1,0}\left(\tilde{p}\right)\\ g_{S-1,1}\left(\tilde{p}\right)\\ .\\ .\\ .\\ g_{S-1,K_{S-1}-1}\left(\tilde{p}\right)\end{array}\right)---\left(16\right)$

Φ为权值矩阵，为：

$\Phi =diag\left(\begin{array}{c}1/{\sigma }_{0,0}^2(\tilde{p};t_m)\\ 1/{\sigma }_{0,1}^2(\tilde{p};t_m)\\ .\\ .\\ .\\ 1/{\sigma }_{0,K_0-1}^2(\tilde{p};t_m)\\ .\\ .\\ .\\ 1/{\sigma }_{i,0}^2(\tilde{p};t_m)\\ 1/{\sigma }_{i,1}^2(\tilde{p};t_m)\\ .\\ .\\ .\\ 1/{\sigma }_{i,K_i-1}^2(\tilde{p};t_m)\\ .\\ .\\ .\\ 1/{\sigma }_{S-1,0}^2(\tilde{p};t_m)\\ 1/{\sigma }_{S-1,1}^2(\tilde{p};t_m)\\ .\\ .\\ .\\ 1/{\sigma }_{S-1,K_{S-1}-1}^2(\tilde{p};t_m)\end{array}\right)---\left(17\right)$

其中，diag表示对角矩阵，即：Φ的第一行第一列为第二行第二列为其他依此类推。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。