一种基于永久散射体干涉的桥梁形变监测方法

1.本发明属于雷达信号处理技术领域，具体涉及一种基于永久散射体干涉的桥梁形变监测方法。


背景技术：

2.合成孔径雷达干涉技术(insar)是一种有效的监测地表形变的手段。相较于利用传感器现场监测的手段，它监测成本低，监测范围大；相较于光学卫星图像监测，它信息更加全面，能够利用相位信息监测毫米尺度的微小形变，且不会被云雨天气影响，具备全天候监测能力。
3.应用干涉技术进行桥梁的长期沉降监测时主要面临如下挑战。在长时间观测的情况下，合成孔径雷达(sar)图像将出现较严重的去相干现象；此外，桥梁的形变形式复杂，存在长期沉降形变、热膨胀形变等多种形变形式。这两点增大了桥梁形变监测的难度。为了提高监测的精度，目前主流的方法是利用永久散射体干涉技术(ps-insar)，它首先在sar图像序列中检测后向散射稳定的像素，称作永久散射体(ps)，随后将ps点连接组网，在此网络内求解参数。
4.关于形变参数求解，现有方法大致分为两类：第一类是相位解缠类方法，在去除地形相位之后，首先在网络内进行相位解缠，再根据解缠相位得到形变值。相位解缠类方法的优点在于不需要对研究区域的形变形式进行假设，因此能够较好地适应非线性形变的环境，但相位解缠法要求满足相位连续性假设这一条件，在ps-insar处理过程中较难满足，因此应用相对受限；第二类是时域模糊度估计方法，该方法首先对去地形相位进行建模，随后在预先设定的解空间中搜索使得时域相干系数最大的参数集合，将该参数值作为最终的形变参数求解结果。该方法回避了相位解缠步骤，且实现较为容易，因此是目前工程实现的主流方法，但它需要预先对形变进行建模，因此对复杂的非线性形变的监测能力较弱，同时解空间需人为设定也降低了算法的鲁棒性。


技术实现要素：

5.有鉴于此，本发明提供了一种基于永久散射体干涉的桥梁形变监测方法，能够更加稳健地监测桁架结构的桥梁形变信息。
6.一种基于永久散射体干涉的桥梁形变监测方法，包括如下步骤：
7.步骤一、选择待监测桥梁的sar图像，选取各sar图像的ps点；
8.步骤二、求解初步的ps点高程参数，具体包括：
9.s201、对去平地相位进行建模
10.将去平地相位后的ps点相位建模为：
[0011][0012]
其中，为去平地后相位，λ为雷达波长，b
⊥
为垂直基线距，r为斜距，θ为雷达下视
角，h为ps点高度相对平地的高度差，w为残余相位；
[0013]
s202、将ps点连接组网并解算参数：
[0014]
使用delaunay三角网络将ps点连接构网，并计算连接弧上的差分相位设定每一个弧连接的两ps点的高度差δh的搜索范围，将搜索范围中每一个待搜值，记作δhm，有：
[0015][0016]
其中γm为时域相干系数，n为sar图像数目，i为选择的sar图像的序号，表示第i幅sar图像中ps点的差分相位，表示第i幅sar图像中ps点的垂直基线距，j为虚数单位；使得γm最大的δhm即为搜索得到的该弧的差分高程估计结果；
[0017]
s203、最小二乘平差、粗差剔除，具体为：
[0018]
设所有ps点相位组成的相位向量为弧上的差分相位向量为k、m分别为ps点数量及连接弧的数量；则弧的连接关系用矩阵b来表示，即：
[0019][0020]
由此，根据搜索得到的ps点间的差分高程估计结果δhm，构成差分高程向量δh＝[δh0,δh1,
…
,δhm]
t
，通过加权最小二乘法得到各ps点高程的估计值，在得到ps点的高程估计值后，再计算ps点间的高程估计值的残差，并将残差大于给定阈值的弧删除；由此更新ps点间的弧的连接关系b，继而更新差分高程向量δh，然后再基于更新后的差分高程向量δh重新各ps点高程的估计值及残差，并判断残差是否大于给定阈值，将大于给定阈值的弧删除；以此类推，直至所有残差均小于或等于给定阈值；
[0021]
步骤三、求解形变参数基于步骤二最后更新的弧连接关系，得到各个ps点的高程值，
[0022]
即为公式(1)中的h，根据该公式计算ps点相位；
[0023]
s301、对剩余相位进行建模
[0024]
首先将基于公式(1)求得的高程相位去除，记剩余相位为对其进行建模：
[0025][0026]
其中，v为长期沉降速率，t为主副图像的时间差，k为热膨胀系数，δt为主副图像的温度差，h
′
为ps点真实高程与步骤二求得的h之差；
[0027]
s302、将ps点连接为高冗余度网络并解算参数
[0028]
首先使用自由连接网络连接ps点，即将所有满足给定距离阈值的ps点对连接起来；随后计算各连接弧的差分相位给定差分高程搜索范围δh
′m、差分形变速率范围δvm以及差分热膨胀系数范围δkm，将该三个参数在各自搜索范围的待搜索值进行组合，计算得到每一种组合参数下的时域相干系数：
[0029]
其中，表示本步骤计算的第i幅sar图像中ps点的差分相位；
[0030]
得到每条连接弧对应的所有参数组合下的时域相干系数γ后，按下述方法建立网络：
[0031]
1)、若时域相干系数γ的最大值γ
max
小于设定的阈值，则直接将该弧删除；
[0032]
2)、对经过1)删除后余下的每条弧，设定一个系数α，记录所有满足时域相干系数γ》αγ
max
的坐标(δh
′
,δv,δk)；
[0033]
3)、对记录的坐标进行dbscan聚类，如聚为一类则保留该弧，聚为多类则删除该弧；最终保留的弧组成的网络即为高冗余度网络；
[0034]
s303、最小二乘平差、粗差剔除
[0035]
与s203类似，基于s302最终得到的高冗余度网络的弧连接关系，根据时域相干系数最大值γ
max
确定搜索得到的差分高程向量δh
′
，差分形变速率向量δv以及差分热膨胀系数向量δk，通过加权最小二乘法分别得到各ps点的高程向量h
′
，形变速率向量v与热膨胀系数向量k；
[0036]
计算有弧连接的两ps点间的高程向量h
′
，形变速率向量v与热膨胀系数向量k残差，将残差大于对应设定阈值的弧删除，更新高冗余度网络；基于更新后高冗余度网络，再根据弧的差分高程向量δh
′
，差分形变速率向量δv以及差分热膨胀系数向量δk，通过加权最小二乘法分别得到各ps点的高程向量h
′
，形变速率向量v与热膨胀系数向量k；如此反复，将高程向量h
′
，形变速率向量v与热膨胀系数向量k进行粗差剔除，最终得到三个参数，完成桥梁形变监测。
[0037]
较佳的，所述步骤一中，选取ps点的方法为：
[0038]
根据轨道基线信息计算各图像之间的理论相干性，并将平均相干性最大的图像作为主图像；而后基于图像复相干系数的滑窗配准方法将所有副图像依次向主图像进行配准；配准完成后计算每幅干涉图的平地相位并去除，随后，计算图像中每个像素的振幅离差值，将小于给定阈值的像素点识别为ps点。
[0039]
本发明具有如下有益效果：
[0040]
本发明首先使用传统的ps-insar方法并只对高度相位进行建模以求解ps点高度信息，然后在第二步中使用改进的结合阈值筛选的自由组合网络，利用第一步初步得到的高度信息求得更精确的高度信息、长期沉降形变信息及周期性热膨胀形变信息。本发明与传统ps-insar方法比较如下：
[0041]
桁架结构桥梁识别出的ps点主要分布于桁架上下两端，因此高度差较大，因此待搜索参数的解空间很大，使用传统ps-insar方法难以得到稳健且正确的参数结果；本发明提出的分步迭代的方式能够优先搜索到相对正确的高程值，以此缩小解空间范围，得到更加稳健的参数结果。
[0042]
由桥梁形状的特点，识别出的ps点在图像中呈长条状分布，传统ps-insar中使用的delaunay三角网络冗余度较低，在后续处理中容易被切断为孤岛，无法得到正确的结果。同时桥梁形变形式较为复杂，长期沉降形变与周期性热膨胀形变同时存在，冗余度较低的delaunay三角网络难以提供足够的相位信息；本发明使用结合了基于聚类的阈值筛选步骤的自由组合网络，在提升网络冗余度，更有效地利用相位信息的同时对错误弧进行剔除，提升了估计结果的准确性。
[0043]
本发明通过分步迭代并结合高冗余度网络的ps-insar方法对桁架桥梁的形变进
行监测，能够获得更加准确且鲁棒的高度及形变参数估计结果，因此本发明的形变监测能力优于传统ps-insar方法。
附图说明
[0044]
图1为本发明的流程图；
[0045]
图2(a)为本发明实施例中桥梁的位置及光学、雷达图像情况，图2(b)为本发明实施例中桥梁ps点选择结果；
[0046]
图3(a)和图3(b)为应保留的弧其相干系数分布图及聚类结果，图3(c)和图3(d)为应删除的弧其相干系数分布图及聚类结果；
[0047]
图4(a)为传统方法在如表2所示情形(1)-(4)下的高程反演结果；图4(b)为本发明在如表2所示情形(1)-(4)下的高程反演结果
[0048]
图5(a)、图5(b)分别为本发明反演得到的长期沉降形变速率、热膨胀系数；图5(c)和图(d)分别为传统ps-insar方法反演得到的长期沉降形变速率、热膨胀系数；
[0049]
图6(a)、图6(b)和图6(c)分别为绘制在谷歌地球上的本发明反演得到的高程结果、长期沉降速率与热膨胀系数。
具体实施方式
[0050]
下面结合附图并举实施实例，对本发明进行详细描述。
[0051]
本发明公开了基于ps-insar的桁架桥梁形变监测方法，在求解高度及形变参数时，首先仅对高程进行建模，使用传统ps-insar方法求解高程初值。然后在第二步中建立高冗余度ps点连接网络，对长期沉降形变和周期性热膨胀形变进行建模，并在第一步求得的高程信息基础上迭代求得所有参数。
[0052]
如图1所示，为本发明的流程示意图，该方法具体包括以下步骤：
[0053]
步骤一、数据预处理及ps点选取：
[0054]
根据轨道基线信息计算各图像之间的理论相干性，并将平均相干性最大的图像作为主图像。而后基于图像复相干系数的滑窗配准方法将所有副图像依次向主图像进行配准。配准完成后计算每幅干涉图的平地相位并去除，随后，计算图像中每个像素的振幅离差值，将小于给定阈值的像素点识别为ps点。
[0055]
步骤二、求解初步的ps点高程参数：
[0056]
本步骤具体为对识别出的ps点实施子步骤s201～s203。
[0057]
s201、对去平地相位进行建模
[0058]
将去平地相位后的ps点相位建模为：
[0059][0060]
其中，为去平地后相位，λ为雷达波长，b
⊥
为垂直基线距，r为斜距，θ为雷达下视角，h为ps点高度相对平地的高度差，w为残余相位。
[0061]
s202、将ps点连接组网并解算参数：
[0062]
使用delaunay三角网络将ps点连接构网，并计算连接弧上的差分相位设定每一个弧连接的两ps点的高度差δh的搜索范围，将搜索范围中每一个待搜值，记作δhm，有：
[0063][0064]
其中γm为时域相干系数，n为sar图像数目，i为选择的sar图像的序号，表示第i幅sar图像中ps点的差分相位，表示第i幅sar图像中ps点的垂直基线距，j为虚数单位。使得γm最大的δhm即为搜索得到的该弧的差分高程估计结果。
[0065]
s203、最小二乘平差、粗差剔除
[0066]
设所有ps点相位组成的相位向量为弧上的差分相位向量为k、m分别为ps点数量及连接弧的数量。则弧的连接关系可以用矩阵b来表示，即：
[0067][0068]
由此，根据搜索得到的ps点间的差分高程估计结果δhm，构成差分高程向量δh＝[δh0,δh1,
…
,δhm]
t
，通过加权最小二乘法(wls)可以得到识别出的各ps点高程的估计值，在得到ps点的高程估计值后，再计算ps点间的高程估计值的残差，并将残差大于给定阈值的弧删除；由此更新ps点间的弧的连接关系b，继而更新差分高程向量δh，然后再基于更新后的差分高程向量δh重新各ps点高程的估计值及残差，并判断残差是否大于给定阈值，将大于给定阈值的弧删除；以此类推，重复该过程直至所有残差均小于或等于给定阈值。
[0069]
步骤三、求解形变参数基于步骤二最后更新的弧连接关系，得到各个ps点的高程值，
[0070]
即为公式(1)中的h，根据该公式计算ps点相位。
[0071]
s301、对剩余相位进行建模
[0072]
首先将基于公式(1)求得的高程相位去除，记剩余相位为对其进行建模：
[0073][0074]
其中，v为长期沉降速率，t为主副图像的时间差，k为热膨胀系数，δt为主副图像的温度差，h
′
为ps点真实高程与步骤二求得的高程之差，其余变量与步骤s201中意义一致。
[0075]
s302、将ps点连接为高冗余度网络并解算参数
[0076]
首先使用自由连接网络连接ps点，即将所有满足给定距离阈值的ps点对连接起来。随后与步骤s202中相同，计算各连接弧的差分相位给定差分高程搜索范围、差分形变速率范围以及差分热膨胀系数范围，将该三个参数在各自搜索范围的待搜索值进行组合，计算得到每一种组合参数下的时域相干系数。
[0077]
此时的时域相干系数计算式如下：
[0078]
得到每条连接弧对应的所有参数组合下的时域相干系数γ后，按下述方法建立网络：
[0079]
1)、若时域相干系数γ的最大值γ
max
小于设定的阈值，则直接将该弧删除；
[0080]
2)、对经过1)删除后余下的每条弧，设定一个系数α，记录所有满足时域相干系数γ》αγ
max
的坐标(δh
′
,δv,δk)；α一般取值为0.85。
[0081]
3)、对记录的坐标进行dbscan聚类，如聚为一类则保留该弧，聚为多类则删除该弧。
[0082]
最终保留的弧组成的网络即为高冗余度网络。
[0083]
由于已经求得各弧的γ
max
，可直接采用以下方法进行平差。
[0084]
s303、最小二乘平差、粗差剔除
[0085]
与s203类似，基于s302最终得到的高冗余度网络的弧连接关系，根据时域相干系数最大值γ
max
确定搜索得到的差分高程向量δh
′
，差分形变速率向量δv以及差分热膨胀系数向量δk，通过加权最小二乘法(wls)可以分别得到各ps点的高程向量h
′
，形变速率向量v与热膨胀系数向量k。
[0086]
计算有弧连接的两ps点间的高程向量h
′
，形变速率向量v与热膨胀系数向量k残差，将残差大于对应设定阈值的弧删除，更新高冗余度网络；基于更新后高冗余度网络，再根据弧的差分高程向量δh
′
，差分形变速率向量δv以及差分热膨胀系数向量δk，通过加权最小二乘法(wls)分别得到各ps点的高程向量h
′
，形变速率向量v与热膨胀系数向量k；如此反复，将高程向量h
′
，形变速率向量v与热膨胀系数向量k进行粗差剔除，最终得到三个参数，完成桥梁形变监测。
[0087]
实施实例：
[0088]
利用本发明对星载sar数据进行处理，以进一步验证所提技术的可行性与有效性。
[0089]
如图2所示，给出了待监测桁架桥梁实例的位置、光学图像、sar图像及sar图像中ps点选择结果，可以看到选择出的ps点均匀分布在桥梁上。
[0090]
表1给出了4种测试情形的详细信息，如图4所示，给出了传统ps-insar方法和本发明所提方法在如表1所示的4种测试情形下的高程估计结果图，可以看出在保留不同ps点数和搜索范围不同的情况下，传统ps-insar方法的估计结果不稳健，易受人工设定的超参数干扰；而本发明在第一步迭代中只建模高程信息的情况下，能够得到更加准确且稳健的估计结果。
[0091]
表1：测试用例信息
[0092][0093]
图3给出了高冗余度网络构建过程中使用聚类方法进行阈值筛选的示意图。图5给出了传统ps-insar方法和本发明所提方法反演出的长期沉降速率及热膨胀系数。可以看出，传统方法构造的低冗余度网络一旦被切断，反演得到的结果会出现跳变，造成错误。表2给出了本发明所提方法和传统ps-insar方法在反演出相关参数后的残余相位标准差，标准差小说明反演出的参数更加准确。两者共同说明了本发明提出的两步参数解算结合高冗余度网络在桁架桥梁形变监测上的优势。
[0094]
表2：残余相位评估结果
[0095][0096]
图5给出了本发明所提方法对待监测桁架桥梁的形变监测结果，结果表明钢结构的桥梁随温度季节性变化的热膨胀形变占据主要部分，长期形变速率与桥梁结构高度相关且处于正常范围内。
[0097]
综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

技术特征：
1.一种基于永久散射体干涉的桥梁形变监测方法，其特征在于，包括如下步骤：步骤一、选择待监测桥梁的sar图像，选取各sar图像的ps点；步骤二、求解初步的ps点高程参数，具体包括：s201、对去平地相位进行建模将去平地相位后的ps点相位建模为：其中，为去平地后相位，λ为雷达波长，b
⊥
为垂直基线距，r为斜距，θ为雷达下视角，h为ps点高度相对平地的高度差，w为残余相位；s202、将ps点连接组网并解算参数：使用delaunay三角网络将ps点连接构网，并计算连接弧上的差分相位设定每一个弧连接的两ps点的高度差δh的搜索范围，将搜索范围中每一个待搜值，记作δh
m
，有：其中γ
m
为时域相干系数，n为sar图像数目，i为选择的sar图像的序号，表示第i幅sar图像中ps点的差分相位，表示第i幅sar图像中ps点的垂直基线距，j为虚数单位；使得γ
m
最大的δh
m
即为搜索得到的该弧的差分高程估计结果；s203、最小二乘平差、粗差剔除，具体为：设所有ps点相位组成的相位向量为弧上的差分相位向量为k、m分别为ps点数量及连接弧的数量；则弧的连接关系用矩阵b来表示，即：由此，根据搜索得到的ps点间的差分高程估计结果δh
m
，构成差分高程向量δh＝[δh0,δh1,
…
,δh
m
]
t
，通过加权最小二乘法得到各ps点高程的估计值，在得到ps点的高程估计值后，再计算ps点间的高程估计值的残差，并将残差大于给定阈值的弧删除；由此更新ps点间的弧的连接关系b，继而更新差分高程向量δh，然后再基于更新后的差分高程向量δh重新各ps点高程的估计值及残差，并判断残差是否大于给定阈值，将大于给定阈值的弧删除；以此类推，直至所有残差均小于或等于给定阈值；步骤三、求解形变参数基于步骤二最后更新的弧连接关系，得到各个ps点的高程值，即为公式(1)中的h，根据该公式计算ps点相位；s301、对剩余相位进行建模首先将基于公式(1)求得的高程相位去除，记剩余相位为对其进行建模：其中，v为长期沉降速率，t为主副图像的时间差，k为热膨胀系数，δt为主副图像的温度差，h
′
为ps点真实高程与步骤二求得的h之差；s302、将ps点连接为高冗余度网络并解算参数
首先使用自由连接网络连接ps点，即将所有满足给定距离阈值的ps点对连接起来；随后计算各连接弧的差分相位给定差分高程搜索范围δh
′
m
、差分形变速率范围δv
m
以及差分热膨胀系数范围δk
m
，将该三个参数在各自搜索范围的待搜索值进行组合，计算得到每一种组合参数下的时域相干系数：其中，表示本步骤计算的第i幅sar图像中ps点的差分相位；得到每条连接弧对应的所有参数组合下的时域相干系数γ后，按下述方法建立网络：1)、若时域相干系数γ的最大值γ
max
小于设定的阈值，则直接将该弧删除；2)、对经过1)删除后余下的每条弧，设定一个系数α，记录所有满足时域相干系数γ>αγ
max
的坐标(δh
′
,δv,δk)；3)、对记录的坐标进行dbscan聚类，如聚为一类则保留该弧，聚为多类则删除该弧；最终保留的弧组成的网络即为高冗余度网络；s303、最小二乘平差、粗差剔除与s203类似，基于s302最终得到的高冗余度网络的弧连接关系，根据时域相干系数最大值γ
max
确定搜索得到的差分高程向量δh
′
，差分形变速率向量δv以及差分热膨胀系数向量δk，通过加权最小二乘法分别得到各ps点的高程向量h
′
，形变速率向量v与热膨胀系数向量k；计算有弧连接的两ps点间的高程向量h
′
，形变速率向量v与热膨胀系数向量k残差，将残差大于对应设定阈值的弧删除，更新高冗余度网络；基于更新后高冗余度网络，再根据弧的差分高程向量δh
′
，差分形变速率向量δv以及差分热膨胀系数向量δk，通过加权最小二乘法分别得到各ps点的高程向量h
′
，形变速率向量v与热膨胀系数向量k；如此反复，将高程向量h
′
，形变速率向量v与热膨胀系数向量k进行粗差剔除，最终得到三个参数，完成桥梁形变监测。2.如权利要求1所述的一种基于永久散射体干涉的桥梁形变监测方法，其特征在于，所述步骤一中，选取ps点的方法为：根据轨道基线信息计算各图像之间的理论相干性，并将平均相干性最大的图像作为主图像；而后基于图像复相干系数的滑窗配准方法将所有副图像依次向主图像进行配准；配准完成后计算每幅干涉图的平地相位并去除，随后，计算图像中每个像素的振幅离差值，将小于给定阈值的像素点识别为ps点。

技术总结
本发明公开了一种基于永久散射体干涉的桥梁形变监测方法，首先使用传统的PS-InSAR方法并只对高度相位进行建模以求解PS点高度信息，然后用改进的结合阈值筛选的自由组合网络，利用第一步初步得到的高度信息求得更精确的高度信息、长期沉降形变信息及周期性热膨胀形变信息；本发明通过分步迭代并结合高冗余度网络的PS-InSAR方法对桁架桥梁的形变进行监测，能够获得更加准确且鲁棒的高度及形变参数估计结果，因此本发明的形变监测能力优于传统PS-InSAR方法。InSAR方法。InSAR方法。


技术研发人员：丁泽刚 李根 贾笑天 李莫凡 胡子浛
受保护的技术使用者：北京理工大学
技术研发日：2022.11.04
技术公布日：2023/9/25