CN111474544A - 一种基于sar数据的滑坡形变监测及预警方法 - Google Patents

一种基于sar数据的滑坡形变监测及预警方法 Download PDF

Info

Publication number
CN111474544A
CN111474544A CN202010143147.1A CN202010143147A CN111474544A CN 111474544 A CN111474544 A CN 111474544A CN 202010143147 A CN202010143147 A CN 202010143147A CN 111474544 A CN111474544 A CN 111474544A
Authority
CN
China
Prior art keywords
deformation
interference
phase
sar
image
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
Application number
CN202010143147.1A
Other languages
English (en)
Other versions
CN111474544B (zh
Inventor
陈建辉
张水华
吴洪平
麦俊义
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Guangdong Mingyuan Survey And Design Co ltd
Original Assignee
Guangdong Mingyuan Survey And Design Co ltd
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Guangdong Mingyuan Survey And Design Co ltd filed Critical Guangdong Mingyuan Survey And Design Co ltd
Priority to CN202010143147.1A priority Critical patent/CN111474544B/zh
Publication of CN111474544A publication Critical patent/CN111474544A/zh
Application granted granted Critical
Publication of CN111474544B publication Critical patent/CN111474544B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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/9004—SAR image acquisition techniques
    • G01S13/9005—SAR image acquisition techniques with optical processing of the SAR signals
    • 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/9021—SAR image post-processing techniques
    • G01S13/9023—SAR image post-processing techniques combined with interferometric techniques
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B7/00—Measuring arrangements characterised by the use of electric or magnetic techniques
    • G01B7/16—Measuring arrangements characterised by the use of electric or magnetic techniques for measuring the deformation in a solid, e.g. by resistance strain gauge

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Signal Processing (AREA)
  • Radar Systems Or Details Thereof (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明涉及一种基于SAR数据的滑坡形变监测及预警方法。该方法包括：获取覆盖区域的SAR影像；公共主影像选取；SAR影像配准；干涉对配对；差分干涉图生成；筛选干涉对；去除平地相位和地形相位；去除大气相位及残差；建立模型；反演沉降速率和时序形变。本发明针对空气中的水气含量会影响SAR卫星获取数据时的传播路径从而带来较大误差的情况，将覆盖实验区域的降雨量数据按一定比例加入到误差影像因子当中，并建立模型利用最小二乘的方法反演出SAR数据的形变速率，得到实验区域在研究时间段内的时间序列形变，减少因暴雨等因素带来的监测误差，提高InSAR技术获取地表形变的精度，为滑坡形变的监测和预警提供数据和决策支持。

Description

一种基于SAR数据的滑坡形变监测及预警方法
技术领域
本发明涉及地质灾害监测技术领域，具体为一种基于SAR数据的滑坡形变监测及预 警方法。
背景技术
近些年来，随着人口的急剧膨胀和城市化进程的不断加快，城市面积的扩张带来了 山区、森林面积的减小和自然、地理环境的破坏，这给我们身边的自然资源、生存环境都带 来了极大的挑战。而地质灾害的频发便是环境恶化的重要标志之一，滑坡灾害便是最为典型 的地质灾害。中国是世界上人口最多的国家，且国土面积居世界第三位，密集的人口和频繁 的人类活动给中国带来了严重的影响，致使中国成为全球地质灾害最多、影响最严重的国家 之一。因此，如何对滑坡灾害影响范围及破坏程度进行评估以及如何进行大范围的监测以做 到滑坡预警是目前各省、市、县政府的一个重要关注方向。
常规的监测手段虽然在精度上能满足滑坡监测的要求，但是均存在缺陷使得这些技 术在地质灾害监测上具有较大的局限性。全站仪测量和水准测量对监测人员现场操作的技能 要求较高，且在容易发生滑坡的山区、陡坡地段监测精度低；GPS测量虽然可以不用专业人 员现场作业，但是其垂直形变精度远低于水平形变精度；深度位移计作为地表以下形变监测 确实具有优势，但是其费用较高，且点位密度小，容易受到灾害体自身形变带来的破坏，难 以实现长时间持续观测。而这些技术共同的缺陷在于它们均只能针对单点或多点进行监测， 很难实现面状观测，而要实行近似面状的观测则需要投入大量的财力、物力、人力。
合成孔径雷达干涉测量技术(InterferometrySyntheticApertureRadar,InSAR)是一门 在新型的空间对地观测技术，它将干涉测量技术与雷达成像技术进行结合，根据电磁波测距 的原理，利用SAR卫星在运动过程中的轨道参数、天线几何参数以及监测区域的数字高程 模型(Digital ElevationModel，DEM)等要素反映卫星在相邻两幅影像获取过程中的真实地 表形变，精度为毫米级精度。与常规的监测方法相比，InSAR技术具有全方位、多层次、宽 领域的优势。它不受云雾雨雪的影响，同时也不需要专业技术人员到现场进行勘测；同时， 它可以实现大范围监测，SAR影像覆盖范围可到50km2以上，减少了数据获取的成本，提 高了安全性，这也使得InSAR技术在滑坡形变监测中得到了重视和广泛应用。
差分合成孔径雷达干涉测量技术(DifferentialInterferometrySyntheticApertureRadar, D-InSAR)经过近20年的发展已经成为了大面积区域地面沉降监测的重要手段，现如今在 地表沉降、地质灾害监测、地下工程以及火山监测等方向都有着广泛的应用和发展空间。然 而D-InSAR对SAR图像的精度及其配准等要求非常高，通常使用传统的D-InSAR技术需要 满足SAR数据的时空基线小、受大气等影响因素小、DEM精度高等条件才能得到精度更高 的结果，并且不能反映出区域在某时间范围内的地面时序形变规律。随着科技的进步，针对 D-InSAR中存在的诸多问题，科研工作者相继提出了时序InSAR技术。其中，具有代表性 的时序InSAR技术有永久散射体干涉测量(Permanent ScatterInterferometry,PS)技术和小 基线集(Small Baseline Subset,SBAS)技术。
Berardina和Lanari等人于2002年首次提出了SBAS技术，它是将同一个时间序列中 得到的SAR图像进行组合，在组合中需要综合考虑研究区背景、SAR影像时间分布以及DEM等数据的限制条件并得到时空基线的最优解。该组合可以得到若干个集合，集合中的SAR图像都是非失相干的，即SAR影像间相干性普遍较好。同时，SBAS技术利用多视处 理降低差分干涉图的相位噪声，并应用奇异值分解(Singular Value Decomposition,SVD)方法对相位进行求解得到形变相位速率，进而解得整个观测时段的形变时间序列。在山区形变 监测中，尽管PS与SBAS技术都可以反映某区域的形变时间序列，但是与PS技术相比较而言，SBAS技术在山区形变监测上具有显著的优势：1、增加了干涉对的数量；2、PS技术 需要至少25景影像，而SBAS需要的影像更少，在某些影像不足的区域也能使用；3、增加 时空基线长短的限制，使得失相关现象减少，并增加了数据间的约束性；4、增加了对大 气、地形相位及其他误差的分析，测量精度高；5、山区中稳定散射体较少，PS点的缺乏会 影响PS技术的结果，而SBAS技术则不需要PS点的结果作为数据支撑。
目前，虽然SBAS技术可以监测到地表时间序列形变，从而对山体滑坡、崩塌、地 裂缝、地面塌陷等引起地表变化的地质灾害进行预警。但由于SAR卫星数据的不连续性， 使得不同卫星和同一卫星在不同轨道下获取的影像不能进行干涉。时序InSAR技术也只能 通过获取相同地区在相同情况下的SAR影像进行处理得到时间序列形变，这也限制了时序InSAR技术的时间分辨率。
发明内容
为了提高时序InSAR技术的时序分辨率，反应出某区域在长时间序列下的形变规律，本发明提供了一种基于SAR数据的滑坡形变监测及预警方法。该方法基于InSAR技术 和覆盖该区域一段时间内的降雨量数据，采用SBAS方法削弱大气路径效应的干扰，获得高 精度的地表形变，采用的技术方案是，包括以下步骤：
步骤一：获取覆盖监测范围的SAR影像数据；
步骤二：公共主影像选择；
步骤三：SAR影像的配准；
步骤四：干涉对配对；
步骤五：生成差分干涉图；
步骤六：差分干涉图的筛选；
步骤七：去除恒定相位和相位跳变；
步骤八：去除大气相位和残余地形相位；
步骤九：反演形变速率和时序形变量；
在步骤二中，公共主影像的选取需要计算由(1)式得出的任意两幅图像生成干涉图的相干系数ρi,j：
Figure BDA0002399802040000031
其中，Ti、Si、Di分别为时间基线、空间垂直基线以及多普勒质心频率的阈值；通过 (2)式求和，计算出第i幅图像作为公共影像时的相干系数ρi，当ρi最大时Ti、Si、Di达 到最佳组合；
Figure BDA0002399802040000032
进一步的，步骤三中，依据图像相关系数的统计特性，以两幅影像相关系数的最大值作为条件进行配准，当达到最大值时，认为同一区域的两幅影像中所有的像元都精确配准，若配准精度无法达到1/8个像元则无法达到精确配准要求，出现该情况的两幅影像无法相互组成干涉对进行干涉，配准步骤具体如下:
确定同名点和搜索窗口之后，使用大小为R×C的窗口在区域内利用(3)式计算相关系数：
Figure BDA0002399802040000033
其中，x1(r,c)和x2(r,c)是窗口在对应位置(k,l)上的两个复数数据，
Figure BDA0002399802040000034
表示为 x2(r,c)的共轭，在区域内计算每一个搜索点的α值后取最大值对应的位置作为配准点，同时α也与信噪比(Signal toNoise,SNR)有着(4)式对应关系：
Figure BDA0002399802040000035
式(4)说明，相关性越高的相对其信噪比越大；相关性越低的干涉图信噪比越小，所 以同样的将SNR作为衡量相关性的标准。
进一步的，步骤四中，针对传感器、观测条件、研究区域条件的不同设定空间基线阈值和时间基线阈值以及平均相干性阈值能够剔除不符合要求的干涉对，生成M对小基线干涉对，干涉对数量应满足式(5)：
Figure BDA0002399802040000041
在地形起伏较大区域，需借助外部DEM数据去除地形起伏导致的地形相位影响，分别得到各差分干涉图中的形变相位。
进一步的，在干涉时，采用Goldstein方法进行滤波以及Delaunay最小费用流法进行 解缠，同时将影像按比例多视，按一定的比值将横向和纵向的分辨率调成近似，以增加干涉 图的信噪比，提供更为可靠的相干性，对于第j景差分干涉图中，方位向坐标A以及距离向 坐标R的像元的干涉相位值用以下的方式表示：
Figure BDA0002399802040000042
式(6)中，j为影像的编号，范围为(1，……，N-1)，λ为景号信号的中心波长， d(tB,A,R)，d(tA,A,R)为A、B两个时刻对于雷达视线方向(LightofSight,LOS)的累计形 变量，Δh为高程误差，B⊥为垂直基线长，θ表示为入射角，
Figure BDA0002399802040000043
为大气延迟 相位，
Figure BDA0002399802040000044
则表示模型的总的噪声分量。
进一步的，剔除掉干涉图中干涉结果较差的部分，消除此类干涉图带来的误差对结 果造成的影响，同时对于认为相干性和干涉条件比较好但并未出现在干涉图集中的干涉对， 通过干涉和手动添加的方式添加的干涉图集中。
进一步的，采用选取地面控制点(Ground Control Point,GCP)，选择地面点进行轨 道精炼需要满足：不能位于残余地形相位上；如果有已知的形变区域，地面点不能位于形变 区域上，需远离形变区域。
进一步的，为了去除DEM误差相位和大气延迟的带来的影响，以获得精确的时间序列形变，假设地表形变的低频部分为：
Figure BDA0002399802040000045
式(7)中
Figure BDA0002399802040000046
分别表示LOS方向形变的平均速率、平均加速度和平均加速度 变化率，将(7)代入(6)得到式(8)：
Figure BDA0002399802040000047
其中，
Figure BDA0002399802040000048
代表大气延迟相位，
Figure BDA0002399802040000049
代表噪声相位，而
Figure BDA00023998020400000410
则表示由于降雨量的大小引起的传播路径误差相位，将这些相位统一规划为残余差分相位Δφ，则 (8)式写成(9)式的矩阵形式：
Ax＝Δφ (9)
其中
Figure BDA0002399802040000051
而系数矩阵A表示成(7)式的形式：
Figure BDA0002399802040000052
使用最小二乘的方法计算得到高程误差Δh和残余差分相位Δφ，而在Δφ中，
Figure BDA0002399802040000053
在时间域上表示为高频信号，在空间域上表示为低频信号；而
Figure BDA0002399802040000054
则表示为高频信号，
Figure BDA0002399802040000055
是外部因素，通过设置降雨量的阈值，同时将当天降雨量与阈值的差值作为误差来源，再乘以降雨量大小对误差影响的比重得出。
进一步的，在去掉残余地形相位、大气相位和相干噪声后，将(6)式简化为(11)式：
Figure BDA0002399802040000056
用矩阵形式表达(11)式写成：
Figure BDA0002399802040000057
式中含有M个观测向量和N个未知数，如果M≥N时，使用最小二乘法解得结果，然 而在差分干涉中，为了抑制时空基线过长带来的去相干，常出现监测时间段内出现多个小基 线集的情况，该情况下干涉对是不连续的，因此使用最小二乘方法解决此类秩亏问题时解不唯一，此时使用奇异值分解法(SingularValueDecomposition，SVD)得到未知参数φ的最小范数解，在使用SVD方法时，将形变的相位表示成影像时间段的速率与时间的乘积，从 而得到形变相位速率的最小范数解，以得到各个时段的形变量。
本发明的有益效果：本发明利用两组或两组以上具有相同覆盖区域、在不同时间段 内的SAR数据，采用SBAS技术分别进行处理。由于在数据选取时不同的SAR卫星参数不同，其时间序列形变结果需要叠加到垂直方向，便于比较和整合。同时，考虑到在重叠的时间域中各组SAR影像中相同位置的点的变化一致性，在各组数据选取地面点时，尽可能的在选取的地面点中包含多组干涉图中均表现为稳定的点，除能提高结果的稳定性外，还能加 深两组数据之间的联系。最后，在得到各组数据的时序形变结果之后，利用数据间重叠的时 间域中显示出的形变速率的联系，将两组不同SAR影像经SBAS处理的结果整合在一起， 从而提升SBAS结果的时间分辨率，得到更长时间序列的地表形变，为山区形变监测、地表 形变变化规律的发现及地质灾害预防提供有效、可靠的技术支持。
附图说明
图1为本发明的处理流程图；
图2为实验区域图；
图3为实验数据覆盖时间范围内的降雨量折线图；
图4为去除大气误差前的解缠图；
图5为去除后的解缠图；
图6为整体时间序列形变图；
图7为滑坡形变监测位置时序列形变图；
图8为滑坡形变监测位置时序列形变图。
具体实施方式
为使本发明的目的、技术方案及优点更加清晰，以下结合附图及实施实例，对本发明进行进一步阐释。应当理解，此处所描述的具体实施事例仅仅用以解释本发明，并不用于限定本发明。
如图1至图8所示，一种基于SAR数据的滑坡形变监测及预警方法主要包括以下几个步骤：
步骤一：获取覆盖监测范围的SAR影像数据
为便于滑坡形变监测，我们选取滑坡灾害频发的山区地区。同时，考虑到滑坡形变的量级较大，为更好的显示滑坡形变监测结果，我们选取L波段的ALOS数据对四川省理县和 周边山区进行了形变监测，其时间跨度为2007年1月4日～2011年3月2日。数据时间及 其他信息详见表1。
表1 ALOS数据信息表
Figure BDA0002399802040000061
步骤二：公共主影像选择
公共主影像的选取需要计算由(1)式得出的任意两幅图像生成干涉图的相干系数ρi,j：
Figure BDA0002399802040000071
其中，Ti、Si、Di分别为时间基线、空间垂直基线以及多普勒质心频率的阈值。并通 过(2)式求和，计算出第i幅图像作为公共影像时的相干系数ρi。当ρi最大时Ti、Si、Di可达到最佳组合；
Figure BDA0002399802040000072
经过时空基线的综合考虑，本次实验将2008年7月9日的影像作为公共主影像。
步骤三：SAR影像的配准
依据图像相关系数的统计特性，以两幅影像相关系数的最大值作为条件进行配准，当达 到最大值时，可以认为同一区域的两幅影像中所有的像元都精确配准。若配准精度无法达到 1/8个像元则无法达到精确配准要求，出现该情况的两幅影像无法相互组成干涉对进行干 涉。配准步骤具体如下:
确定同名点和搜索窗口之后，使用大小为R×C的窗口在区域内利用(3)式计算相关系数：
Figure BDA0002399802040000073
其中，x1(r,c)和x2(r,c)是窗口在对应位置(k,l)上的两个复数数据，
Figure BDA0002399802040000074
表示为x2(r,c)的共轭。在区域内计算每一个搜索点的α值后取最大值对应的位置作为配准点。同时 α也与信噪比(Signal toNoise,SNR)有着(4)式对应关系：
Figure BDA0002399802040000075
式(4)说明，相关性越高的相对其信噪比越大；相关性越低的干涉图信噪比越小。所 以同样的可以将SNR作为衡量相关性的标准。
步骤四：干涉对配对
针对传感器、观测条件、研究区域等条件的不同设定空间基线阈值和时间基线阈值以及 平均相干性阈值可以剔除不符合要求的干涉对，生成M对小基线干涉对，干涉对数量应满 足式(5)：
Figure BDA0002399802040000076
在地形起伏较大区域，需借助外部DEM数据去除地形起伏导致的地形相位影响，分别得到各差分干涉图中的形变相位。
步骤五：生成差分干涉图
在干涉时，采用Goldstein方法进行滤波以及Delaunay最小费用流法进行解缠，同时将影像按比例多视，按一定的比值将横向和纵向的分辨率调成近似，以增加干涉图的信噪比，提供更为可靠的相干性。对于第j景差分干涉图中，方位向坐标A以及距离向坐标R的像 元的干涉相位值可以用以下的方式表示：
Figure BDA0002399802040000081
式中，j为影像的编号，范围为(1，……，N-1)，λ为景号信号的中心波长， d(tB,A,R)，d(tA,A,R)为A、B两个时刻对于雷达视线方向(LightofSight,LOS)的累计形 变量。Δh为高程误差，B⊥为垂直基线长，θ表示为入射角，
Figure BDA0002399802040000082
为大气延迟相 位，
Figure BDA0002399802040000083
则表示模型的总的噪声分量。
步骤六：差分干涉图的筛选
需要剔除掉干涉图中干涉结果较差的部分，该部分的出现回影响到整个干涉图集的相干 性，从而让连续的相干图在反演后变成离散的散点图。其中影响最终干涉的结果包括：1、 轨道数据不精确；2、相干性低；3、配准参数不正确；4、解缠图显示了严重的不连续性。剔除该部分干涉图，可以消除此类干涉图带来的误差对结果造成的影响，同时对于认为相干 性和干涉条件比较好但并未出现在干涉图集中的干涉对，也可以通过干涉和手动添加的方式 添加的干涉图集中，进一步提高干涉结果的精度和准确性。
经过筛选，本次实验参加干涉的干涉对共有36对。
步骤七：去除恒定相位和相位跳变
通过选取地面控制点(Ground Control Point,GCP)来实现。该步骤是形变监测的关键步骤，SBAS技术得到的结果个地面点的选取有很大的关系。由于SBAS技术中需要对所有 像对的解缠结果进行精炼，所以地面点需要在绝大多数干涉图中都是稳定的。所以选择地面点进行轨道精炼需要注意以下几点：1、不能位于残余地形相位上；2、如果有已知的形变区域，地面点不能位于形变区域上，需远离形变区域。
步骤八：去除大气相位和残余地形相位
为了去除DEM误差相位和大气延迟的带来的影响，以获得精确的时间序列形变，可假设地表形变的低频部分为：
Figure BDA0002399802040000084
式(7)中
Figure BDA0002399802040000091
分别表示LOS方向形变的平均速率、平均加速度和平均加速度 变化率。将(7)代入(6)可以得到式(8)：
Figure BDA0002399802040000092
其中，
Figure BDA0002399802040000093
代表大气延迟相位，
Figure BDA0002399802040000094
代表噪声相位，而
Figure BDA0002399802040000095
则表示由于降雨量的大小引起的传播路径误差相位，可将这些相位统一规划为残余差分相位Δφ， 则(8)式可以写成(9)式的矩阵形式：
Ax＝Δφ (9)
其中
Figure BDA0002399802040000096
而系数矩阵A可以表示成(7)式的形式：
Figure BDA0002399802040000097
使用最小二乘的方法可以计算得到高程误差Δh和残余差分相位Δφ，而在Δφ中，
Figure BDA0002399802040000098
在时间域上表示为高频信号，在空间域上表示为低频信号；而
Figure BDA0002399802040000099
则表示为高频信号。
Figure BDA00023998020400000910
是外部因素，通过设置降雨量的阈值，同时将当天降雨量与阈值的差值作为误差来源，再乘以降雨量大小对误差影响的比重得出。
本次实验时间范围中，由于选取的是ALOS PALSAR-1数据，其重访周期为46天， 因此在一个月中至多只能获取到一幅SAR影像，所以我们将SAR影像获取时所在月的当月 平均降雨量数据作为观测时的降雨量。同时，将公共主影像2008年7月9日的降雨量作为 降雨量的参考阈值。同时，根据根据干涉相对之间的时间间隔、降雨量与总降雨量比值的大 小作为降雨量对大气误差影响中定权的依据，系数可用下式定义：
Figure RE-RE-GDA00025491051500000911
改正前与改正后的图如图4和图5所示。
步骤九：反演形变速率和时序形变量
在去掉残余地形相位、大气相位和相干噪声后，可以将(6)式简化为(11)式
Figure BDA00023998020400000912
用矩阵形式表达(11)式可写成：
Figure BDA00023998020400000913
式中含有M个观测向量和N个未知数，如果M≥N时，可以使用最小二乘法解得结 果，然而在差分干涉中，为了抑制时空基线过长带来的去相干，常出现监测时间段内出现多 个小基线集的情况，该情况下干涉对是不连续的，因此使用最小二乘方法解决此类秩亏问题 时解不唯一，此时应使用奇异值分解法(SingularValue Decomposition，SVD)得到未知参 数φ的最小范数解。在使用SVD方法时，通常将形变的相位表示成影像时间段的速率与时 间的乘积，从而得到形变相位速率的最小范数解，以得到各个时段的形变量。
本发明在利用SAR影像进行SBAS技术处理时，根据选取的SAR影像时间和空间的特性，选取满足要求的覆盖监测范围的降雨量数据(日均、月均或年均)，考虑到降雨量大小对传播路径的影响，将降雨量数据按比例计入到大气误差当中去。与现有技术进行比较，本发明引入了其他可以消除误差的因素，从而提高结果的精度，干涉结果的质量更好，更能减少暴雨带来的失真效果，更能反应地表真实性变。

Claims (8)

1.一种基于SAR数据的滑坡形变监测及预警方法，其特征在于，包括以下步骤：
步骤一：获取覆盖监测范围的SAR影像数据；
步骤二：公共主影像选择；
步骤三：SAR影像的配准；
步骤四：干涉对配对；
步骤五：生成差分干涉图；
步骤六：差分干涉图的筛选；
步骤七：去除恒定相位和相位跳变；
步骤八：去除大气相位和残余地形相位；
步骤九：反演形变速率和时序形变量；
在步骤二中，公共主影像的选取需要计算由(1)式得出的任意两幅图像生成干涉图的相干系数ρi,j：
Figure FDA0002399802030000011
其中，Ti、Si、Di分别为时间基线、空间垂直基线以及多普勒质心频率的阈值；通过(2)式求和，计算出第i幅图像作为公共影像时的相干系数ρi，当ρi最大时Ti、Si、Di达到最佳组合；
Figure FDA0002399802030000012
2.根据权利要求1所述的基于SAR数据的滑坡形变监测及预警方法，其特征在于，步骤三中，依据图像相关系数的统计特性，以两幅影像相关系数的最大值作为条件进行配准，当达到最大值时，认为同一区域的两幅影像中所有的像元都精确配准，若配准精度无法达到1/8个像元则无法达到精确配准要求，出现该情况的两幅影像无法相互组成干涉对进行干涉，配准步骤具体如下:
确定同名点和搜索窗口之后，使用大小为R×C的窗口在区域内利用(3)式计算相关系数：
Figure FDA0002399802030000013
其中，x1(r,c)和x2(r,c)是窗口在对应位置(k,l)上的两个复数数据，
Figure FDA0002399802030000021
表示为x2(r,c)的共轭，在区域内计算每一个搜索点的α值后取最大值对应的位置作为配准点，同时α也与信噪比(Signal toNoise,SNR)有着(4)式对应关系：
Figure FDA0002399802030000022
式(4)说明，相关性越高的相对其信噪比越大；相关性越低的干涉图信噪比越小，所以同样的将SNR作为衡量相关性的标准。
3.根据权利要求1所述的基于SAR数据的滑坡形变监测及预警方法，其特征在于，步骤四中，针对传感器、观测条件、研究区域条件的不同设定空间基线阈值和时间基线阈值以及平均相干性阈值能够剔除不符合要求的干涉对，生成M对小基线干涉对，干涉对数量应满足式(5)：
Figure FDA0002399802030000023
在地形起伏较大区域，需借助外部DEM数据去除地形起伏导致的地形相位影响，分别得到各差分干涉图中的形变相位。
4.根据权利要求1所述的基于SAR数据的滑坡形变监测及预警方法，其特征在于，在干涉时，采用Goldstein方法进行滤波以及Delaunay最小费用流法进行解缠，同时将影像按比例多视，按一定的比值将横向和纵向的分辨率调成近似，以增加干涉图的信噪比，提供更为可靠的相干性，对于第j景差分干涉图中，方位向坐标A以及距离向坐标R的像元的干涉相位值用以下的方式表示：
Figure FDA0002399802030000024
式(6)中，j为影像的编号，范围为(1，……，N-1)，λ为景号信号的中心波长，d(tB,A,R)，d(tA,A,R)为A、B两个时刻对于雷达视线方向(LightofSight,LOS)的累计形变量，Δh为高程误差，B⊥为垂直基线长，θ表示为入射角，
Figure FDA0002399802030000025
为大气延迟相位，
Figure FDA0002399802030000026
则表示模型的总的噪声分量。
5.根据权利要求1所述的基于SAR数据的滑坡形变监测及预警方法，其特征在于，剔除掉干涉图中干涉结果较差的部分，消除此类干涉图带来的误差对结果造成的影响，同时对于认为相干性和干涉条件比较好但并未出现在干涉图集中的干涉对，通过干涉和手动添加的方式添加的干涉图集中。
6.根据权利要求1所述的基于SAR数据的滑坡形变监测及预警方法，其特征在于，采用选取地面控制点(Ground Control Point,GCP)，选择地面点进行轨道精炼需要满足：不能位于残余地形相位上；如果有已知的形变区域，地面点不能位于形变区域上，需远离形变区域。
7.根据权利要求4所述的基于SAR数据的滑坡形变监测及预警方法，其特征在于，为了去除DEM误差相位和大气延迟的带来的影响，以获得精确的时间序列形变，假设地表形变的低频部分为：
Figure FDA0002399802030000031
式(7)中
Figure FDA0002399802030000032
分别表示LOS方向形变的平均速率、平均加速度和平均加速度变化率，将(7)代入(6)得到式(8)：
Figure FDA0002399802030000033
其中，
Figure FDA0002399802030000034
代表大气延迟相位，
Figure FDA0002399802030000035
代表噪声相位，而
Figure FDA0002399802030000036
则表示由于降雨量的大小引起的传播路径误差相位，将这些相位统一规划为残余差分相位Δφ，则(8)式写成(9)式的矩阵形式：
Ax＝Δφ (9)
其中
Figure FDA0002399802030000037
而系数矩阵A表示成(7)式的形式：
Figure FDA0002399802030000038
使用最小二乘的方法计算得到高程误差Δh和残余差分相位Δφ，而在Δφ中，
Figure FDA0002399802030000039
在时间域上表示为高频信号，在空间域上表示为低频信号；而
Figure FDA00023998020300000310
则表示为高频信号，
Figure FDA00023998020300000311
是外部因素，通过设置降雨量的阈值，同时将当天降雨量与阈值的差值作为误差来源，再乘以降雨量大小对误差影响的比重得出。
8.根据权利要求4所述的基于SAR数据的滑坡形变监测及预警方法，其特征在于，在去掉残余地形相位、大气相位和相干噪声后，将(6)式简化为(11)式：
Figure FDA00023998020300000312
用矩阵形式表达(11)式写成：
Figure FDA00023998020300000313
式中含有M个观测向量和N个未知数，如果M≥N时，使用最小二乘法解得结果，然而在差分干涉中，为了抑制时空基线过长带来的去相干，常出现监测时间段内出现多个小基线集的情况，该情况下干涉对是不连续的，因此使用最小二乘方法解决此类秩亏问题时解不唯一，此时使用奇异值分解法(Singular Value Decomposition，SVD)得到未知参数φ的最小范数解，在使用SVD方法时，将形变的相位表示成影像时间段的速率与时间的乘积，从而得到形变相位速率的最小范数解，以得到各个时段的形变量。
CN202010143147.1A 2020-03-04 2020-03-04 一种基于sar数据的滑坡形变监测及预警方法 Active CN111474544B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010143147.1A CN111474544B (zh) 2020-03-04 2020-03-04 一种基于sar数据的滑坡形变监测及预警方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010143147.1A CN111474544B (zh) 2020-03-04 2020-03-04 一种基于sar数据的滑坡形变监测及预警方法

Publications (2)

Publication Number Publication Date
CN111474544A true CN111474544A (zh) 2020-07-31
CN111474544B CN111474544B (zh) 2022-11-18

Family

ID=71748027

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010143147.1A Active CN111474544B (zh) 2020-03-04 2020-03-04 一种基于sar数据的滑坡形变监测及预警方法

Country Status (1)

Country Link
CN (1) CN111474544B (zh)

Cited By (25)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112068136A (zh) * 2020-09-14 2020-12-11 广东省核工业地质局测绘院 一种基于幅度偏移量的方位向形变监测方法
CN112146622A (zh) * 2020-10-23 2020-12-29 湖南航天智远科技有限公司 输电线路沿线地质沉降监测方法
CN112363165A (zh) * 2020-11-19 2021-02-12 中南大学 一种林下地形反演方法、装置、设备及介质
CN112505068A (zh) * 2020-11-03 2021-03-16 桂林理工大学 一种基于gnss-ir的地表土壤湿度多星组合反演方法
CN112698328A (zh) * 2020-11-30 2021-04-23 四川大学 一种用于大坝及滑坡变形gb-sar监测的相位解缠方法及系统
CN112946647A (zh) * 2021-02-02 2021-06-11 河海大学 大气误差改正InSAR干涉图堆叠地质灾害普查方法和装置
CN112986993A (zh) * 2021-02-07 2021-06-18 同济大学 一种基于空间约束的InSAR形变监测方法
CN113281744A (zh) * 2021-03-11 2021-08-20 中南大学 基于假设检验和自适应形变模型的时序InSAR方法
CN113281742A (zh) * 2021-06-02 2021-08-20 西南交通大学 一种基于滑坡形变信息和气象数据的sar滑坡预警方法
CN113790667A (zh) * 2021-11-18 2021-12-14 中大检测(湖南)股份有限公司 基于雷达的大坝变形检测方法
CN113848551A (zh) * 2021-09-24 2021-12-28 成都理工大学 一种利用InSAR升降轨形变数据的滑坡深度反演方法
CN114001703A (zh) * 2021-10-09 2022-02-01 四川轻化工大学 一种滑坡变形数据实时过滤方法
CN114111654A (zh) * 2021-12-06 2022-03-01 国网湖南省电力有限公司 基于DS-InSAR技术监测输电通道附近滑坡的方法及系统
CN114236541A (zh) * 2021-12-08 2022-03-25 电子科技大学 基于Sentinel-1卫星SAR图像的大面积地表三维形变计算方法
CN114440758A (zh) * 2022-01-09 2022-05-06 西北大学 一种区域尺度上滑坡对降雨响应的分析方法
CN114594479A (zh) * 2022-05-07 2022-06-07 中国测绘科学研究院 一种全散射体FS-InSAR方法及系统
CN114812491A (zh) * 2022-01-19 2022-07-29 国网甘肃省电力公司经济技术研究院 基于长时间序列分析的输电线路地表形变预警方法及装置
CN115201825A (zh) * 2022-09-16 2022-10-18 眉山环天智慧科技有限公司 一种InSAR震间形变监测中的大气延迟校正方法
CN115512531A (zh) * 2022-09-28 2022-12-23 重庆地质矿产研究院 一种基于形变有序性的滑坡灾害多监测点融合预警方法
CN115856889A (zh) * 2023-02-03 2023-03-28 北京数慧时空信息技术有限公司 一种误差自动校正的InSAR时序变形监测方法
CN115993601A (zh) * 2023-03-22 2023-04-21 四川省公路规划勘察设计研究院有限公司 一种强盐渍土区域公路变形的时序InSAR监测方法
CN116148855A (zh) * 2023-04-04 2023-05-23 之江实验室 时序InSAR大气相位去除和形变解算的方法及系统
CN116258005A (zh) * 2023-02-23 2023-06-13 哈尔滨工业大学（深圳） 一种基于InSAR数据的桥梁形变预测方法及相关装置
CN117031425A (zh) * 2023-06-26 2023-11-10 昆明理工大学 基于星载sar的大梯度滑坡失稳时间预测方法及系统
CN117274342A (zh) * 2023-11-21 2023-12-22 中铁水利水电规划设计集团有限公司 一种基于卫星数据的水利工程形变监测方法

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2413158A1 (en) * 2010-07-26 2012-02-01 Consorci Institut de Geomatica A method for monitoring terrain and man-made feature displacements using ground-based synthetic aperture radar (GBSAR) data
CN106950556A (zh) * 2017-05-03 2017-07-14 三亚中科遥感研究所 基于分布式散射体时序干涉sar技术的遗产地形变监测方法
CN107132539A (zh) * 2017-05-03 2017-09-05 中国地质科学院探矿工艺研究所 一种基于小基线集的时间序列InSAR的滑坡早期识别方法
CN108387899A (zh) * 2018-04-17 2018-08-10 南京师范大学 合成孔径雷达干涉测量中地面控制点自动选取方法
CN108627832A (zh) * 2018-05-11 2018-10-09 电子科技大学 一种基于多时序sar图像提取输电通道地表形变的方法
CN108663017A (zh) * 2018-08-13 2018-10-16 伟志股份公司 一种监测城市地铁沿线地表沉降的方法
CN108957456A (zh) * 2018-08-13 2018-12-07 伟志股份公司 基于多数据源sbas技术的滑坡监测和早期识别方法
CN110174044A (zh) * 2019-04-16 2019-08-27 中国矿业大学 一种基于psi技术的桥梁纵向位移形变监测的方法
CN110673145A (zh) * 2019-10-24 2020-01-10 中国地质大学（北京） 一种基于间断相干的InSAR地表形变监测方法及系统

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2413158A1 (en) * 2010-07-26 2012-02-01 Consorci Institut de Geomatica A method for monitoring terrain and man-made feature displacements using ground-based synthetic aperture radar (GBSAR) data
CN106950556A (zh) * 2017-05-03 2017-07-14 三亚中科遥感研究所 基于分布式散射体时序干涉sar技术的遗产地形变监测方法
CN107132539A (zh) * 2017-05-03 2017-09-05 中国地质科学院探矿工艺研究所 一种基于小基线集的时间序列InSAR的滑坡早期识别方法
CN108387899A (zh) * 2018-04-17 2018-08-10 南京师范大学 合成孔径雷达干涉测量中地面控制点自动选取方法
CN108627832A (zh) * 2018-05-11 2018-10-09 电子科技大学 一种基于多时序sar图像提取输电通道地表形变的方法
CN108663017A (zh) * 2018-08-13 2018-10-16 伟志股份公司 一种监测城市地铁沿线地表沉降的方法
CN108957456A (zh) * 2018-08-13 2018-12-07 伟志股份公司 基于多数据源sbas技术的滑坡监测和早期识别方法
CN110174044A (zh) * 2019-04-16 2019-08-27 中国矿业大学 一种基于psi技术的桥梁纵向位移形变监测的方法
CN110673145A (zh) * 2019-10-24 2020-01-10 中国地质大学（北京） 一种基于间断相干的InSAR地表形变监测方法及系统

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
刘广全: "基于SBAS-InSAR 的丹巴县滑坡探测与监测", 《中国优秀博硕士学位论文全文数据库(硕士)基础科学辑》 *
周琦: "基于SBAS-InSAR技术的北京平原区地面沉降监测研究", 《中国优秀博硕士学位论文全文数据库(硕士)基础科学辑》 *
姜德才: "重轨卫星InSAR技术在多云雨山区DEM生成及更新中的应用研究", 《中国优秀博硕士学位论文全文数据库(硕士)基础科学辑》 *

Cited By (37)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112068136A (zh) * 2020-09-14 2020-12-11 广东省核工业地质局测绘院 一种基于幅度偏移量的方位向形变监测方法
CN112146622A (zh) * 2020-10-23 2020-12-29 湖南航天智远科技有限公司 输电线路沿线地质沉降监测方法
CN112505068B (zh) * 2020-11-03 2023-08-11 桂林理工大学 一种基于gnss-ir的地表土壤湿度多星组合反演方法
CN112505068A (zh) * 2020-11-03 2021-03-16 桂林理工大学 一种基于gnss-ir的地表土壤湿度多星组合反演方法
CN112363165A (zh) * 2020-11-19 2021-02-12 中南大学 一种林下地形反演方法、装置、设备及介质
CN112363165B (zh) * 2020-11-19 2022-06-24 中南大学 一种林下地形反演方法、装置、设备及介质
CN112698328A (zh) * 2020-11-30 2021-04-23 四川大学 一种用于大坝及滑坡变形gb-sar监测的相位解缠方法及系统
CN112698328B (zh) * 2020-11-30 2021-08-10 四川大学 一种用于大坝及滑坡变形gb-sar监测的相位解缠方法及系统
CN112946647A (zh) * 2021-02-02 2021-06-11 河海大学 大气误差改正InSAR干涉图堆叠地质灾害普查方法和装置
CN112986993A (zh) * 2021-02-07 2021-06-18 同济大学 一种基于空间约束的InSAR形变监测方法
CN113281744A (zh) * 2021-03-11 2021-08-20 中南大学 基于假设检验和自适应形变模型的时序InSAR方法
CN113281742A (zh) * 2021-06-02 2021-08-20 西南交通大学 一种基于滑坡形变信息和气象数据的sar滑坡预警方法
CN113848551A (zh) * 2021-09-24 2021-12-28 成都理工大学 一种利用InSAR升降轨形变数据的滑坡深度反演方法
CN113848551B (zh) * 2021-09-24 2023-03-21 成都理工大学 一种利用InSAR升降轨形变数据的滑坡深度反演方法
CN114001703A (zh) * 2021-10-09 2022-02-01 四川轻化工大学 一种滑坡变形数据实时过滤方法
CN114001703B (zh) * 2021-10-09 2023-07-28 四川轻化工大学 一种滑坡变形数据实时过滤方法
CN113790667A (zh) * 2021-11-18 2021-12-14 中大检测(湖南)股份有限公司 基于雷达的大坝变形检测方法
CN114111654A (zh) * 2021-12-06 2022-03-01 国网湖南省电力有限公司 基于DS-InSAR技术监测输电通道附近滑坡的方法及系统
CN114236541A (zh) * 2021-12-08 2022-03-25 电子科技大学 基于Sentinel-1卫星SAR图像的大面积地表三维形变计算方法
CN114236541B (zh) * 2021-12-08 2023-05-16 电子科技大学 基于Sentinel-1卫星SAR图像的大面积地表三维形变计算方法
CN114440758A (zh) * 2022-01-09 2022-05-06 西北大学 一种区域尺度上滑坡对降雨响应的分析方法
CN114812491B (zh) * 2022-01-19 2023-10-20 国网甘肃省电力公司经济技术研究院 基于长时间序列分析的输电线路地表形变预警方法及装置
CN114812491A (zh) * 2022-01-19 2022-07-29 国网甘肃省电力公司经济技术研究院 基于长时间序列分析的输电线路地表形变预警方法及装置
CN114594479B (zh) * 2022-05-07 2022-07-26 中国测绘科学研究院 一种全散射体FS-InSAR方法及系统
CN114594479A (zh) * 2022-05-07 2022-06-07 中国测绘科学研究院 一种全散射体FS-InSAR方法及系统
CN115201825B (zh) * 2022-09-16 2023-01-17 眉山环天智慧科技有限公司 一种InSAR震间形变监测中的大气延迟校正方法
CN115201825A (zh) * 2022-09-16 2022-10-18 眉山环天智慧科技有限公司 一种InSAR震间形变监测中的大气延迟校正方法
CN115512531A (zh) * 2022-09-28 2022-12-23 重庆地质矿产研究院 一种基于形变有序性的滑坡灾害多监测点融合预警方法
CN115856889A (zh) * 2023-02-03 2023-03-28 北京数慧时空信息技术有限公司 一种误差自动校正的InSAR时序变形监测方法
CN116258005B (zh) * 2023-02-23 2023-11-21 哈尔滨工业大学（深圳） 一种基于InSAR数据的桥梁形变预测方法及相关装置
CN116258005A (zh) * 2023-02-23 2023-06-13 哈尔滨工业大学（深圳） 一种基于InSAR数据的桥梁形变预测方法及相关装置
CN115993601A (zh) * 2023-03-22 2023-04-21 四川省公路规划勘察设计研究院有限公司 一种强盐渍土区域公路变形的时序InSAR监测方法
CN116148855A (zh) * 2023-04-04 2023-05-23 之江实验室 时序InSAR大气相位去除和形变解算的方法及系统
CN117031425A (zh) * 2023-06-26 2023-11-10 昆明理工大学 基于星载sar的大梯度滑坡失稳时间预测方法及系统
CN117031425B (zh) * 2023-06-26 2024-06-07 昆明理工大学 基于星载sar的大梯度滑坡失稳时间预测方法及系统
CN117274342A (zh) * 2023-11-21 2023-12-22 中铁水利水电规划设计集团有限公司 一种基于卫星数据的水利工程形变监测方法
CN117274342B (zh) * 2023-11-21 2024-02-13 中铁水利水电规划设计集团有限公司 一种基于卫星数据的水利工程形变监测方法

Also Published As

Publication number Publication date
CN111474544B (zh) 2022-11-18

Similar Documents

Publication Publication Date Title
CN111474544B (zh) 一种基于sar数据的滑坡形变监测及预警方法
Zebker et al. On the derivation of coseismic displacement fields using differential radar interferometry: The Landers earthquake
Hagberg et al. Repeat-pass SAR interferometry over forested terrain
CN110888130A (zh) 一种基于升降轨时序InSAR的煤矿区地表形变监测方法
JP2003500658A (ja) 市街地領域及び地滑り地帯の運動に関するレーダー測定のための手順
CN112284332B (zh) 基于高分辨率insar的高层建筑沉降监测结果三维定位方法
Chen et al. Monitoring earth surface deformations with InSAR technology: principles and some critical issues
CN112444188B (zh) 一种多视角InSAR海堤高精度三维形变测量方法
Zhang et al. Deformations monitoring in complicated-surface areas by adaptive distributed Scatterer InSAR combined with land cover: Taking the Jiaju landslide in Danba, China as an example
De Oliveira et al. Assessment of radargrammetric DSMs from TerraSAR-X Stripmap images in a mountainous relief area of the Amazon region
Nitti et al. On the use of COSMO/SkyMed data and Weather Models for interferometric DEM generation
CN113238228B (zh) 基于水准约束的三维地表形变获取方法、系统及装置
Mao et al. Estimation and compensation of ionospheric phase delay for multi-aperture InSAR: An azimuth split-spectrum interferometry approach
Michelini et al. Introduction to IBIS-ArcSAR: a circular scanning GB-SAR system for deformation monitoring
Chang et al. InSAR atmospheric distortions mitigation: GPS observations and NCEP FNL data
Wang et al. Bistatic InSAR
Trevoho et al. Monitoring of displacements and deformations of the earth’s surface near the Stebnyk city using radar images of Sentinel-1
Wieczorek EVALUATION OF DEFORMATIONS IN THE URBAN AREA OF OLSZTYN USING SENTINEL-1 SAR INTERFEROMETRY.
Kumar et al. Spaceborne InSAR technique for study of Himalayan glaciers using ENVISAT ASAR and ERS data
Arief et al. A fusion of digital elevation model based on interferometry SAR technique from ascending and descending path in urban area
Wang et al. Deformation monitoring and evaluation of mountain slope stability combined with ground-based radar and spaceborne InSAR methods
Kang et al. Technology on High-accuracy DEM Extraction From Airborne Interferometric SAR
Balaji Estimation and correction of tropospheric and ionospheric effects on differential SAR interferograms
Zamiri Aghdam et al. Monitoring of Urmia Lake Bridge Subsidence during 2014-2021 Using DInSAR-SBAS Method and GPS Data
Rouyet et al. Spaceborne radar interferometry (InSAR) for natural hazards, landslides and infrastructure: limitations and potential/Satellittbasert radarinterferometri (InSAR) for naturfare, skred og infrastruktur: begrensninger og muligheter (5/2015)

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