CN112711022A - 一种GNSS层析技术辅助的InSAR大气延迟改正方法 - Google Patents

一种GNSS层析技术辅助的InSAR大气延迟改正方法 Download PDF

Info

Publication number
CN112711022A
CN112711022A CN202011502257.9A CN202011502257A CN112711022A CN 112711022 A CN112711022 A CN 112711022A CN 202011502257 A CN202011502257 A CN 202011502257A CN 112711022 A CN112711022 A CN 112711022A
Authority
CN
China
Prior art keywords
gnss
refractive index
delay
atmospheric
total refractive
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
CN202011502257.9A
Other languages
English (en)
Other versions
CN112711022B (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.)
China University of Mining and Technology CUMT
Original Assignee
China University of Mining and Technology CUMT
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 China University of Mining and Technology CUMT filed Critical China University of Mining and Technology CUMT
Priority to CN202011502257.9A priority Critical patent/CN112711022B/zh
Publication of CN112711022A publication Critical patent/CN112711022A/zh
Application granted granted Critical
Publication of CN112711022B publication Critical patent/CN112711022B/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/9021—SAR image post-processing techniques
    • G01S13/9023—SAR image post-processing techniques combined with interferometric techniques
    • 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/86—Combinations of radar systems with non-radar systems, e.g. sonar, direction finder

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Investigating Or Analysing Materials By Optical Means (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明提供了一种GNSS层析技术辅助的InSAR大气延迟改正方法，包括如下步骤：S10重构研究区域内大气水汽的三维总折射率场；S20以研究区域的无线电探空数据为参考值，对层析结果进行精度评定，若精度符合要求，进行步骤S30；S30利用空间插值算法获取精细化三维总折射率场；S40利用Ray‑tracing算法估计主影像和辅影像中SAR信号的斜路径总延迟量；S50生成干涉SAR影像的大气延迟相位图
Figure DDA0002843952630000011
S60生成SAR解缠干涉相位图；S70获取大气改正的InSAR相位图。本发明的一种GNSS层析技术辅助的InSAR大气延迟改正方法，利用GNSS层析技术和空间插值算法获取精细化三维总折射率场，并结合Ray‑tracing算法精确估计三维SAR信号的斜路径大气延迟量，提高SAR影像大气延迟相位的反演精度，实现InSAR大气延迟的高精度改正。

Description

一种GNSS层析技术辅助的InSAR大气延迟改正方法
技术领域
本发明涉及InSAR大气延迟改正领域，具体涉及一种GNSS层析技术辅助的InSAR大气延迟改正方法。
背景技术
合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar,InSAR)作为近年来极具应用价值的空间大地测量手段，具有精度高、范围广、高分辨率等优势，在地表高程重建、地面沉降监测、山体滑坡及冰川移动监测等领域具有广阔的应用前景。在实际测量中，SAR卫星信号受大气水汽影响，引起传播路径延迟，即大气延迟效应，这不仅降低了InSAR技术的测量精度，还在一定程度上限制了InSAR技术的应用范围。
InSAR测量中的大气延迟是指两次SAR成像时刻大气折射引起的传播路径的延迟的差值，由于水汽的高时空变化特性，两次成像时刻的大气延迟量并不一致，导致InSAR结果的错误解译。研究表明20％的大气相对湿度变化会导致10-14cm的形变测量误差，因此大气延迟成为制约InSAR精度的主要误差源，准确估计InSAR大气延迟量成为InSAR测量中亟待解决的问题之一。
目前有大量学者致力于InSAR大气延迟改正的研究，主要包括两类方法：一类是利用外部气象数据直接估计大气延迟量，如MERIS和MODIS水汽产品、ERA-Interim再分析资料、WRF数值预报模型，该类方法面临观测时间不一致，空间分辨率不匹配、易受云的影响等问题，降低了其自身的适用性。另一类是利用InSAR数据自身进行大气延迟改正，如PS-InSAR、SBAS-InSAR技术等，该类方法须对三维大气的时空分布构造经验假设模型，适用于垂直分层延迟占主导地位的大气改正，无法消除由雷电、暴风雨等突发气象事件引起的大气延迟影响。此外，随着GNSS技术的快速发展，部分学者利用GNSS数据估计研究区域的对流层总延迟，对InSAR大气延迟进行改正。由于GNSS站点的密度远远低于InSAR影像的分辨率，获取的对流层总延迟空间分辨率较低，无法满足高分辨率InSAR大气改正的需求。
发明内容
为了解决上述问题，本发明提供一种GNSS层析技术辅助的InSAR大气延迟改正方法，利用全天候的GNSS数据，实时提供高时空分辨率的三维总折射率场，有利于获取与SAR影像同步的大气延迟相位图，解决传统方法观测时间不一致的问题；利用空间插值算法获取精细化三维总折射率场，并结合Ray-tracing算法精确估计三维SAR信号的斜路径大气延迟量，避免了信号投影带来的估计误差，提高SAR影像大气延迟相位的精度，实现InSAR大气延迟的高精度改正。
为了实现以上目的，本发明采取的一种技术方案是：
一种GNSS层析技术辅助的InSAR大气延迟改正方法，包括如下步骤：S10根据研究区域内GNSS信号的斜路径总延迟观测量，采用GNSS层析技术重构研究区域内大气水汽的三维总折射率场；S20以研究区域的无线电探空数据为参考值，对层析结果进行精度评定，若精度符合要求，进行步骤S30，否则，重复步骤S10，优化GNSS层析模型的配置与解算策略，重新反演三维总折射率场；S30利用空间插值算法获取精细化三维总折射率场：在对流层上下两部分以不同的空间分辨率对初始三维总折射率场进行插值加密，在水平方向采用双线性插值算法，在垂直方向采用三次样条插值算法计算出加密后每个小体素块的总折射率值，构造精细化三维总折射率场；S40利用Ray-tracing算法估计主影像和辅影像中SAR信号的斜路径总延迟量：根据SAR卫星的轨道信息与SAR影像中每个像元的经纬度坐标信息可构造出多条SAR卫星观测信号，每一个像元与SAR卫星之间的连线相当于一条观测信号，针对SAR卫星观测信号，利用Ray-tracing算法估计出每一条三维SAR信号的斜路径总延迟量；S50将步骤S40计算的SAR主影像和辅影像中SAR信号的斜路径总延迟记为
Figure BDA0002843952610000021
和
Figure BDA0002843952610000022
利用公式
Figure BDA0002843952610000023
分别计算出主影像和辅影像的大气延迟相位
Figure BDA0002843952610000024
和
Figure BDA0002843952610000025
将两者相减生成干涉SAR影像的大气延迟相位图
Figure BDA0002843952610000026
S60对主影像和辅影像进行干涉处理，生成差分干涉图，利用滤波算法对生成的干涉图进行滤波，并对其进行相位解缠；生成SAR解缠干涉相位图；S70获取大气改正的InSAR相位图：对干涉SAR影像的大气延迟相位图进行滤波处理，削弱差分延迟相位中的噪声以及操作误差，利用公式
Figure BDA0002843952610000031
剔除大气延迟效应，获得大气改正后干涉相位图，其中，
Figure BDA0002843952610000032
表示未进行改正的解缠干涉相位图，
Figure BDA0002843952610000033
表示大气延迟改正后的InSAR相位图。
进一步地，所述步骤S10包括如下步骤：S11在研究区域设置多个GNSS接收机获取卫星发射的电磁波信号，所述电磁波信号包括伪距和载波相位观测数据，利用高精度GNSS数据处理软件GAMIT/GLOBK对所述观测数据进行解算，获得GNSS信号的位置信息以及GNSS测站处的气象参数值，其中所述位置信息包括高度角ele、方位角azi；所述气象参数值包括测站天顶对流层延迟ZTD、测站上空东西方向的大气梯度GEW和南北方向的大气梯度GNS；S12采用公式(1)计算GNSS测站处的天顶大气干延迟，其中，ZHD为天顶大气干延迟，PS表示GNSS测站处的地表气压，
Figure BDA0002843952610000034
表示GNSS测站的纬度，H表示GNSS测站的高程；
Figure BDA0002843952610000035
S13根据公式ZWD＝ZTD-ZHD计算出GNSS测站处的天顶大气湿延迟ZWD，同时将东西方向和南北方向的大气梯度分解为两个方向上的干延迟梯度
Figure BDA0002843952610000036
和湿延迟梯度
Figure BDA0002843952610000037
S14采用公式(2)分别计算出GNSS观测信号的斜路径干延迟和斜路径湿延迟，并将两者相加获得斜路径总延迟，STD＝SHD+SWD，STD表示斜路径总延迟；
Figure BDA0002843952610000041
其中，mfh和mfw分别表示干映射函数和湿映射函数，mfg表示水平梯度函数；S15综合GNSS测站的分布情况以及SAR影像的观测范围确定GNSS层析模型的层析区域，将层析区域上空的三维空间进行离散化，在水平方向上按照一定的经度差和纬度差进行划分，在垂直方向上采用非均匀划分方式，形成p个均匀分布的体素块，构成GNSS三维体素块层析模型，其中，p＝n·m·t，n、m、t分别表示体素块的行数、列数和层数；S16建立GNSS三维层析模型的观测方程组和约束方程组：根据每一条GNSS信号与三维体素块的空间位置关系可以构造GNSS层析观测方程
Figure BDA0002843952610000042
其中，STDq表示卫星信号线q的斜路径总延迟量，aijk表示该信号在第k层、第i行、第j列的体素块中的截距，Nijk表示该体素块的总折射率值；多条GNSS信号组成的层析观测方程可用矩阵形式表达STD＝A·Ν，其中STD表示由多条GNSS信号的斜路径总延迟组成的观测值向量，A表示多条信号对应的观测方程组的系数矩阵，N表示总折射率未知参数组成的列向量；同时，添加水平约束条件和垂直约束条件来解算层析观测方程组的病态性问题，所述水平约束条件选用基于高斯加权函数的约束模型，假设同一层面内某一体素块的总折射率与其周围体素块的总折射率存在式(3)的数量关系：
w1N1+w2N2+…+wi-1Ni-1-Ni+wi+1Ni+1+…+wnmNnm＝0 (3)
其中，N1，N2，…，Nnm表示同一层内nm个体素块的总折射率，w1，w2，…，wnm表示对应每个体素块的权系数，
Figure BDA0002843952610000043
dk表示同一层内第k个体素块与体素块i之间的距离，σ为平滑因子；所述垂直约束条件选用基于指数递减函数的约束模型，不同高度层的体素块的总折射率之间存在式(4)所示的指数型关系：
Figure BDA0002843952610000051
其中，hj和hi分别表示第j层和第i层体素块的高度，N(hj)和N(hi)则表示对应体素块的总折射率，Hscale表示研究区域的水汽标高。S17结合GNSS层析观测方程组、水平约束方程和垂直约束方程构建式(5)所示的GNSS三维层析模型，
Figure BDA0002843952610000052
其中，Aobs，AH，AV分别表示GNSS观测信号、水平约束、垂直约束对应的系数矩阵，Δobs，ΔH，ΔV分别表示对应的噪声，X为总折射率值；S18通过代数重构算法(6)对公式(5)进行解算，获得GNSS层析的三维总折射率场，其中，yk和yk+1分别表示第k次和第k+1次迭代的近似解，λ为迭代松弛因子，bi为GNSS水汽层析模型左侧列向量的第i行元素，Ai表示层析模型系数矩阵A的第i行元素组成的行向量，
Figure BDA0002843952610000053
进一步地，所述GNSS接收机的分布密度为15km2～20km2。
进一步地，所述层析区域的范围略大于所述SAR影像的观测范围，所述体素块在水平方向的边长为10～15km，所述体素块的总层数t为10～20层。
进一步地，所述步骤S20中，选用研究区域内部的探空站的观测数据采用式(7)计算出所述探空数据在不同高度处的总折射率的参考值，利用层析解算的结果确定同一高度、同一位置处的总折射率值，然后分别根据式(8)、式(9)计算层析估值与探空参考值的平均偏差和均方根偏差；
Figure BDA0002843952610000061
Figure BDA0002843952610000062
Figure BDA0002843952610000063
其中，p为无线电探空获取的在某一高度的气压，e和T分别表示该高度处的水汽压和温度，Zi代表GNSS层析模型解算出来的体素块的总折射率，Zi'代表探空数据计算出来的体素块的总折射率值，s表示参与精度评定的体素块个数。
进一步地，当平均偏差小于0.5g/m3，均方根偏差小于1.5g/m3判定层析结果质量好，满足SAR信号大气延迟估计的要求。
进一步地，所述步骤S30中，所述双线性插值算法如下式(10)所示：
Figure BDA0002843952610000064
其中，(b,l)为待插值点的坐标，Nwet(b,l)表示待插值点的总折射率值，(b11,l11)，(b21,l21)，(b12,l12)，(b22,l22)分别表示待插值点四周临近体素块的中心坐标，Nwet(b11,l11)，Nwet(b21,l21)，Nwet(b12,l12)，Nwet(b22,l22)则分别对应四周体素块的总折射率值；所述三次样条插值算法如下式(11)所示：
Figure BDA0002843952610000065
其中，h表示待插值点的高度，Nwet(h)表示待插值点的总折射率值，hi，hi-1分别表示插值点所在高度层的上下两个层面的高度，Δhi为该层的层间高度，即hi-hi-1，Nwet(hi)，Nwet(hi-1)为对应两个高度的总折射率值，Nwet(h1)，Nwet(h2)，…，Nwet(hk)则分别表示层析范围内每一高度层处的总折射率值，Di,:，Di-1,:分别表示转换矩阵D的第i列，第i-1列对应的列向量，式(11)中，D＝-A-1B表示转换矩阵，只与垂直分层有关，更详细地，其中A，B的计算公式如式(12)及式(13)所示
Figure BDA0002843952610000071
Figure BDA0002843952610000072
进一步地，所述步骤S40中，所述Ray-tracing算法如下式(14)所示：
Figure BDA0002843952610000073
其中，n＝1+N×10-6表示大气总折射率指数，si为SAR信号穿过精细化三维总折射率场中第i层的距离。
本发明的上述技术方案相比现有技术具有以下优点：
本发明的一种GNSS层析技术辅助的InSAR大气延迟改正方法，利用全天候的GNSS数据，实时提供高时空分辨率的三维总折射率场，有利于获取与SAR影像同步的大气延迟相位图，解决传统方法观测时间不一致的问题；利用空间插值算法获取精细化三维总折射率场，并结合Ray-tracing算法精确估计三维SAR信号的斜路径大气延迟量，避免了信号投影带来的估计误差，提高SAR影像大气延迟相位的精度，实现InSAR大气延迟的高精度改正。
附图说明
下面结合附图，通过对本发明的具体实施方式详细描述，将使本发明的技术方案及其有益效果显而易见。
下面结合附图，通过对本发明的具体实施方式详细描述，将使本发明的技术方案及其有益效果显而易见。
图1所示为本发明一实施例的一种GNSS层析技术辅助的InSAR大气延迟改正方法流程图；
图2所示为本发明一实施例的GNSS层析技术辅助的InSAR大气延迟改正方法流程简图；
图3所示为本发明一实施例的初始三维总折射率场与精细化三维总折射率场示意图；
图4所示为本发明一实施例的GNSS与SAR信号空间二维分布示意图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本实施例提供了一种GNSS层析技术辅助的InSAR大气延迟改正方法，如图1～2所示，包括如下步骤：S10根据研究区域内GNSS信号的斜路径总延迟观测量，采用GNSS层析技术重构研究区域内大气水汽的三维总折射率场。S20以研究区域的无线电探空数据为参考值，对层析结果进行精度评定，若精度符合要求，进行步骤S30，否则，重复步骤S10，优化GNSS层析模型的配置与解算策略，重新反演三维总折射率场。S30利用空间插值算法获取精细化三维总折射率场：在对流层上下两部分以不同的空间分辨率对初始三维总折射率场进行插值加密，在水平方向采用双线性插值算法，在垂直方向采用三次样条插值算法计算出加密后每个小体素块的总折射率值，构造精细化三维总折射率场。S40利用Ray-tracing算法估计主影像和辅影像中SAR信号的斜路径总延迟量：根据SAR卫星的轨道信息与SAR影像中每个像元的经纬度坐标信息可构造出多条SAR卫星观测信号，每一个像元与SAR卫星之间的连线相当于一条观测信号，针对SAR卫星观测信号，利用Ray-tracing算法估计出每一条三维SAR信号的斜路径总延迟量。S50将步骤S40计算的SAR主影像和辅影像中SAR信号的斜路径总延迟记为
Figure BDA0002843952610000091
和
Figure BDA0002843952610000092
利用公式
Figure BDA0002843952610000093
分别计算出主影像和辅影像的大气延迟相位
Figure BDA0002843952610000094
和
Figure BDA0002843952610000095
将两者相减生成干涉SAR影像的大气延迟相位图
Figure BDA0002843952610000096
S60对主影像和辅影像进行干涉处理，生成差分干涉图，利用滤波算法对生成的干涉图进行滤波，并对其进行相位解缠；生成SAR解缠干涉相位图。S70获取大气改正的InSAR相位图：对干涉SAR影像的大气延迟相位图进行滤波处理，削弱差分延迟相位中的噪声以及操作误差，利用公式
Figure BDA0002843952610000097
剔除大气延迟效应，获得大气改正后干涉相位图，其中，
Figure BDA0002843952610000098
表示未进行改正的解缠干涉相位图，
Figure BDA0002843952610000099
表示大气延迟改正后的InSAR相位图。
所述步骤S10包括如下步骤：S11在研究区域设置多个GNSS接收机获取卫星发射的电磁波信号，所述电磁波信号包括伪距和载波相位观测数据，利用高精度GNSS数据处理软件GAMIT/GLOBK对所述观测数据进行解算，获得GNSS信号的位置信息以及GNSS测站处的气象参数值，其中所述位置信息包括高度角ele、方位角azi；所述气象参数值包括测站天顶对流层延迟ZTD、测站上空东西方向的大气梯度GEW和南北方向的大气梯度GNS。所述GNSS接收机的分布密度为15km2～20km2。
S12采用公式(1)计算GNSS测站处的天顶大气干延迟，其中，ZHD为天顶大气干延迟，PS表示GNSS测站处的地表气压，
Figure BDA0002843952610000101
表示GNSS测站的纬度，H表示GNSS测站的高程；
Figure BDA0002843952610000102
S13根据公式ZWD＝ZTD-ZHD计算出GNSS测站处的天顶大气湿延迟ZWD，同时将东西方向和南北方向的大气梯度分解为两个方向上的干延迟梯度
Figure BDA0002843952610000103
和湿延迟梯度
Figure BDA0002843952610000104
S14采用公式(2)分别计算出GNSS观测信号的斜路径干延迟和斜路径湿延迟，并将两者相加获得斜路径总延迟，STD＝SHD+SWD，STD表示斜路径总延迟。
Figure BDA0002843952610000105
其中，mfh和mfw分别表示干映射函数和湿映射函数，mfg表示水平梯度函数。
S15综合GNSS测站的分布情况以及SAR影像的观测范围确定GNSS层析模型的层析区域，将层析区域上空的三维空间进行离散化，在水平方向上按照一定的经度差和纬度差进行划分，在垂直方向上采用非均匀划分方式，形成p个均匀分布的体素块，构成GNSS三维体素块层析模型，其中，p＝n·m·t，n、m、t分别表示体素块的行数、列数和层数。优选所述层析区域的范围略大于所述SAR影像的观测范围，所述体素块在水平方向的边长为10～15km，所述体素块的总层数t为10～20层。
S16建立GNSS三维层析模型的观测方程组和约束方程组：根据每一条GNSS信号与三维体素块的空间位置关系可以构造GNSS层析观测方程
Figure BDA0002843952610000111
其中，STDq表示卫星信号线q的斜路径总延迟量，aijk表示该信号在第k层、第i行、第j列的体素块中的截距，Nijk表示该体素块的总折射率值；多条GNSS信号组成的层析观测方程可用矩阵形式表达STD＝A·Ν，其中STD表示由多条GNSS信号的斜路径总延迟组成的观测值向量，A表示多条信号对应的观测方程组的系数矩阵，N表示总折射率未知参数组成的列向量；同时，添加水平约束条件和垂直约束条件来解算层析观测方程组的病态性问题，所述水平约束条件选用基于高斯加权函数的约束模型，假设同一层面内某一体素块的总折射率与其周围体素块的总折射率存在式(3)的数量关系：
w1N1+w2N2+…+wi-1Ni-1-Ni+wi+1Ni+1+…+wnmNnm＝0 (3)
其中，N1，N2，…，Nnm表示同一层内nm个体素块的总折射率，w1，w2，…，wnm表示对应每个体素块的权系数，
Figure BDA0002843952610000112
dk表示同一层内第k个体素块与体素块i之间的距离，σ为平滑因子；
所述垂直约束条件选用基于指数递减函数的约束模型，不同高度层的体素块的总折射率之间存在式(4)所示的指数型关系：
Figure BDA0002843952610000113
其中，hj和hi分别表示第j层和第i层体素块的高度，N(hj)和N(hi)则表示对应体素块的总折射率，Hscale表示研究区域的水汽标高。
S17结合GNSS层析观测方程组、水平约束方程和垂直约束方程构建式(5)所示的GNSS三维层析模型，
Figure BDA0002843952610000121
其中，Aobs，AH，AV分别表示GNSS观测信号、水平约束、垂直约束对应的系数矩阵，Δobs，ΔH，ΔV分别表示对应的噪声，X为总折射率值。
S18通过代数重构算法(6)对公式(5)进行解算，获得GNSS层析的三维总折射率场，其中，yk和yk+1分别表示第k次和第k+1次迭代的近似解，λ为迭代松弛因子，bi为GNSS水汽层析模型左侧列向量的第i行元素，Ai表示层析模型系数矩阵A的第i行元素组成的行向量，
Figure BDA0002843952610000122
所述步骤S20中，选用研究区域内部的探空站的观测数据采用式(7)计算出所述探空数据在不同高度处的总折射率的参考值，利用层析解算的结果确定同一高度、同一位置处的总折射率值，然后分别根据式(8)、式(9)计算层析估值与探空参考值的平均偏差和均方根偏差；
Figure BDA0002843952610000123
Figure BDA0002843952610000124
Figure BDA0002843952610000125
其中，p为无线电探空获取的在某一高度的气压，e和T分别表示该高度处的水汽压和温度，Zi代表GNSS层析模型解算出来的体素块的总折射率，Zi'代表探空数据计算出来的体素块的总折射率值，s表示参与精度评定的体素块个数。当平均偏差小于0.5g/m3，均方根偏差小于1.5g/m3判定层析结果质量好，满足SAR信号大气延迟估计的要求。
如图3所示(a)图代表步骤S20中获取的初始总折射率场，在水平和垂直方向分别进行插值加密，形成图(b)所示的精细化三维总折射率场。所述步骤S30中，所述双线性插值算法如下式(10)所示：
Figure BDA0002843952610000131
其中，(b,l)为待插值点的坐标，Nwet(b,l)表示待插值点的总折射率值，(b11,l11)，(b21,l21)，(b12,l12)，(b22,l22)分别表示待插值点四周临近体素块的中心坐标，Nwet(b11,l11)，Nwet(b21,l21)，Nwet(b12,l12)，Nwet(b22,l22)则分别对应四周体素块的总折射率值；
所述三次样条插值算法如下式(11)所示：
Figure BDA0002843952610000132
其中，h表示待插值点的高度，Nwet(h)表示待插值点的总折射率值，hi，hi-1分别表示插值点所在高度层的上下两个层面的高度，Δhi为该层的层间高度，即hi-hi-1，Nwet(hi)，Nwet(hi-1)为对应两个高度的总折射率值，Nwet(h1)，Nwet(h2)，…，Nwet(hk)则分别表示层析范围内每一高度层处的总折射率值，Di,:，Di-1,:分别表示转换矩阵D的第i列，第i-1列对应的列向量，式(11)中，D＝-A-1B表示转换矩阵，只与垂直分层有关，更详细地，其中A，B的计算公式如式(12)及式(13)所示
Figure BDA0002843952610000141
Figure BDA0002843952610000142
如图4所示，展示了GNSS信号(虚线)和SAR信号(实线)在精细化三维总折射率场的二维示意图，根据SAR信号在每个精细化体素块的距离信息和该体素块的总折射率值，采用Ray-tracing算法计算SAR信号的大气延迟量。所述步骤S40中，所述Ray-tracing算法如下式(14)所示：
Figure BDA0002843952610000143
其中，n＝1+N×10-6表示大气总折射率指数，si为SAR信号穿过精细化三维总折射率场中第i层的距离。
以上所述仅为本发明的示例性实施例，并非因此限制本发明专利保护范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。

Claims (8)

1.一种GNSS层析技术辅助的InSAR大气延迟改正方法，其特征在于，包括如下步骤：
S10根据研究区域内GNSS信号的斜路径总延迟观测量，采用GNSS层析技术重构研究区域内大气水汽的三维总折射率场；
S20以研究区域的无线电探空数据为参考值，对层析结果进行精度评定，若精度符合要求，进行步骤S30，否则，重复步骤S10，优化GNSS层析模型的配置与解算策略，重新反演三维总折射率场；
S30利用空间插值算法获取精细化三维总折射率场：在对流层上下两部分以不同的空间分辨率对初始三维总折射率场进行插值加密，在水平方向采用双线性插值算法，在垂直方向采用三次样条插值算法计算出加密后每个小体素块的总折射率值，构造精细化三维总折射率场；
S40利用Ray-tracing算法估计主影像和辅影像中SAR信号的斜路径总延迟量：根据SAR卫星的轨道信息与SAR影像中每个像元的经纬度坐标信息可构造出多条SAR卫星观测信号，每一个像元与SAR卫星之间的连线相当于一条观测信号，针对SAR卫星观测信号，利用Ray-tracing算法估计出每一条三维SAR信号的斜路径总延迟量；
S50将步骤S40计算的SAR主影像和辅影像中SAR信号的斜路径总延迟记为
Figure FDA0002843952600000011
和
Figure FDA0002843952600000012
利用公式
Figure FDA0002843952600000013
分别计算出主影像和辅影像的大气延迟相位
Figure FDA0002843952600000014
和
Figure FDA0002843952600000015
将两者相减生成干涉SAR影像的大气延迟相位图
Figure FDA0002843952600000016
S60对主影像和辅影像进行干涉处理，生成差分干涉图，利用滤波算法对生成的干涉图进行滤波，并对其进行相位解缠；生成SAR解缠干涉相位图；
S70获取大气改正的InSAR相位图：对干涉SAR影像的大气延迟相位图进行滤波处理，削弱差分延迟相位中的噪声以及操作误差，利用公式
Figure FDA0002843952600000021
剔除大气延迟效应，获得大气改正后干涉相位图，其中，
Figure FDA0002843952600000022
表示未进行改正的解缠干涉相位图，
Figure FDA0002843952600000023
表示大气延迟改正后的InSAR相位图。
2.根据权利要求1所述的GNSS层析技术辅助的InSAR大气延迟改正方法，
其特征在于，所述步骤S10包括如下步骤：
S11在研究区域设置多个GNSS接收机获取卫星发射的电磁波信号，所述电磁波信号包括伪距和载波相位观测数据，利用高精度GNSS数据处理软件GAMIT/GLOBK对所述观测数据进行解算，获得GNSS信号的位置信息以及GNSS测站处的气象参数值，其中所述位置信息包括高度角ele、方位角azi；所述气象参数值包括测站天顶对流层延迟ZTD、测站上空东西方向的大气梯度GEW和南北方向的大气梯度GNS；
S12采用公式(1)计算GNSS测站处的天顶大气干延迟，其中，ZHD为天顶大气干延迟，PS表示GNSS测站处的地表气压，
Figure FDA0002843952600000024
表示GNSS测站的纬度，H表示GNSS测站的高程；
Figure FDA0002843952600000025
S13根据公式ZWD＝ZTD-ZHD计算出GNSS测站处的天顶大气湿延迟ZWD，同时将东西方向和南北方向的大气梯度分解为两个方向上的干延迟梯度
Figure FDA0002843952600000026
和湿延迟梯度
Figure FDA0002843952600000027
S14采用公式(2)分别计算出GNSS观测信号的斜路径干延迟和斜路径湿延迟，并将两者相加获得斜路径总延迟，STD＝SHD+SWD，STD表示斜路径总延迟；
Figure FDA0002843952600000031
其中，mfh和mfw分别表示干映射函数和湿映射函数，mfg表示水平梯度函数；
S15综合GNSS测站的分布情况以及SAR影像的观测范围确定GNSS层析模型的层析区域，将层析区域上空的三维空间进行离散化，在水平方向上按照一定的经度差和纬度差进行划分，在垂直方向上采用非均匀划分方式，形成p个均匀分布的体素块，构成GNSS三维体素块层析模型，其中，p＝n·m·t，n、m、t分别表示体素块的行数、列数和层数；
S16建立GNSS三维层析模型的观测方程组和约束方程组：根据每一条GNSS信号与三维体素块的空间位置关系可以构造GNSS层析观测方程
Figure FDA0002843952600000032
其中，STDq表示卫星信号线q的斜路径总延迟量，aijk表示该信号在第k层、第i行、第j列的体素块中的截距，Nijk表示该体素块的总折射率值；多条GNSS信号组成的层析观测方程可用矩阵形式表达STD＝A·Ν，其中STD表示由多条GNSS信号的斜路径总延迟组成的观测值向量，A表示多条信号对应的观测方程组的系数矩阵，N表示总折射率未知参数组成的列向量；同时，添加水平约束条件和垂直约束条件来解算层析观测方程组的病态性问题，所述水平约束条件选用基于高斯加权函数的约束模型，假设同一层面内某一体素块的总折射率与其周围体素块的总折射率存在式(3)的数量关系：
w1N1+w2N2+…+wi-1Ni-1-Ni+wi+1Ni+1+…+wnmNnm＝0 (3)
其中，N1，N2，…，Nnm表示同一层内nm个体素块的总折射率，w1，w2，…，wnm表示对应每个体素块的权系数，
Figure FDA0002843952600000033
dk表示同一层内第k个体素块与体素块i之间的距离，σ为平滑因子；
所述垂直约束条件选用基于指数递减函数的约束模型，不同高度层的体素块的总折射率之间存在式(4)所示的指数型关系：
Figure FDA0002843952600000041
其中，hj和hi分别表示第j层和第i层体素块的高度，N(hj)和N(hi)则表示对应体素块的总折射率，Hscale表示研究区域的水汽标高。
S17结合GNSS层析观测方程组、水平约束方程和垂直约束方程构建式(5)所示的GNSS三维层析模型，
Figure FDA0002843952600000042
其中，Aobs，AH，AV分别表示GNSS观测信号、水平约束、垂直约束对应的系数矩阵，Δobs，ΔH，ΔV分别表示对应的噪声，X为总折射率值；
S18通过代数重构算法(6)对公式(5)进行解算，获得GNSS层析的三维总折射率场，其中，yk和yk+1分别表示第k次和第k+1次迭代的近似解，λ为迭代松弛因子，bi为GNSS水汽层析模型左侧列向量的第i行元素，Ai表示层析模型系数矩阵A的第i行元素组成的行向量，
Figure FDA0002843952600000043
3.根据权利要求2所述的GNSS层析技术辅助的InSAR大气延迟改正方法，其特征在于，所述GNSS接收机的分布密度为15km2～20km2。
4.根据权利要求3所述的GNSS层析技术辅助的InSAR大气延迟改正方法，其特征在于，所述层析区域的范围略大于所述SAR影像的观测范围，所述体素块在水平方向的边长为10～15km，所述体素块的总层数t为10～20层。
5.根据权利要求4所述的GNSS层析技术辅助的InSAR大气延迟改正方法，其特征在于，所述步骤S20中，选用研究区域内部的探空站的观测数据采用式(7)计算出所述探空数据在不同高度处的总折射率的参考值，利用层析解算的结果确定同一高度、同一位置处的总折射率值，然后分别根据式(8)、式(9)计算层析估值与探空参考值的平均偏差和均方根偏差；
Figure FDA0002843952600000051
Figure FDA0002843952600000052
Figure FDA0002843952600000053
其中，p为无线电探空获取的在某一高度的气压，e和T分别表示该高度处的水汽压和温度，Zi代表GNSS层析模型解算出来的体素块的总折射率，Zi'代表探空数据计算出来的体素块的总折射率值，s表示参与精度评定的体素块个数。
6.根据权利要求5所述的GNSS层析技术辅助的InSAR大气延迟改正方法，其特征在于，当平均偏差小于0.5g/m3，均方根偏差小于1.5g/m3判定层析结果质量好，满足SAR信号大气延迟估计的要求。
7.根据权利要求5所述的GNSS层析技术辅助的InSAR大气延迟改正方法，其特征在于，所述步骤S30中，
所述双线性插值算法如下式(10)所示：
Figure FDA0002843952600000054
其中，(b,l)为待插值点的坐标，Nwet(b,l)表示待插值点的总折射率值，(b11,l11)，(b21,l21)，(b12,l12)，(b22,l22)分别表示待插值点四周临近体素块的中心坐标，Nwet(b11,l11)，Nwet(b21,l21)，Nwet(b12,l12)，Nwet(b22,l22)则分别对应四周体素块的总折射率值；
所述三次样条插值算法如下式(11)所示：
Figure FDA0002843952600000061
其中，h表示待插值点的高度，Nwet(h)表示待插值点的总折射率值，hi，hi-1分别表示插值点所在高度层的上下两个层面的高度，Δhi为该层的层间高度，即hi-hi-1，Nwet(hi)，Nwet(hi-1)为对应两个高度的总折射率值，Nwet(h1)，Nwet(h2)，…，Nwet(hk)则分别表示层析范围内每一高度层处的总折射率值，Di,:，Di-1,:分别表示转换矩阵D的第i列，第i-1列对应的列向量，式(11)中，D＝-A-1B表示转换矩阵，只与垂直分层有关，更详细地，其中A，B的计算公式如式(12)及式(13)所示
Figure FDA0002843952600000062
Figure FDA0002843952600000071
8.根据权利要求1所述的GNSS层析技术辅助的InSAR大气延迟改正方法，其特征在于，所述步骤S40中，所述Ray-tracing算法如下式(14)所示：
Figure FDA0002843952600000072
其中，n＝1+N×10-6表示大气总折射率指数，si为SAR信号穿过精细化三维总折射率场中第i层的距离。
CN202011502257.9A 2020-12-18 2020-12-18 一种GNSS层析技术辅助的InSAR大气延迟改正方法 Active CN112711022B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011502257.9A CN112711022B (zh) 2020-12-18 2020-12-18 一种GNSS层析技术辅助的InSAR大气延迟改正方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011502257.9A CN112711022B (zh) 2020-12-18 2020-12-18 一种GNSS层析技术辅助的InSAR大气延迟改正方法

Publications (2)

Publication Number Publication Date
CN112711022A true CN112711022A (zh) 2021-04-27
CN112711022B CN112711022B (zh) 2022-08-30

Family

ID=75544450

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011502257.9A Active CN112711022B (zh) 2020-12-18 2020-12-18 一种GNSS层析技术辅助的InSAR大气延迟改正方法

Country Status (1)

Country Link
CN (1) CN112711022B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113534213A (zh) * 2021-07-26 2021-10-22 中国电子科技集团公司第五十四研究所 一种公里级区域大气相位不一致性高精度建模与修正方法
CN116340710A (zh) * 2023-05-30 2023-06-27 中国科学院精密测量科学与技术创新研究院 基于分层快速三维射线追踪的中性大气斜延迟计算方法

Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7446705B1 (en) * 2007-10-24 2008-11-04 Wisconsin Alumni Research Foundation Method and apparatus for determining parameters for a parametric expression characterizing the phase of an acquired signal
US20120319893A1 (en) * 2011-06-20 2012-12-20 California Institute Of Technology Damage proxy map from interferometric synthetic aperture radar coherence
CN103885046A (zh) * 2012-12-20 2014-06-25 河南省电力勘测设计院 基于GPS的InSAR大气延迟改正方法
CN105842692A (zh) * 2016-03-17 2016-08-10 中国科学院遥感与数字地球研究所 一种insar测量中的大气校正方法
CN109782282A (zh) * 2019-03-13 2019-05-21 武汉大学 一种集成对流层大气延迟改正的时间序列InSAR分析方法
CN110031841A (zh) * 2019-04-01 2019-07-19 中国科学院遥感与数字地球研究所 基于ECMWF的InSAR大气延迟改正的方法和系统
CN110058236A (zh) * 2019-05-21 2019-07-26 中南大学 一种面向三维地表形变估计的InSAR和GNSS定权方法
US20200258296A1 (en) * 2019-02-08 2020-08-13 Ursa Space Systems Inc. Satellite sar artifact suppression for enhanced three-dimensional feature extraction, change detection, and visualizations
CN111998766A (zh) * 2020-08-31 2020-11-27 同济大学 一种基于时序InSAR技术的地表形变反演方法
CN112050725A (zh) * 2020-09-14 2020-12-08 广东省核工业地质局测绘院 一种融合InSAR和GPS的地表形变监测方法

Patent Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7446705B1 (en) * 2007-10-24 2008-11-04 Wisconsin Alumni Research Foundation Method and apparatus for determining parameters for a parametric expression characterizing the phase of an acquired signal
US20120319893A1 (en) * 2011-06-20 2012-12-20 California Institute Of Technology Damage proxy map from interferometric synthetic aperture radar coherence
CN103885046A (zh) * 2012-12-20 2014-06-25 河南省电力勘测设计院 基于GPS的InSAR大气延迟改正方法
CN105842692A (zh) * 2016-03-17 2016-08-10 中国科学院遥感与数字地球研究所 一种insar测量中的大气校正方法
US20200258296A1 (en) * 2019-02-08 2020-08-13 Ursa Space Systems Inc. Satellite sar artifact suppression for enhanced three-dimensional feature extraction, change detection, and visualizations
CN109782282A (zh) * 2019-03-13 2019-05-21 武汉大学 一种集成对流层大气延迟改正的时间序列InSAR分析方法
CN110031841A (zh) * 2019-04-01 2019-07-19 中国科学院遥感与数字地球研究所 基于ECMWF的InSAR大气延迟改正的方法和系统
CN110058236A (zh) * 2019-05-21 2019-07-26 中南大学 一种面向三维地表形变估计的InSAR和GNSS定权方法
CN111998766A (zh) * 2020-08-31 2020-11-27 同济大学 一种基于时序InSAR技术的地表形变反演方法
CN112050725A (zh) * 2020-09-14 2020-12-08 广东省核工业地质局测绘院 一种融合InSAR和GPS的地表形变监测方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
MARION HEUBLEIN ET AL.: "COMPRESSIVE SENSING FOR NEUTROSPHERIC WATER VAPOR TOMOGRAPHY USING GNSS AND INSAR OBSERVATIONS", 《IEEE》 *
MARION HEUBLEIN ET AL.: "Compressive sensing reconstruction of 3D wet refractivity based on GNSS and InSAR observations", 《JOURNAL OF GEODESY》 *
WENYUAN ZHANG ET AL.: "A Tropospheric Tomography Method with a Novel Height Factor Model Including Two Parts:Isotropic and Anisotropic Height Factors", 《REMOTE SENSING 》 *
WENYUAN ZHANG ET AL.: "AN IMPROVED TROPOSPHERIC TOMOGRAPHY METHOD BASED ON THE DYNAMIC NODE PARAMETERIZED ALGORITHM", 《ACTA GEODYNAMICA ET GEOMATERLALIA》 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113534213A (zh) * 2021-07-26 2021-10-22 中国电子科技集团公司第五十四研究所 一种公里级区域大气相位不一致性高精度建模与修正方法
CN113534213B (zh) * 2021-07-26 2022-06-10 中国电子科技集团公司第五十四研究所 一种公里级区域大气相位不一致性高精度建模与修正方法
CN116340710A (zh) * 2023-05-30 2023-06-27 中国科学院精密测量科学与技术创新研究院 基于分层快速三维射线追踪的中性大气斜延迟计算方法
CN116340710B (zh) * 2023-05-30 2023-09-12 中国科学院精密测量科学与技术创新研究院 基于分层快速三维射线追踪的中性大气斜延迟计算方法

Also Published As

Publication number Publication date
CN112711022B (zh) 2022-08-30

Similar Documents

Publication Publication Date Title
CN109782282B (zh) 一种集成对流层大气延迟改正的时间序列InSAR分析方法
Mateus et al. Sentinel-1 interferometric SAR mapping of precipitable water vapor over a country-spanning area
Benevides et al. Bridging InSAR and GPS tomography: a new differential geometrical constraint
Mateus et al. Can spaceborne SAR interferometry be used to study the temporal evolution of PWV?
CN112711022B (zh) 一种GNSS层析技术辅助的InSAR大气延迟改正方法
Abdelfattah et al. Topographic SAR interferometry formulation for high-precision DEM generation
CN113671505B (zh) 一种基于系统几何误差补偿的合成孔径雷达立体定位方法
CN114689015B (zh) 一种提高光学卫星立体影像dsm高程精度的方法
Xia et al. Assessing water vapor tomography in Hong Kong with improved vertical and horizontal constraints
Tang et al. Atmospheric correction in time-series SAR interferometry for land surface deformation mapping–A case study of Taiyuan, China
Zhang et al. GNSS-RS tomography: Retrieval of tropospheric water vapor fields using GNSS and RS observations
Tang et al. High-spatial-resolution mapping of precipitable water vapour using SAR interferograms, GPS observations and ERA-Interim reanalysis
CN114357770A (zh) 一种对流层层析方法
Zhang et al. A new method for tropospheric tomography using GNSS and Fengyun-4A data
Feng et al. An improved geometric calibration model for spaceborne SAR systems with a case study of large-scale Gaofen-3 images
CN113009531A (zh) 一种小尺度高精度的低空对流层大气折射率模型
CN115980317B (zh) 基于修正相位的地基gnss-r数据土壤水分估算方法
CN115980751A (zh) 一种幂律模型InSAR对流层延迟改正方法
Izanlou et al. Enhanced Troposphere Tomography: Integration of GNSS and Remote Sensing Data with Optimal Vertical Constraints
CN111126466A (zh) 一种多源pwv数据融合方法
CN116466376A (zh) 数值预报模式辅助的实时ppp改善方法
CN115755115A (zh) 基于gnss对流层层析技术的ppp改善方法
CN115755103A (zh) 一种抗差自适应的gnss水汽层析方法
CN113534213B (zh) 一种公里级区域大气相位不一致性高精度建模与修正方法
CN113341410B (zh) 一种大范围林下地形估计方法、装置、设备及介质

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