CN114660602A - 一种广域InSAR形变速率自适应拼接融合方法、装置、设备及介质 - Google Patents
一种广域InSAR形变速率自适应拼接融合方法、装置、设备及介质 Download PDFInfo
- Publication number
- CN114660602A CN114660602A CN202210297371.5A CN202210297371A CN114660602A CN 114660602 A CN114660602 A CN 114660602A CN 202210297371 A CN202210297371 A CN 202210297371A CN 114660602 A CN114660602 A CN 114660602A
- Authority
- CN
- China
- Prior art keywords
- deformation
- deformation rate
- sar image
- sar
- east
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000007500 overflow downdraw method Methods 0.000 title claims description 9
- 230000003044 adaptive effect Effects 0.000 title description 12
- 238000012544 monitoring process Methods 0.000 claims abstract description 75
- 238000000034 method Methods 0.000 claims abstract description 59
- 230000000007 visual effect Effects 0.000 claims abstract description 24
- 238000003384 imaging method Methods 0.000 claims abstract description 15
- 230000001174 ascending effect Effects 0.000 claims description 16
- 239000011159 matrix material Substances 0.000 claims description 12
- 238000004590 computer program Methods 0.000 claims description 11
- 238000012937 correction Methods 0.000 claims description 11
- 238000004422 calculation algorithm Methods 0.000 claims description 9
- 238000006073 displacement reaction Methods 0.000 claims description 9
- 230000008569 process Effects 0.000 claims description 8
- NQLVQOSNDJXLKG-UHFFFAOYSA-N prosulfocarb Chemical compound CCCN(CCC)C(=O)SCC1=CC=CC=C1 NQLVQOSNDJXLKG-UHFFFAOYSA-N 0.000 claims description 6
- 241000690776 Hassar Species 0.000 claims description 5
- 238000012952 Resampling Methods 0.000 claims description 4
- 238000005457 optimization Methods 0.000 claims description 4
- 102100029469 WD repeat and HMG-box DNA-binding protein 1 Human genes 0.000 claims description 3
- 101710097421 WD repeat and HMG-box DNA-binding protein 1 Proteins 0.000 claims description 3
- 230000015572 biosynthetic process Effects 0.000 claims description 3
- 238000005755 formation reaction Methods 0.000 claims description 3
- 238000007781 pre-processing Methods 0.000 claims description 3
- 230000004927 fusion Effects 0.000 abstract description 9
- 238000005516 engineering process Methods 0.000 description 6
- 238000010586 diagram Methods 0.000 description 3
- 238000012545 processing Methods 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000001427 coherent effect Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000013500 data storage Methods 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 230000006870 function Effects 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 230000000630 rising effect Effects 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/88—Radar or analogous systems specially adapted for specific applications
- G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
- G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
- G01S13/9094—Theoretical aspects
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Radar, Positioning & Navigation (AREA)
- Physics & Mathematics (AREA)
- Electromagnetism (AREA)
- Computer Networks & Wireless Communication (AREA)
- General Physics & Mathematics (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明公开了一种广域InSAR形变速率自适应拼接融合方法、装置、设备及介质，方法为：获取广域目标区域内所有SAR图幅，解算各图幅沿雷达视线向的形变速率；利用各SAR图幅的成像信息，获得各图幅的空间覆盖范围和观测几何信息；基于各SAR图幅的覆盖范围、观测几何和形变速率，将各SAR图幅沿雷达视线向的形变速率转换到垂直向和东西向；利用相邻SAR图幅间影像重叠部分的垂直向和东西向的形变速率，改正各SAR图幅参考基准之间的差异；利用各SAR图幅基于统一参考基准的垂直向和东西向的形变速率，生成得到待研究广域目标区域的全域形变监测结果。本发明实现广域InSAR形变速率的自适应拼接融合，提升广域形变结果精度。
Description
技术领域
本发明属于InSAR领域，特别涉及一种广域InSAR形变速率自适应拼接融合方法、装置、设备及介质。
背景技术
对大区域范围内的地面变形信息进行高精度监测，可以获得形变场的整体信息，提高人们对地球物理现象的认识，具有重大科学意义。随着近年来合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar，InSAR)技术的飞速发展和越来越多在轨/预发射SAR卫星，InSAR技术以其面域、无接触、高精度、高空间分辨率以及不受云雨天气影响等优势已经成为地面形变监测的一种不可或缺的大地测量技术手段。而且，随着越来越丰富的InSAR数据，为大区域范围形变监测提供了可能。利用InSAR技术对广域范围内的地面形变进行监测时，往往需要利用多个图幅才能对待监测区域实现全覆盖。然而，由于InSAR侧视成像的观测几何特征和各图幅独立解算的形变监测策略，各图幅得到的InSAR形变结果为地表真实三维形变在各图幅雷达视线向的一维投影。而各图幅范围内每个监测点处的雷达入射角和飞行方位角等不相同，且各图幅形变解算的起始参考点之间也存在偏差，不绝对为零。这些因素阻碍了广域InSAR高精度形变产品的生成。
目前已有的大范围InSAR形变监测项目或研究的结果生成过程中通常不考虑上述问题，而直接将各图幅结果进行简单叠加或者做一些简单的常数偏移改正，未考虑各图幅结果之间观测几何的差异，且未利用平差思想对各图幅结果进行整体优化、统一参考基准，不能达到广域精密形变成果生成的精度要求，也没有一套完整的广域InSAR成果标准格式输出系统。目前，单图幅InSAR形变解算技术趋于成熟，越来越多SAR卫星在轨运行，以及高时空分辨率小图幅SAR卫星星座的建设，均为广域InSAR形变监测和多图幅形变结果的拼接融合提供应用需求和发展契机。同时，多图幅间参考基准的不统一和各相邻图幅间观测几何参数的不一致等问题阻碍着广域InSAR形变成果的高精度获取实际应用。
发明内容
本发明提供一种广域InSAR形变速率自适应拼接融合方法、装置、设备及介质，充分利用监测时间内覆盖大范围研究区域的多图幅InSAR数据集的成像和监测信息，实现广域InSAR形变速率的自适应拼接融合，可有效规避拼接融合误差，显著提升广域InSAR形变结果精度。
为实现上述技术目的，本发明采用如下技术方案：
一种广域InSAR形变速率自适应拼接融合方法，包括：
S1：获取监测时间内覆盖待研究的广域目标区域内所有SAR图幅的数据集，并利用时序InSAR算法解算得到各图幅数据集内沿雷达视线向的形变速率结果；
S2：利用各SAR图幅的成像信息，获得各图幅的空间覆盖范围和其内每个监测点的观测几何信息；
S3：基于各SAR图幅的覆盖范围、观测几何和形变速率信息，将各SAR图幅沿雷达视线向的形变速率转换到垂直向和东西向的形变速率；
S4：利用相邻SAR图幅间影像重叠部分的垂直向和东西向的形变速率，采用广域最小二乘平差方法改正各SAR图幅参考基准之间的差异，得到各SAR图幅基于统一参考基准的垂直向和东西向的形变速率；
S5：利用各SAR图幅的空间覆盖范围以及基于统一参考基准的垂直向和东西向的形变速率，生成得到待研究广域目标区域的全域形变监测结果。
进一步的，S1中所述雷达视线向的形变速率，为地表沿垂直、东西和南北三个真实位移方向的三维形变速率在雷达视线向的一维投影，表示为：
dLOS＝cosθ·dU+sinθ·sinα·dN-sinθ·cosα·dE
式中，dLOS为雷达视线向形变速率，dU,dN,dE分别为沿垂直、南北和东西三个真实位移方向的三维形变速率，θ,α分别表示SAR卫星雷达的入射角和飞行方位角。
进一步的，S3针对各SAR图幅的具体步骤包括：
S31：根据SAR图幅的空间覆盖范围和观测几何信息，将SAR图幅分类为升轨或降轨；
S32：针对广域目标区域内的地表监测点，分两种情况解算其在垂直、南北和东西三个方向的形变速率：
(1)若地表监测点同时存在升轨和降轨的SAR图幅覆盖，则忽略南北向形变分量对雷达视线向形变的贡献，依据升轨和降轨两种SAR图幅的形变速率和观测几何，解算得到地表监测点在垂直向和东西向的形变速率：
式中
分别为升轨和降轨SAR图幅对应的地表监测点形变速率，θasc,θdes分别为升轨和降轨SAR图幅对应的雷达入射角，αasc,αdes分别为升轨和降轨SAR图幅对应的雷达飞行方位角；dU,dN,dE分别为地表监测点在垂直、南北和东西三个真实位移方向的三维形变速率；
(2)若地表监测点只存在升轨或者降轨的SAR图幅覆盖，则忽略南北向和东西向形变分量对雷达视线向形变的贡献，利用以下投影几何关系将雷达视线向形变速率直接转换为垂直向形变速率：
式中，dLOS为SAR图幅对应的地表监测点形变速率，θ表示SAR卫星雷达的入射角。
进一步的，S4的具体步骤包括：
S41：根据SAR图幅的空间覆盖范围，得到各相邻SAR图幅之间的影像重叠区域范围；设步骤S1的数据集中共获取有M个SAR图幅，共有N个影像重叠区域；
S42：针对形成第i个影像重叠区域的两个相邻SAR图幅，计算每个SAR图幅在该影像重叠区域范围内所有对应的地表监测点在各个方向上的形变速率的平均值，记为该SAR图幅在该影像重叠区域各个方向上的形变速率均值，并记两幅SAR图幅在第i个影像重叠区域各个方向上的形变速率均值分别为defi 1,defi 2；
S43：以影像重叠区域i的两个SAR图幅在该影像重叠区域的形变均值之差
为观测值，以所有影像重叠区域的
之和最小为目标，采用广域最小二乘平差方法解算得到所有M个SAR图幅各方向上的形变速率改正量；其中
S44：使用M个SAR图幅各方向上的形变速率改正量，对S3得到的对应形变速率进行改正，得到各SAR图幅基于统一参考基准的垂直向和东西向的形变速率。
进一步的，步骤S43采用广域最小二乘平差方法对所有M个SAR图幅在各方向上的形变速率进行修正的表达式为：
V＝A·X-α
X＝(ATPA)-1·ATP·V
式中，X＝[X1,X2,…,XM]为M个SAR图幅形变速率的改正量Xi构成的矩阵，
为N个影像重叠区域的
构成的矩阵；A是大小为N×M的系数矩阵，每一行代表一个重叠区域，每一列代表一个SAR图幅，其中，在每一行，存在重叠区域的两个SAR图幅索引处对应的元素值分别为1和-1，其余元素值均为0；α为残余误差；P为各影像重叠区域观测值
的权重所构成的权重矩阵。
进一步的，S5生成得到待研究广域目标区域的全域形变监测结果，是以度为单位的标准图幅GeoTIFF和GRD格式文件，具体生成过程为：
S51：根据待监测的广域目标区域范围，设置其在经纬度地理坐标下外接矩形的最大和最小边界范围；
S52：利用探测到的边界范围，以度为单位设置待输出GeoTIFF文件的参数文件；其中，根据输出成果的具体要求设置参数文件的空间分辨率；
S53：基于各SAR图幅的空间覆盖范围和设定的GeoTIFF参数文件，通过探测各SAR图幅在每个输出文件范围内的覆盖情况，将其依次重采样到GeoTIFF文件中；
其中，对于每一个GeoTIFF文件，当待采样点处存在多个SAR图幅监测结果时，取多个结果的平均值为该点处的结果；
S54：基于生成的每个GeoTIFF文件的速率文件和参数文件，生成对应的GRD格式文件，并分别保存GeoTIFF和GRD格式两类文件至物理介质。
一种基于上述任一项所述广域InSAR形变速率自适应拼接融合方法方法的装置，包括：
数据预处理单元：获取监测时间内覆盖待研究的广域目标区域内所有SAR图幅的数据集，并利用时序InSAR算法解算得到各图幅数据集内沿雷达视线向的形变速率结果；
成像信息获取单元：利用各SAR图幅的成像信息，获得各图幅的空间覆盖范围和其内每个监测点的观测几何信息；
形变分量优化解算单元：基于各SAR图幅的覆盖范围、观测几何和形变速率信息，将各SAR图幅沿雷达视线向的形变速率转换到垂直向和东西向的形变速率；
统一参考基准单元：利用相邻SAR图幅间影像重叠部分的垂直向和东西向的形变速率，采用广域最小二乘平差方法改正各SAR图幅参考基准之间的差异，得到各SAR图幅基于统一参考基准的垂直向和东西向的形变速率；
成果输出单元：利用各SAR图幅的空间覆盖范围以及基于统一参考基准的垂直向和东西向的形变速率，生成得到待研究广域目标区域的全域形变监测结果；
一种电子设备，包括存储器及处理器，所述存储器中存储有计算机程序，所述计算机程序被所述处理器执行时，使得所述处理器实现上述任一项技术方案所述的广域InSAR形变速率自适应拼接融合方法。
一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现上述任一项技术方案所述的广域InSAR形变速率自适应拼接融合方法。
有益效果
本发明提供了一种广域InSAR形变速率自适应拼接融合方法、装置、设备及介质，通过收集监测时间内覆盖大范围研究区域的所有SAR图幅的数据集，利用时序InSAR算法解算得到每个图幅数据集内沿雷达视线向的形变速率结果；利用各SAR图幅的成像信息，获得各图幅的空间覆盖范围和其内每个监测点的观测几何信息；基于各图幅的形变速率和观测几何等信息，采用一种优化的垂直向/东西向形变解算方法将各图幅沿雷达视线向的形变速率转换到垂直向/南北向/东西向；利用各相邻图幅间影像重叠部分的垂直向/南北向/东西向的形变速率，采用一种广域最小二乘平差方法改正各图幅参考基准之间的差异；利用各图幅空间覆盖和改正后的垂直向/东西向形变速率，采用一种自适应影像重采样方法生成广域范围内以度为单位的标准图幅GeoTIFF和GRD格式文件。该方法计算简单，实现方便。
本发明的技术效果主要体现在以下几点：
第一、相对于目前利用最为广泛的单图幅InSAR形变监测技术，本发明结合多个图幅的InSAR形变监测结果进行拼接融合，可以获得大区域范围整体形变结果，提升了对广域形变的整体科学认识，突破多图幅InSAR结果高精度拼接融合的技术瓶颈；
第二、相对于目前InSAR大范围形变场获取所采取的多图幅结果的简单调整和叠加，本发明数据处理效率较高，可有效规避拼接融合误差，显著提升广域InSAR形变结果精度；
第三、随着目前可获取的SAR卫星数据越来越丰富，在一定监测范围内进行广域InSAR形变监测已经成为一种客观需求和趋势，本发明拓展了InSAR技术的应用前景，跨越了现有InSAR形变监测技术传统思维，具有重大科学和工程应用价值和意义。
附图说明
图1是本发明实施例所述方法的流程示意图；
图2是本发明实施例模拟的多图幅InSAR监测结果以及模拟的基准差异示意图；
图3是本发明实施例模拟和恢复出的多图幅形变和基准差异调整结果。
具体实施方式
下面对本发明的实施例作详细说明，本实施例以本发明的技术方案为依据开展，给出了详细的实施方式和具体的操作过程，对本发明的技术方案作进一步解释说明。
实施例1
如图1所示，本实施例提供一种广域InSAR形变速率自适应拼接融合方法，包括以下步骤：
S1：获取监测时间内覆盖待研究的广域目标区域内所有SAR图幅的数据集，并利用时序InSAR算法解算得到各图幅数据集内沿雷达视线向的形变速率结果。
所述时序InSAR算法为已有算法，如永久散射体InSAR(PS-InSAR)、短基线集InSAR(SBAS-InSAR)、相干点目标分析(IPTA)等时序InSAR数据处理技术；
所述雷达视线向的形变速率，为地表沿垂直、南北和东西三个真实位移方向的三维形变速率在雷达视线向的一维投影，表示为：
dLOS＝cosθ·dU+sinθ·sinα·dN-sinθ·cosα·dE
式中，dLOS为雷达视线向形变速率，dU,dN,dE分别为沿垂直、南北和东西三个真实位移方向的三维形变速率，θ,α分别表示SAR卫星雷达的入射角和飞行方位角。
由于SAR卫星近极轨飞行侧视成像的观测模式，其雷达视线向监测结果对南北向形变分量的敏感度很低，且地面垂直向形变分量为视线向监测结果的主要成分。
S2：利用各SAR图幅的成像信息，获得各图幅的空间覆盖范围和其内每个监测点的观测几何信息；观测几何信息包括每个监测点位处的雷达入射角和飞行方位角等。
S3：基于各SAR图幅的覆盖范围、观测几何和形变速率信息，将各SAR图幅沿雷达视线向的形变速率转换到垂直向和东西向的形变速率；
针对各SAR图幅的具体步骤包括：
S31：根据SAR图幅的空间覆盖范围和观测几何信息，将SAR图幅分类为升轨或降轨；
S32：针对广域目标区域内的地表监测点，分两种情况解算其在垂直、南北和东西三个方向的形变速率：
(1)若地表监测点同时存在升轨和降轨的SAR图幅覆盖，则忽略南北向形变分量对雷达视线向形变的贡献(由于SAR雷达卫星沿近极轨侧视成像的特点，其沿雷达视线向的监测结果对南北向形变分量不敏感，因此可忽略南北向分量，记dN＝0)，依据升轨和降轨两种SAR图幅的形变速率和观测几何，解算得到地表监测点在垂直向和东西向的形变速率：
式中
分别为升轨和降轨SAR图幅对应的地表监测点形变速率，θasc,θdes分别为升轨和降轨SAR图幅对应的雷达入射角，αasc,αdes分别为升轨和降轨SAR图幅对应的雷达飞行方位角；dU,dN,dE分别为地表监测点在垂直、南北和东西三个真实位移方向的三维形变速率；
(2)若地表监测点只存在升轨或者降轨的SAR图幅覆盖，则忽略南北向和东西向形变分量对雷达视线向形变的贡献(由于SAR雷达卫星沿近极轨侧视成像的特点，其沿雷达视线向的监测结果对南北向形变分量不敏感，且垂直向分量为主分量，因此可忽略南北向分量和东西向量，记dN＝0，dE＝0)，利用以下投影几何关系将雷达视线向形变速率直接转换为垂直向形变速率：
式中，dLOS为SAR图幅对应的地表监测点形变速率，θ表示SAR卫星雷达的入射角。
本实施例在监测点同时存在升轨和降轨多个观测几何的SAR图幅数据集时，将视线向形变速率的东西向分量考虑进来，可以得到更高精度的垂直向形变速率和东西向形变速率。
S4：利用相邻SAR图幅间影像重叠部分的垂直向和东西向的形变速率，采用广域最小二乘平差方法改正各SAR图幅参考基准之间的差异，得到各SAR图幅基于统一参考基准的垂直向和东西向的形变速率；S4的具体步骤包括：
S41：根据SAR图幅的空间覆盖范围，得到各相邻SAR图幅之间的影像重叠区域范围；设步骤S1的数据集中共获取有M个SAR图幅，共有N个影像重叠区域；
S42：针对形成第i个影像重叠区域的两个相邻SAR图幅，计算每个SAR图幅在该影像重叠区域范围内所有对应的地表监测点在各个方向上的形变速率的平均值，记为该SAR图幅在该影像重叠区域各个方向上的形变速率均值，并记两幅SAR图幅在第i个影像重叠区域各个方向上的形变速率均值分别为defi 1,defi 2；
由于本实施例在步骤S3的两种情况均未考虑南北向分量，因此本步骤中的形变速率仅包括垂直向的形变速率，或者垂直向和东西向两个维度的数据；
S43：以影像重叠区域i的两个SAR图幅在该影像重叠区域的形变均值之差
为观测值，以所有影像重叠区域的
之和最小为目标，采用广域最小二乘平差方法解算得到所有M个SAR图幅各方向上的形变速率改正量；其中
步骤S43采用广域最小二乘平差方法对所有M个SAR图幅在各方向上的形变速率进行修正的表达式为：
V＝A·X-α
X＝(ATPA)-1·ATP·V
式中，X＝[X1,X2,…,XM]为M个SAR图幅形变速率的改正量Xi构成的矩阵，
为N个影像重叠区域的
构成的矩阵；A是大小为N×M的系数矩阵，每一行代表一个重叠区域，每一列代表一个SAR图幅，其中，在每一行，存在重叠区域的两个SAR图幅索引处对应的元素值分别为1和-1，其余元素值均为0；α为残余误差；P为各影像重叠区域观测值
的权重所构成的权重矩阵，通常由各重叠区域面积和其内观测点质量和数量确定；
S44：使用M个SAR图幅各方向上的形变速率改正量，对S3得到的对应形变速率进行改正，得到各SAR图幅基于统一参考基准的垂直向和东西向的形变速率。
S5：利用各SAR图幅的空间覆盖范围以及基于统一参考基准的垂直向和东西向的形变速率，生成得到待研究广域目标区域的全域形变监测结果。
本实施例中，生成的全域形变监测结果，是以度为单位的标准图幅GeoTIFF和GRD格式文件，具体生成过程为：
S51：根据待监测的广域目标区域范围，设置其在经纬度地理坐标下外接矩形的边界范围(取整数)；
S52：利用探测到的边界范围，以度为单位设置待输出GeoTIFF文件的参数文件；其中，根据输出成果的具体要求设置参数文件的空间分辨率；
S53：基于各SAR图幅的空间覆盖范围和设定的GeoTIFF参数文件，通过探测各SAR图幅在每个输出文件范围内的覆盖情况，将其依次重采样到GeoTIFF文件中；
其中，对于每一个GeoTIFF文件，当待采样点处存在多个SAR图幅监测结果时，取多个结果的平均值为该点处的结果；
S54：基于生成的每个GeoTIFF文件的速率文件和参数文件，生成对应的GRD格式文件，并分别保存GeoTIFF和GRD格式两类文件至物理介质，如硬盘等。
由于本实施例所述方法涉及广域InSAR结果多图幅之间的重叠区域及其参数的自动识别，雷达视线向形变的优化分解，多图幅结果之间参考基准统一，以及广域成果生成等多个步骤，有些步骤利用模拟实验较难表达。因此，下述将结合部分图来进一步解释本发明提供的所述方法及其部分特性和优势。
以某地区的多图幅InSAR形变监测为例，假设图幅数量M＝6，各图幅的覆盖范围如图2(a)所示。通过在各图幅中随机加入均值为0，均方根误差为1cm的高斯随机噪声(如图2(b)所示)，可以模拟得到参考基准存在差异的多图幅形变速率结果(如图2(c)所示)。将这些参考基准存在差异的影像进行简单叠加，可以看到各图幅之间存在加大的整体差异(图3(a))。利用本发明所述广域最小二乘整体平差方法对各图幅范围重叠部分进行处理，构建平差改正模型，即可对各图幅结果基准之间的差异进行弥补(图3(c))，进而得到参考基准统一的广域InSAR形变结果(图3(b))。此外，我们展示了图2(a)和图3(a-b)中点划线所示剖面线处的形变速率真值、模拟值和平差解算值，如图3(d)所示。
从结果图中可以看出，由于参考基准之间存在差异，利用多个图幅进行广域InSAR形变测量时，各图幅之间参考基准的差异将对生成广域结果造成困扰。而本发明可以有效消除该因素的影响。计算出真值与结果值之间的均方根误差为0.011cm。从而验证了本发明所述方法的可行性和可靠性。
实施例2
本实施例提供一种基于实施例1所述方法的广域InSAR形变速率自适应拼接融合装置，包括：
数据预处理单元：获取监测时间内覆盖待研究的广域目标区域内所有SAR图幅的数据集，并利用时序InSAR算法解算得到各图幅数据集内沿雷达视线向的形变速率结果；
成像信息获取单元：利用各SAR图幅的成像信息，获得各图幅的空间覆盖范围和其内每个监测点的观测几何信息；
形变分量优化解算单元：基于各SAR图幅的覆盖范围、观测几何和形变速率信息，将各SAR图幅沿雷达视线向的形变速率转换到垂直向和东西向的形变速率；
统一参考基准单元：利用相邻SAR图幅间影像重叠部分的垂直向和东西向的形变速率，采用广域最小二乘平差方法改正各SAR图幅参考基准之间的差异，得到各SAR图幅基于统一参考基准的垂直向和东西向的形变速率；
成果输出单元：利用各SAR图幅的空间覆盖范围以及基于统一参考基准的垂直向和东西向的形变速率，生成得到待研究广域目标区域的全域形变监测结果；
数据存储单元：用于将上述生成的GeoTIFF和GRD格式的文件存储到固态/机械/移动硬盘等物理介质中。
上述各模块的具体实现过程请参照实施例1所述方法的对应过程。应当理解，上述单元模块的具体实现过程参照方法内容，本发明在此不进行具体的赘述，且上述功能模块单元的划分仅仅是一种逻辑功能的划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。同时，上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。
应当理解，本发明各个实施例中的功能单元模块可以集中在一个处理单元中，也可以是各个单元模块单独物理存在，也可以是两个或两个以上的单元模块集成在一个单元模块中，可以采用硬件或软件的形式来实现。
实施例3
本实施例提供一种电子设备，包括存储器及处理器，所述存储器中存储有计算机程序，所述计算机程序被所述处理器执行时，使得所述处理器实现实施例1所述的广域InSAR形变速率自适应拼接融合方法。
实施例4
本实施例提供一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现上述任一项技术方案所述的广域InSAR形变速率自适应拼接融合方法。
以上实施例为本申请的优选实施例，本领域的普通技术人员还可以在此基础上进行各种变换或改进，在不脱离本申请总的构思的前提下，这些变换或改进都应当属于本申请要求保护的范围之内。
Claims (9)
1.一种广域InSAR形变速率自适应拼接融合方法，其特征在于，包括：
S1：获取监测时间内覆盖待研究的广域目标区域内所有SAR图幅的数据集，并利用时序InSAR算法解算得到各图幅数据集内沿雷达视线向的形变速率结果；
S2：利用各SAR图幅的成像信息，获得各图幅的空间覆盖范围和其内每个监测点的观测几何信息；
S3：基于各SAR图幅的覆盖范围、观测几何和形变速率信息，将各SAR图幅沿雷达视线向的形变速率转换到垂直向和东西向的形变速率；
S4：利用相邻SAR图幅间影像重叠部分的垂直向和东西向的形变速率，采用广域最小二乘平差方法改正各SAR图幅参考基准之间的差异，得到各SAR图幅基于统一参考基准的垂直向和东西向的形变速率；
S5：利用各SAR图幅的空间覆盖范围以及基于统一参考基准的垂直向和东西向的形变速率，生成得到待研究广域目标区域的全域形变监测结果。
2.根据权利要求1所述的方法，其特征在于，S1中所述雷达视线向的形变速率，为地表沿垂直、东西和南北三个真实位移方向的三维形变速率在雷达视线向的一维投影，表示为：
dLOS＝cosθ·dU+sinθ·sinα·dN-sinθ·cosα·dE
式中，dLOS为雷达视线向形变速率，dU,dN,dE分别为沿垂直、南北和东西三个真实位移方向的三维形变速率，θ,α分别表示SAR卫星雷达的入射角和飞行方位角。
3.根据权利要求1所述的方法，其特征在于，S3针对各SAR图幅的具体步骤包括：
S31：根据SAR图幅的空间覆盖范围和观测几何信息，将SAR图幅分类为升轨或降轨；
S32：针对广域目标区域内的地表监测点，分两种情况解算其在垂直、南北和东西三个方向的形变速率：
(1)若地表监测点同时存在升轨和降轨的SAR图幅覆盖，则忽略南北向形变分量对雷达视线向形变的贡献，依据升轨和降轨两种SAR图幅的形变速率和观测几何，解算得到地表监测点在垂直向和东西向的形变速率：
式中
分别为升轨和降轨SAR图幅对应的地表监测点形变速率，θasc,θdes分别为升轨和降轨SAR图幅对应的雷达入射角，αasc,αdes分别为升轨和降轨SAR图幅对应的雷达飞行方位角；dU,dN,dE分别为地表监测点在垂直、南北和东西三个真实位移方向的三维形变速率；
(2)若地表监测点只存在升轨或者降轨的SAR图幅覆盖，则忽略南北向和东西向形变分量对雷达视线向形变的贡献，利用以下投影几何关系将雷达视线向形变速率直接转换为垂直向形变速率：
式中，dLOS为SAR图幅对应的地表监测点形变速率，θ表示SAR卫星雷达的入射角。
4.根据权利要求1所述的方法，其特征在于，S4的具体步骤包括：
S41：根据SAR图幅的空间覆盖范围，得到各相邻SAR图幅之间的影像重叠区域范围；设步骤S1的数据集中共获取有M个SAR图幅，共有N个影像重叠区域；
S42：针对形成第i个影像重叠区域的两个相邻SAR图幅，计算每个SAR图幅在该影像重叠区域范围内所有对应的地表监测点在各个方向上的形变速率的平均值，记为该SAR图幅在该影像重叠区域各个方向上的形变速率均值，并记两幅SAR图幅在第i个影像重叠区域各个方向上的形变速率均值分别为defi 1,defi 2；
S43：以影像重叠区域i的两个SAR图幅在该影像重叠区域的形变均值之差
为观测值，以所有影像重叠区域的
之和最小为目标，采用广域最小二乘平差方法解算得到所有M个SAR图幅各方向上的形变速率改正量；其中
S44：使用M个SAR图幅各方向上的形变速率改正量，对S3得到的对应形变速率进行改正，得到各SAR图幅基于统一参考基准的垂直向和东西向的形变速率。
5.根据权利要求4所述的方法，其特征在于，步骤S43采用广域最小二乘平差方法对所有M个SAR图幅在各方向上的形变速率进行修正的表达式为：
V＝A·X-α
X＝(ATPA)-1·ATP·V
式中，X＝[X1,X2,…,XM]为M个SAR图幅形变速率的改正量Xi构成的矩阵，
为N个影像重叠区域的
构成的矩阵；A是大小为N×M的系数矩阵，每一行代表一个重叠区域，每一列代表一个SAR图幅，其中，在每一行，存在重叠区域的两个SAR图幅索引处对应的元素值分别为1和-1，其余元素值均为0；α为残余误差；P为各影像重叠区域观测值
的权重所构成的权重矩阵。
6.根据权利要求1所述的方法，其特征在于，S5生成得到待研究广域目标区域的全域形变监测结果，是以度为单位的标准图幅GeoTIFF和GRD格式文件，具体生成过程为：
S51：根据待监测的广域目标区域范围，设置其在经纬度地理坐标下外接矩形的最大和最小边界范围；
S52：利用探测到的边界范围，以度为单位设置待输出GeoTIFF文件的参数文件；其中，根据输出成果的具体要求设置参数文件的空间分辨率；
S53：基于各SAR图幅的空间覆盖范围和设定的GeoTIFF参数文件，通过探测各SAR图幅在每个输出文件范围内的覆盖情况，将其依次重采样到GeoTIFF文件中；
其中，对于每一个GeoTIFF文件，当待采样点处存在多个SAR图幅监测结果时，取多个结果的平均值为该点处的结果；
S54：基于生成的每个GeoTIFF文件的速率文件和参数文件，生成对应的GRD格式文件，并分别保存GeoTIFF和GRD格式两类文件至物理介质。
7.一种基于权利要求1～6任一项所述方法的装置，其特征在于，包括：
数据预处理单元：获取监测时间内覆盖待研究的广域目标区域内所有SAR图幅的数据集，并利用时序InSAR算法解算得到各图幅数据集内沿雷达视线向的形变速率结果；
成像信息获取单元：利用各SAR图幅的成像信息，获得各图幅的空间覆盖范围和其内每个监测点的观测几何信息；
形变分量优化解算单元：基于各SAR图幅的覆盖范围、观测几何和形变速率信息，将各SAR图幅沿雷达视线向的形变速率转换到垂直向和东西向的形变速率；
统一参考基准单元：利用相邻SAR图幅间影像重叠部分的垂直向和东西向的形变速率，采用广域最小二乘平差方法改正各SAR图幅参考基准之间的差异，得到各SAR图幅基于统一参考基准的垂直向和东西向的形变速率；
成果输出单元：利用各SAR图幅的空间覆盖范围以及基于统一参考基准的垂直向和东西向的形变速率，生成得到待研究广域目标区域的全域形变监测结果；
8.一种电子设备，包括存储器及处理器，所述存储器中存储有计算机程序，其特征在于，所述计算机程序被所述处理器执行时，使得所述处理器实现如权利要求1～6中任一项所述的方法。
9.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现如权利要求1～6中任一项所述的方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210297371.5A CN114660602B (zh) | 2022-03-24 | 2022-03-24 | 一种广域InSAR形变速率自适应拼接融合方法、装置、设备及介质 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210297371.5A CN114660602B (zh) | 2022-03-24 | 2022-03-24 | 一种广域InSAR形变速率自适应拼接融合方法、装置、设备及介质 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114660602A true CN114660602A (zh) | 2022-06-24 |
| CN114660602B CN114660602B (zh) | 2024-05-24 |
Family
ID=82030821
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210297371.5A Active CN114660602B (zh) | 2022-03-24 | 2022-03-24 | 一种广域InSAR形变速率自适应拼接融合方法、装置、设备及介质 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114660602B (zh) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115469308A (zh) * | 2022-08-25 | 2022-12-13 | 中国地震局地质研究所 | 多轨道InSAR震间形变速率场拼接方法、装置、设备及介质 |
| CN115951354A (zh) * | 2023-02-14 | 2023-04-11 | 中国铁道科学研究院集团有限公司铁道建筑研究所 | 一种融合升降轨的D-InSAR形变监测方法 |
| CN116299245A (zh) * | 2023-05-11 | 2023-06-23 | 中山大学 | 一种时序InSAR形变速率结果自适应镶嵌校正方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109029344A (zh) * | 2018-07-10 | 2018-12-18 | 湖南中科星图信息技术有限公司 | 一种基于高分影像和升降轨InSAR的堤坝沉降监测方法 |
| CN109738892A (zh) * | 2019-01-24 | 2019-05-10 | 中南大学 | 一种矿区地表高时空分辨率三维形变估计方法 |
| CN111650579A (zh) * | 2020-06-12 | 2020-09-11 | 中南大学 | 一种岩移参数自适应获取的InSAR矿区三维形变估计方法、装置及介质 |
| WO2021227423A1 (zh) * | 2020-05-13 | 2021-11-18 | 深圳大学 | 基于动态基线的InSAR数字高程模型构建方法及系统 |
-
2022
- 2022-03-24 CN CN202210297371.5A patent/CN114660602B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109029344A (zh) * | 2018-07-10 | 2018-12-18 | 湖南中科星图信息技术有限公司 | 一种基于高分影像和升降轨InSAR的堤坝沉降监测方法 |
| CN109738892A (zh) * | 2019-01-24 | 2019-05-10 | 中南大学 | 一种矿区地表高时空分辨率三维形变估计方法 |
| WO2021227423A1 (zh) * | 2020-05-13 | 2021-11-18 | 深圳大学 | 基于动态基线的InSAR数字高程模型构建方法及系统 |
| CN111650579A (zh) * | 2020-06-12 | 2020-09-11 | 中南大学 | 一种岩移参数自适应获取的InSAR矿区三维形变估计方法、装置及介质 |
Non-Patent Citations (3)
| Title |
|---|
| LEI ZHANG; XIAOLI DING; ZHONG LU; HYUNG-SUP JUNG; JUN HU; GUANGCAI FENG: "A Novel Multitemporal InSAR Model for Joint Estimation of Deformation Rates and Orbital Errors", IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, vol. 52, no. 06, 6 November 2013 (2013-11-06) * |
| 敖萌;张路;廖明生;张丽;: "基于方差分量估计的多源InSAR数据自适应融合形变测量", 地球物理学报, no. 08, 3 August 2020 (2020-08-03) * |
| 胡俊;李志伟;朱建军;任小冲;丁晓利;: "融合升降轨SAR干涉相位和幅度信息揭示地表三维形变场的研究", 中国科学:地球科学, no. 03, 15 March 2010 (2010-03-15) * |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115469308A (zh) * | 2022-08-25 | 2022-12-13 | 中国地震局地质研究所 | 多轨道InSAR震间形变速率场拼接方法、装置、设备及介质 |
| CN115951354A (zh) * | 2023-02-14 | 2023-04-11 | 中国铁道科学研究院集团有限公司铁道建筑研究所 | 一种融合升降轨的D-InSAR形变监测方法 |
| CN116299245A (zh) * | 2023-05-11 | 2023-06-23 | 中山大学 | 一种时序InSAR形变速率结果自适应镶嵌校正方法 |
| CN116299245B (zh) * | 2023-05-11 | 2023-07-28 | 中山大学 | 一种时序InSAR形变速率结果自适应镶嵌校正方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114660602B (zh) | 2024-05-24 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN114660602A (zh) | 一种广域InSAR形变速率自适应拼接融合方法、装置、设备及介质 | |
| CN111708038B (zh) | 基于姿态传感器和gnss的无人船激光雷达点云数据校正方法 | |
| CN108983232B (zh) | 一种基于邻轨数据的InSAR二维地表形变监测方法 | |
| EP3132231B1 (en) | A method and system for estimating information related to a vehicle pitch and/or roll angle | |
| Xie et al. | Design and data processing of China's first spaceborne laser altimeter system for earth observation: GaoFen-7 | |
| CN109613583B (zh) | 基于单星与地面站测向及联合测时差的无源目标定位方法 | |
| CN106373159A (zh) | 一种简化的无人机多目标定位方法 | |
| JP2008506167A (ja) | 画像に関連するロケーションを確定する方法および装置 | |
| CN103697883B (zh) | 一种基于天际线成像的飞行器水平姿态确定方法 | |
| CN103129752A (zh) | 一种基于地面导航的光学遥感卫星姿态角误差动态补偿方法 | |
| CN102346033A (zh) | 基于卫星观测角误差估计的直接定位方法及系统 | |
| CN113409400A (zh) | 一种基于自动跟踪的机载光电系统目标地理定位方法 | |
| CN116401833A (zh) | 一种基于sgp4模型的卫星轨道计算方法 | |
| US7768631B1 (en) | Method and system for providing a known reference point for an airborne imaging platform | |
| Brunel et al. | Operational AVHRR navigation results | |
| CN110017833B (zh) | 基于像元类地模型的全屏像点地理坐标定位方法 | |
| Goncalves et al. | Accuracy analysis of DEMs derived from ASTER imagery | |
| CN113483739B (zh) | 海上目标位置测量方法 | |
| CN115077560A (zh) | 一种船载可见光及中波红外系统光轴平行度动态检测方法 | |
| CN111486850B (zh) | 一种对地观测卫星在轨帧频调整策略 | |
| Sheng | Theoretical analysis of the iterative photogrammetric method to determining ground coordinates from photo coordinates and a DEM | |
| CN103217145B (zh) | 一种火星dem制作和航带法空中三角测量方法 | |
| CN112526624A (zh) | 重力卫星东西方向差分观测数据构建及反演方法和系统 | |
| CN113936061B (zh) | 海上动态目标定位系统及其定位方法 | |
| JP7451461B2 (ja) | 広域撮像方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |