融合GNSS数据与InSAR技术的桥梁变形监测方法与流程

融合gnss数据与insar技术的桥梁变形监测方法
技术领域
1.本发明属于桥梁变形监测技术领域，具体涉及一种融合gnss数据与insar技术的桥梁变形监测方法。


背景技术：

2.桥梁作为交通基础设施的重要组成部分，已经成为城市经济发展不可缺少的一部分。但桥梁本身重量大、荷载复杂多变，地质条件变化、超载、恶劣天气等因素会对桥梁产生诸多不良影响，给桥梁安全带来隐患，影响行车安全，甚至威胁人们生命财产安全。因此，加强对桥梁的变形监测，对桥梁安全运营具有重要的意义。
3.传统桥梁形变测量方式，如gnss、水准测量等方式，虽然精度较高，但需要现场布设设备，耗费人力物力、成本较高，且监测结果为稀疏点形式，不利于对桥梁开展整体分析。合成孔径雷达干涉测量技术，利用sar影像的相位信息实现对地表变形高精度的监测，其形变监测精度可以达到厘米甚至毫米量级；同时，sar卫星全天时全天候的影像获取能力，以及insar大范围、高密度、长时序、可回溯的技术特点，相比传统gnss、水准测量等方式具有独特的优势。
4.但是，在insar技术应用过程中，获取的往往是相对于未知形变情况的参考点的相对形变量；同时，由于桥梁的变形影响因素较多，其在时间上往往呈现出非线性的特征，而insar技术中通常应用线性形变模型求解形变，导致结果可能会引入较大的误差。


技术实现要素：

5.针对现有技术存在的缺陷，本发明提供一种融合gnss数据与insar技术的桥梁变形监测方法，可有效解决上述问题。
6.本发明采用的技术方案如下：
7.本发明提供一种融合gnss数据与insar技术的桥梁变形监测方法，包括以下步骤：
8.步骤1，确定桥梁监测对象；获取覆盖桥梁监测对象的g幅高分辨率sar影像；
9.其中，g幅高分辨率sar影像为不同监测时刻得到的sar影像；每幅高分辨率sar影像具有成像时间属性、影像后向散射信息和成像轨道信息；
10.步骤2，在各幅高分辨率sar影像中，确定一幅高分辨率sar影像作为sar主影像，其他g
‑
1幅高分辨率sar影像作为sar从影像；
11.步骤3，根据各幅高分辨率sar影像的成像轨道信息和影像后向散射信息，对各幅sar从影像进行旋转缩放采样操作，实现各幅sar从影像与sar主影像的影像配准，由此得到g
‑
1幅配准后的sar从影像；
12.步骤4，对于sar主影像和g
‑
1幅配准后的sar从影像，任意两幅sar影像形成一个干涉像对，由此生成所有干涉像对；
13.对于每个干涉像对，其包含的两幅sar影像共轭相乘，生成一幅干涉相位图；
14.因此，对于所有干涉像对，生成多幅干涉相位图，由此组成干涉相位图集合；
15.步骤5，对于干涉相位图集合中的每幅干涉相位图，进行预处理、掩膜处理和相位解缠处理，得到解缠差分干涉相位图；
16.由此将干涉相位图集合，转化为解缠差分干涉相位图集合；
17.步骤6，对于步骤5得到的解缠差分干涉相位图集合，假设共有m幅解缠差分干涉相位图，分别表示为：解缠差分干涉相位图tu1，解缠差分干涉相位图tu2，...，解缠差分干涉相位图tu
m
；
18.其中：对于任意解缠差分干涉相位图tu
j
，其中，j＝1，2，...，m，对应的干涉像对表示为：sar影像i
j0
和sar影像i
j1
，即：sar影像i
j0
和sar影像i
j1
共轭相乘后，经过预处理、掩膜处理和相位解缠处理后，得到解缠差分干涉相位图tu
j
；
19.对于m幅解缠差分干涉相位图，共对应m个干涉像对，即：不考虑干涉像对包含共同sar影像时，m个干涉像对包含2m幅sar影像；再对2m幅sar影像去重处理，假设最终得到n 1幅sar影像；
20.将n 1幅sar影像，按成像时间排序，依次表示为：sar影像i0，sar影像i1，...，sar影像i
n
，对应的成像时间依次为：t0，t1，...，t
n
；即：sar影像i0的成像时间为t0，sar影像i1的成像时间为t1，...，sar影像i
n
的成像时间为t
n
；
21.因此，将m幅解缠差分干涉相位图中的任意解缠差分干涉相位图，表示为解缠差分干涉相位图tu
j
(i
c
，i
d
)，其中，c＞d，c≠d，c＝1，2，...，n；d＝0，1，2，...，n
‑
1，含义为：解缠差分干涉相位图tu
j
(i
c
，i
d
)，是由sar影像i
c
和sar影像i
d
共轭相乘后，经过预处理、掩膜处理和相位解缠处理后生成，解缠差分干涉相位图tu
j
(i
c
，i
d
)中，任意一个像元的解缠相位差，可直接通过解缠差分干涉相位图读取获得；
22.步骤7，m幅解缠差分干涉相位图和n 1幅sar影像的图像大小完全相同，具有相同数量的像元；
23.对于任意解缠差分干涉相位图中的像元(l，p)，其中，l，p分别为像元的行号和列号，执行步骤7.1
‑
步骤7.3，计算得到像元(l，p)对应的在地理坐标系下的insar雷达视线向形变量相对位移变化序列：
24.步骤7.1，建立以下方程：
[0025][0026]
其中：
[0027]
分别表示sar影像i1相对于sar影像i0在像元(l，p)位置的解缠相位差，sar影像i2相对于sar影像i0在像元(l，p)位置的解缠相位差，...，sar影像i
n
相对于sar影像i0在像元(l，p)位置的解缠相位差；为未知待求值；
[0028]
分别表示解缠差分干涉相位图tu1在像元(l，p)位置的解缠相位差，解缠差分干涉相位图tu2在像元(l，p)位置的解缠相位差，...，解缠差分干涉相
位图tu
m
在像元(l，p)位置的解缠相位差；可通过对应的解缠差分干涉相位图直接读取到；
[0029]
为m*n的系数矩阵a，为已知值，通过以下方式确定：
[0030]
假设对于系数矩阵a中的任意第j行，其中，j＝1，2，...，m，则有以下关系式：
[0031][0032]
为解缠差分干涉相位图tu
j
(i
c
，i
d
)在像元(l，p)位置的解缠相位差；解缠差分干涉相位图tu
j
(i
c
，i
d
)对应sar影像i
c
和sar影像i
d
；
[0033]
如果d＝0，则：在中，令前面的系数a
jc
为1，其他n
‑
1个系数均为0，由此确定了系数a
j1
，...，a
jn
的值；
[0034]
如果d≠0，则：在中，令前面的系数a
jc
为1，令前面的系数a
jd
为
‑
1，其他n
‑
2个系数均为0，由此确定了系数a
j1
，...，a
jn
的值；
[0035]
步骤7.2，求解步骤7.1建立的方程，得到桥梁在像元(l，p)位置在图像坐标系下的insar雷达视线向形变量相对相位变化序列即：桥梁在像元(l，p)位置，相对于t0时刻，分别在t1，...，t
n
时刻在图像坐标系下的insar雷达视线向形变量相对相位变化值；
[0036]
步骤7.3，依据sar成像几何参数和波长信息，将桥梁在像元(l，p)位置在图像坐标系下的insar雷达视线向形变量相对相位变化序列转换为图像坐标系下insar雷达视线向形变量相对位移变化序列并将图像坐标系下insar雷达视线向形变量相对位移变化序列地理编码，得到地理坐标系下的insar雷达视线向形变量相对位移变化序列；表示为：
[0037]
通过以上方式，假设任意解缠差分干涉相位图中，桥梁监测对象覆盖n0个像元，则t1时刻，各个像元在地理坐标系下的insar雷达视线向形变量相对位移变化值，绘制形成桥梁形变图gr1；
[0038]
t2时刻，各个像元在地理坐标系下的insar雷达视线向形变量相对位移变化值，绘制形成桥梁形变图gr2；
[0039]
依此类推
[0040]
t
n
时刻，各个像元在地理坐标系下的insar雷达视线向形变量相对位移变化值，绘制形成桥梁形变图gr
n
；
[0041]
步骤8，在桥梁监测对象布置num个gnss测点；
[0042]
从初始化时刻开始，各个gnss测点按采样间隔，可持续测得测点所在位置的三维
形变位移值；
[0043]
因此，获得t
i
时刻各个gnss测点测得的三维形变位移值，其中，i＝1，2，...，n，采用t
i
时刻各个gnss测点测得的三维形变位移值，校正桥梁形变图gr
i
，得到校正后的桥梁形变图gr’i
，表示t
i
时刻相对于初始化时刻，各个像元在地理坐标系下的insar雷达视线向形变量位移变化值。
[0044]
优选的，步骤5具体为：
[0045]
步骤5.1，预处理：
[0046]
对每幅干涉相位图进行预处理，包括平地相位去除、地形相位去除和滤波处理，得到处理后的差分干涉相位图；
[0047]
步骤5.2，掩膜处理：
[0048]
步骤5.2.1，制作桥梁监测对象掩膜；其中，所述桥梁监测对象掩膜是指：在掩膜图像中，桥梁监测对象的像素点为1，其他所有位置的像素点为0；
[0049]
步骤5.2.2，采用桥梁监测对象掩膜，对每幅处理后的差分干涉相位图进行掩膜处理，得到掩膜处理后的差分干涉相位图；其中，在掩膜处理后的差分干涉相位图中，桥梁监测对象的像素点相位保持不变，其他位置的像素点相位统一变为0；
[0050]
步骤5.3，相位解缠处理：
[0051]
对于每幅掩膜处理后的差分干涉相位图，均采用多种相位解缠算法进行相位解缠运算，由此得到多幅解缠差分干涉相位图，并采用筛选算法，剔除存在解缠误差的差分干涉相位图，保留相位解缠结果正确的解缠差分干涉相位图；
[0052]
将多种相位解缠算法获取的解缠结果正确的解缠差分干涉相位图合并，最终得到多个解缠差分干涉相位图。
[0053]
优选的，步骤5.1中，采用以下方法，进行地形相位去除：
[0054]
获取覆盖桥梁监测对象的高精度dsm数据；
[0055]
基于高精度dsm数据，进行地形相位去除。
[0056]
优选的，步骤7.2中，采用最小二乘或奇异值分解方法，求解步骤7.1建立的方程。
[0057]
优选的，当m≥n，且矩阵a的秩为n，通过最小二乘法求解；
[0058]
当解缠差分干涉相位图集分为l组，并且，矩阵a的秩小于n，为n
‑
l 1，使用奇异值分解方法求解。
[0059]
优选的，步骤8中，具体采用以下方法，对桥梁形变图gr
i
进行校正：
[0060]
步骤8.1，获得t
i
时刻各个gnss测点测得的三维形变位移值；其中，对于任意第z个gnss测点，z＝1，2，...，num，其在t
i
时刻的三维形变位移值包括：d
xi
、d
yi
、d
ui
，分别代表东西向、南北向和垂向上的位移值；
[0061]
步骤8.2，将第z个gnss测点测得的三维形变位移值，采用以下公式，转化为雷达视线向形变测量值：
[0062][0063]
其中：
[0064]
θ为步骤2中sar主影像对应的雷达波局部入射角；
[0065]
α为步骤2中sar主影像的卫星轨道方位角；
[0066]
在t
i
时刻，第z个gnss测点测得的雷达视线向形变测量值；
[0067]
步骤8.3，对于第z个gnss测点，查找桥梁形变图gr
i
，获得第z个gnss测点所在位置的像元的insar雷达视线向形变量相对位移变化值
[0068]
步骤8.4，采用下式，得到校正偏移量
[0069][0070]
步骤8.5，使桥梁形变图gr
i
整体发生形变得到校正后的桥梁形变图gr’i
。
[0071]
本发明提供的融合gnss数据与insar技术的桥梁变形监测方法具有以下优点：
[0072]
本发明的目的在于克服传统gnss测量等桥梁形变监测手段稀疏布设导致的形变监测空间分辨率低、insar技术手段获取的多是相对未知形变情况参考点的形变量以及insar常用的线性形变模型在桥梁形变监测的误差较大的不足，提出一种融合gnss数据与insar技术的桥梁变形监测方法，将gnss的已有测量结果和insar高空间分辨率形变监测能力结合，极大改善了桥梁基础设施形变监测的空间分辨率。
附图说明
[0073]
图1为本发明提供的融合gnss数据与insar技术的桥梁变形监测方法的流程示意图；
[0074]
图2为sar影像的干涉像对组合图；
[0075]
图3为处理后的差分干涉相位图的示意图；
[0076]
图4为解缠结果正确的解缠差分干涉相位图示意图；
[0077]
图5是桥梁insar形变分布图；
[0078]
图6是系数矩阵确定方式的原理图。
具体实施方式
[0079]
为了使本发明所解决的技术问题、技术方案及有益效果更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
[0080]
本发明的目的在于克服传统gnss测量等桥梁形变监测手段稀疏布设导致的形变监测空间分辨率低、insar技术手段获取的多是相对未知形变情况参考点的形变量以及insar常用的线性形变模型在桥梁形变监测的误差较大的不足，提出一种融合gnss数据与insar技术的桥梁变形监测方法，将gnss的已有测量结果和insar高空间分辨率形变监测能力结合，极大改善了桥梁基础设施形变监测的空间分辨率。
[0081]
参考图1，本发明提供的融合gnss数据与insar技术的桥梁变形监测方法，包括以下步骤：
[0082]
步骤1，确定桥梁监测对象；获取覆盖桥梁监测对象的g幅高分辨率sar(synthetic aperture radar，即合成孔径雷达)影像；
[0083]
其中，g幅高分辨率sar影像为不同监测时刻得到的sar影像；每幅高分辨率sar影像具有成像时间属性、影像后向散射信息和成像轨道信息；
[0084]
实际应用中，高分辨率sar影像要求sar数据地面像元优于3米，可以为cosmo
‑
skymed雷达卫星星座、terrasar
‑
x/tandem
‑
x雷达卫星、radarsat
‑
2雷达卫星等高分辨率雷达卫星系统的聚束模式，高分辨率成像的条带模式，为获取桥梁重复轨道的高分辨率星载sar数据集提供数据源；同时，选取的sar影像之间的垂直基线需要尽量短，以terrasar
‑
x为例，sar数据集垂直基线分布要求在300米以内，以尽量降低dsm高程误差带来的影响。
[0085]
步骤2，在各幅高分辨率sar影像中，确定一幅高分辨率sar影像作为sar主影像，其他g
‑
1幅高分辨率sar影像作为sar从影像；
[0086]
步骤3，根据各幅高分辨率sar影像的成像轨道信息和影像后向散射信息，对各幅sar从影像进行旋转缩放采样操作，实现各幅sar从影像与sar主影像的影像配准，由此得到g
‑
1幅配准后的sar从影像；
[0087]
具体的，可以基于sar影像成像轨道信息求解主从sar影像的配准偏移量，对主从sar影像进行粗配准；再使用强度归一化互相关算法对主从sar影像进行精确配准，配准精度需要优于1/8个像素。
[0088]
步骤4，对于sar主影像和g
‑
1幅配准后的sar从影像，任意两幅sar影像形成一个干涉像对，由此生成所有干涉像对；
[0089]
对于每个干涉像对，其包含的两幅sar影像共轭相乘，生成一幅干涉相位图；
[0090]
因此，对于所有干涉像对，生成多幅干涉相位图，由此组成干涉相位图集合；
[0091]
具体的，忽略时间基线和空间基线，对所有高分辨率sar影像进行两两干涉组合，生成尽可能多的干涉像对，对于每个干涉像对，其包含的两幅sar影像共轭相乘，由此组成干涉相位图集合；参照附图2，为sar影像的干涉像对组合图。
[0092]
步骤5，对于干涉相位图集合中的每幅干涉相位图，进行预处理、掩膜处理和相位解缠处理，得到解缠差分干涉相位图；
[0093]
由此将干涉相位图集合，转化为解缠差分干涉相位图集合；
[0094]
步骤5具体为：
[0095]
步骤5.1，预处理：
[0096]
对每幅干涉相位图进行预处理，包括平地相位去除、地形相位去除和滤波处理，得到处理后的差分干涉相位图；
[0097]
本步骤中，采用以下方法，进行地形相位去除：
[0098]
获取覆盖桥梁监测对象的高精度dsm数据；基于高精度dsm数据，进行地形相位去除。
[0099]
其中：高精度dsm数据要求现势性强，其提供的桥梁监测对象的高程信息在桥梁建成以后；高精度dsm数据的精度和空间分辨率要求高于aw3d30 dsm数据。实际应用中，高分七号立体像对、worldview
‑
3立体像对、terrasar
‑
x/tandem
‑
x双星系统、无人机倾斜摄影测量等可为研究区提供高分辨率高精度的dsm数据。
[0100]
因此，利用sar影像轨道信息和高精度dsm数据去除平地相位和地形相位，并采用goldstein等滤波方法对干涉相位图进行滤波，削弱噪声的影响，参照附图3，为处理后的差分干涉相位图的示意图。
[0101]
步骤5.2，掩膜处理：
[0102]
步骤5.2.1，制作桥梁监测对象掩膜；其中，所述桥梁监测对象掩膜是指：在掩膜图像中，桥梁监测对象的像素点为1，其他所有位置的像素点为0；
[0103]
步骤5.2.2，采用桥梁监测对象掩膜，对每幅处理后的差分干涉相位图进行掩膜处理，得到掩膜处理后的差分干涉相位图；其中，在掩膜处理后的差分干涉相位图中，桥梁监测对象的像素点相位保持不变，其他位置的像素点相位统一变为0；
[0104]
具体的，通过制作桥梁监测对象掩膜，对每幅处理后的差分干涉相位图进行掩膜处理，使差分干涉相位图仅保留桥梁监测对象的影像，降低桥梁外其他背景影像的影响。
[0105]
步骤5.3，相位解缠处理：
[0106]
对于每幅掩膜处理后的差分干涉相位图，均采用多种相位解缠算法进行相位解缠运算，由此得到多幅解缠差分干涉相位图，并采用筛选算法，剔除存在解缠误差的差分干涉相位图，保留相位解缠结果正确的解缠差分干涉相位图；
[0107]
将多种相位解缠算法获取的解缠结果正确的解缠差分干涉相位图合并，最终得到多个解缠差分干涉相位图。
[0108]
具体的，对于每幅掩膜处理后的差分干涉相位图，在桥梁相位稳定区域选取参考点，分别采用最小费用流法、枝切法等多种相位解缠算法，对所有差分干涉相位图进行相位解缠；
[0109]
然后，通过人工目视解译或剖面线分析等方法检查解缠差分干涉相位图是否存在解缠误差，比如明显的相位跳跃，剔除存在解缠误差的差分干涉相位图，保留相位解缠结果正确的解缠差分干涉相位图；
[0110]
然后，将多种相位解缠算法获取的解缠结果正确的解缠差分干涉相位图合并，相同差分干涉组合保留一幅，最大化相位解缠结果正确的解缠差分干涉相位图，参照附图4，为解缠结果正确的解缠差分干涉相位图示意图。
[0111]
步骤6，对于步骤5得到的解缠差分干涉相位图集合，假设共有m幅解缠差分干涉相位图，分别表示为：解缠差分干涉相位图tu1，解缠差分干涉相位图tu2，...，解缠差分干涉相位图tu
m
；
[0112]
其中：对于任意解缠差分干涉相位图tu
j
，其中，j＝1，2，...，m，对应的干涉像对表示为：sar影像i
j0
和sar影像i
j1
，即：sar影像i
j0
和sar影像i
j1
共轭相乘后，经过预处理、掩膜处理和相位解缠处理后，得到解缠差分干涉相位图tu
j
；
[0113]
对于m幅解缠差分干涉相位图，共对应m个干涉像对，即：不考虑干涉像对包含共同sar影像时，m个干涉像对包含2m幅sar影像；再对2m幅sar影像去重处理，假设最终得到n 1幅sar影像；
[0114]
将n 1幅sar影像，按成像时间排序，依次表示为：sar影像i0，sar影像i1，...，sar影像i
n
，对应的成像时间依次为：t0，t1，...，t
n
；即：sar影像i0的成像时间为t0，sar影像i1的成像时间为t1，...，sar影像i
n
的成像时间为t
n
；
[0115]
因此，将m幅解缠差分干涉相位图中的任意解缠差分干涉相位图，表示为解缠差分干涉相位图tu
j
(i
c
，i
d
)，其中，c＞d，c≠d，c＝1，2，...，n；d＝0，1，2，...，n
‑
1，含义为：解缠差分干涉相位图tu
j
(i
c
，i
d
)，是由sar影像i
c
和sar影像i
d
共轭相乘后，经过预处理、掩膜处理和相位解缠处理后生成，解缠差分干涉相位图tu
j
(i
c
，i
d
)中，任意一个像元的解缠相位差，
可直接通过解缠差分干涉相位图读取获得；
[0116]
步骤7，m幅解缠差分干涉相位图和n 1幅sar影像的图像大小完全相同，具有相同数量的像元；
[0117]
对于任意解缠差分干涉相位图中的像元(l，p)，其中，l，p分别为像元的行号和列号，执行步骤7.1
‑
步骤7.3，计算得到像元(l，p)对应的在地理坐标系下的insar雷达视线向形变量相对位移变化序列：
[0118]
步骤7.1，建立以下方程：
[0119][0120]
其中：
[0121]
分别表示sar影像i1相对于sar影像i0在像元(l，p)位置的解缠相位差，sar影像i2相对于sar影像i0在像元(l，p)位置的解缠相位差，...，sar影像in相对于sar影像i0在像元(l，p)位置的解缠相位差；为未知待求值；
[0122]
分别表示解缠差分干涉相位图tu1在像元(l，p)位置的解缠相位差，解缠差分干涉相位图tu2在像元(l，p)位置的解缠相位差，...，解缠差分干涉相位图tu
m
在像元(l，p)位置的解缠相位差；可通过对应的解缠差分干涉相位图直接读取到；
[0123]
为m*n的系数矩阵a，为已知值，通过以下方式确定：
[0124]
假设对于系数矩阵a中的任意第j行，其中，j＝1，2，...，m，则有以下关系式：
[0125][0126]
为解缠差分干涉相位图tu
j
(i
c
，i
d
)在像元(l，p)位置的解缠相位差；解缠差分干涉相位图tu
j
(i
c
，i
d
)对应sar影像i
c
和sar影像i
d
；
[0127]
如果d＝0，则：在中，令前面的系数a
jc
为1，其他n
‑
1个系数均为0，由此确定了系数a
j1
，...，a
jn
的值；
[0128]
如果d≠0，则：在中，令前面的系数a
jc
为1，令前面的系数a
jd
为
‑
1，其他n
‑
2个系数均为0，由此确定了系数a
j1
，...，a
jn
的值；
[0129]
步骤7.2，求解步骤7.1建立的方程，得到桥梁在像元(l，p)位置在图像坐标系下的insar雷达视线向形变量相对相位变化序列即：桥梁在像
元(l，p)位置，相对于t0时刻，分别在t1，...，t
n
时刻在图像坐标系下的insar雷达视线向形变量相对相位变化值；
[0130]
为方便对步骤6和步骤7进行理解，下面列举一个实施例：
[0131]
参考图6：
[0132]
假设经步骤3处理后，共有4幅sar影像，分别为sar影像i0，sar影像i1，sar影像i2，sar影像i3。
[0133]
在步骤4中，任意两幅sar影像形成一个干涉像对，由此生成6个干涉像对。对于每个干涉像对，生成一幅干涉相位图，因此，共得到6幅干涉相位图，
[0134]
在步骤5中，通过对每幅干涉相位图进行预处理、掩膜处理和相位解缠处理，剔除存在解缠误差的差分干涉相位图，保留相位解缠结果正确的解缠差分干涉相位图，因此，本步骤中，最后得到3幅解缠差分干涉相位图，分别表示为：
[0135]
解缠差分干涉相位图tu3，对应的干涉像对为：sar影像i0和sar影像i3；其任意像元的解缠相位差表示为为已知值；
[0136]
解缠差分干涉相位图tu4，对应的干涉像对为：sar影像i0和sar影像i2；其任意像元的解缠相位差表示为为已知值；
[0137]
解缠差分干涉相位图tu5，对应的干涉像对为：sar影像i1和sar影像i3；其任意像元的解缠相位差表示为为已知值；
[0138]
因此，共得到m＝3幅解缠差分干涉相位图，对应的sar影像分别为i0、i1、i2和i3，为4幅，即：n 1＝4，n＝3。
[0139]
假设4幅sar影像，按成像时间排序，即为：sar影像i0、sar影像i1、sar影像i2和sar影像i3。
[0140]
由此建立以下方程：
[0141][0142]
其中：
[0143]
为未知待求值；
[0144]
上述方程可转化为以下方程组：
[0145][0146][0147][0148]
其中：
[0149]
解缠差分干涉相位图tu3(i
c
，i
d
)＝tu3(i3，i0)，d＝0，因此，a
13
＝1，a
11
和a
12
为0；
[0150]
解缠差分干涉相位图tu4(i
c
，i
d
)＝tu3(i2，i1)，d≠0，因此，a
22
＝1，a
21
＝
‑
1，a
23
为0；
[0151]
解缠差分干涉相位图tu5(i
c
，i
d
)＝tu3(i3，i1)，d≠0，因此，a
33
＝1，a
31
＝
‑
1，a
32
为0。
[0152]
由此确定系数矩阵。
[0153]
然后，根据方程组中未知待求值与已知值之间的数量关系，选择最小二乘或奇异值分解方法求解。
[0154]
上述系数确定的原理可以理解为：
[0155]
假设有m幅解缠差分干涉相位图解缠正确，涉及到的sar影像n 1幅。解缠差分干涉相位图上像元位置可用方位向坐标l和距离向坐标p来描述，第j幅解缠差分干涉相位图在像元(l，p)的解缠相位可表示为：
[0156][0157]
其中：
[0158]
第j幅解缠差分干涉相位图，是由成像时间为t
b
的sar影像和成像时间为t
a
的sar影像，经过共轭、预处理、掩膜处理和相位解缠处理后生成。
[0159]
j∈[1，m]，d(t
b
，l，p)和d(t
a
，l，p)为像元(l，p)在t
b
和t
a
时间相对于sar数据集起始时间雷达视线向的形变量；
[0160]
λ为雷达波长，r为像元(l，p)与卫星之间的距离，θ为像元局部入射角，b
j
⊥
为干涉对的时间基线，δz为高程残差，和为像元(l，p)在t
b
和t
a
时间的大气延迟相位，δn
j
为噪声相位。
[0161]
当使用的dsm高程精度较高时，且sar影像干涉像对垂直基线较短时，由高程残差引起的高程残差相位可以忽略；同时，由于桥梁覆盖范围往往较小，可认为该范围内大气分布相对较均匀，因而解缠差分干涉相位图中的大气延迟相位也可以忽略；噪声相位可通过滤波进行削弱，则针对桥梁第j幅差分干涉相位图解缠后的相位可优化为：
[0162][0163]
因此，t1...t
n
时间点上各sar影像某像元的相位用表示，表示为各差分干涉相位图解缠后该像元的相位值，则有：
[0164][0165]
知可组合为：其中为未知量，为已知量，a为一个m
×
n的矩阵，每一行与相乘对应着一幅解缠差分干涉相位图。为m个差分干涉对，包含m个方程，包含n个未知的sar影像对应的原始相位值。
[0166]
本步骤具体为：采用最小二乘或奇异值分解方法，求解步骤7.1建立的方程。
[0167]
具体的，，当m≥n，且矩阵a的秩为n，通过最小二乘法求解；
[0168][0169]
其中：
[0170][0171]
上标t代表矩阵的转秩。
[0172]
当解缠差分干涉相位图集分为l组，并且，矩阵a的秩小于n，为n
‑
l 1，使用奇异值分解方法求解。
[0173]
步骤7.3，依据sar成像几何参数和波长信息，将桥梁在像元(l，p)位置在图像坐标系下的insar雷达视线向形变量相对相位变化序列转换为图像坐标系下insar雷达视线向形变量相对位移变化序列δd(t1)，δd(t2)，...，δd(t
n
)，并将图像坐标系下insar雷达视线向形变量相对位移变化序列地理编码，得到地理坐标系下的insar雷达视线向形变量相对位移变化序列；表示为：
[0174]
通过以上方式，假设任意解缠差分干涉相位图中，桥梁监测对象覆盖n0个像元，则t1时刻，各个像元在地理坐标系下的insar雷达视线向形变量相对位移变化值，绘制形成桥梁形变图gr1，如图5所示，为桥梁形变图。
[0175]
t2时刻，各个像元在地理坐标系下的insar雷达视线向形变量相对位移变化值，绘制形成桥梁形变图gr2；
[0176]
依此类推
[0177]
t
n
时刻，各个像元在地理坐标系下的insar雷达视线向形变量相对位移变化值，绘制形成桥梁形变图gr
n
；
[0178]
步骤8，在桥梁监测对象布置num个gnss测点；
[0179]
从初始化时刻开始，各个gnss测点按采样间隔，可持续测得测点所在位置的三维形变位移值；
[0180]
因此，获得t
i
时刻各个gnss测点测得的三维形变位移值，其中，i＝1，2，...，n，采用t
i
时刻各个gnss测点测得的三维形变位移值，校正桥梁形变图gr
i
，得到校正后的桥梁形变图gr’i
，表示t
i
时刻相对于初始化时刻，各个像元在地理坐标系下的insar雷达视线向形变量位移变化值。
[0181]
具体采用以下方法，对桥梁形变图gr
i
进行校正：
[0182]
步骤8.1，获得ti时刻各个gnss测点测得的三维形变位移值；其中，对于任意第z个gnss测点，z＝1，2，...，num，其在t
i
时刻的三维形变位移值包括：d
xi
、d
yi
、d
ui
，分别代表东西向、南北向和垂向上的位移值；
[0183]
步骤8.2，将第z个gnss测点测得的三维形变位移值，采用以下公式，转化为雷达视线向形变测量值：
[0184][0185]
其中：
[0186]
θ为步骤2中sar主影像对应的雷达波局部入射角；
[0187]
α为步骤2中sar主影像的卫星轨道方位角；
[0188]
在t
i
时刻，第z个gnss测点测得的雷达视线向形变测量值；
[0189]
步骤8.3，对于第z个gnss测点，查找桥梁形变图gr
i
，获得第z个gnss测点所在位置的像元的insar雷达视线向形变量相对位移变化值
[0190]
步骤8.4，采用下式，得到校正偏移量
[0191][0192]
步骤8.5，使桥梁形变图gr
i
整体发生形变得到校正后的桥梁形变图gr’i
。
[0193]
步骤8的构思为：
[0194]
提取gnss设备所在位置，即gnss测点所在位置的桥梁insar时序同名点的雷达视线向形变测量值，求解得到雷达视线向形变测量值对桥梁各个时间的insar桥梁形变图的校正量值，并对桥梁各时间的insar桥梁形变图进行校正，获取insar校正后的桥梁形变分布和形变历史。
[0195]
基于以上步骤，可以获取融合了gnss监测数据的桥梁形变分布和形变历史。
[0196]
本发明提供的融合gnss数据与insar技术的桥梁变形监测方法具有以下优点：
[0197]
本发明为一种融合gnss数据与insar技术的桥梁变形监测方法，利用桥梁在长时间段相干性仍保留较好的特点，采用高精度dsm和高分辨率sar影像获取桥梁时间序列形变分布，利用gnss同期形变监测结果对桥梁insar时序结果进行校正，克服了传统gnss测量、水准测量等桥梁形变监测手段稀疏布设导致的形变监测空间分辨率低、insar技术手段获取的是相对未知形变情况的参考点的变形量和线性形变模型在桥梁形变监测适用性较差的不足，极大改善了桥梁基础设施形变监测的空间分辨率，为桥梁安全运营提供了更加全面的数据支撑，在桥梁形变监测领域有着广阔的应用前景。
[0198]
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。