CN116485857A - 一种基于多源遥感数据的高时间分辨率冰川厚度反演方法 - Google Patents
一种基于多源遥感数据的高时间分辨率冰川厚度反演方法 Download PDFInfo
- Publication number
- CN116485857A CN116485857A CN202310497009.7A CN202310497009A CN116485857A CN 116485857 A CN116485857 A CN 116485857A CN 202310497009 A CN202310497009 A CN 202310497009A CN 116485857 A CN116485857 A CN 116485857A
- Authority
- CN
- China
- Prior art keywords
- glacier
- image pair
- resolution
- remote sensing
- flow velocity
- 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.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 75
- 238000006073 displacement reaction Methods 0.000 claims abstract description 44
- 230000003287 optical effect Effects 0.000 claims abstract description 34
- 238000009499 grossing Methods 0.000 claims abstract description 8
- 238000001914 filtration Methods 0.000 claims abstract description 5
- 239000011159 matrix material Substances 0.000 claims description 24
- 230000008859 change Effects 0.000 claims description 14
- 238000004364 calculation method Methods 0.000 claims description 11
- 238000004590 computer program Methods 0.000 claims description 10
- 238000004422 calculation algorithm Methods 0.000 claims description 9
- 230000006870 function Effects 0.000 claims description 7
- 230000033001 locomotion Effects 0.000 claims description 5
- 230000008569 process Effects 0.000 claims description 5
- 230000000007 visual effect Effects 0.000 claims description 4
- 208000019901 Anxiety disease Diseases 0.000 claims description 3
- 238000012952 Resampling Methods 0.000 claims description 3
- 238000010276 construction Methods 0.000 claims description 3
- 238000012545 processing Methods 0.000 claims description 3
- 239000002245 particle Substances 0.000 claims description 2
- 230000000694 effects Effects 0.000 description 4
- 230000002159 abnormal effect Effects 0.000 description 3
- 238000012937 correction Methods 0.000 description 3
- 238000012544 monitoring process Methods 0.000 description 3
- 230000007812 deficiency Effects 0.000 description 2
- 238000013461 design Methods 0.000 description 2
- 239000000463 material Substances 0.000 description 2
- 230000009466 transformation Effects 0.000 description 2
- 238000002679 ablation Methods 0.000 description 1
- 230000002411 adverse Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 230000004927 fusion Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000002093 peripheral effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 238000011144 upstream manufacturing Methods 0.000 description 1
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/86—Combinations of radar systems with non-radar systems, e.g. sonar, direction finder
- G01S13/867—Combination of radar systems with cameras
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
- G01B21/00—Measuring arrangements or details thereof, where the measuring technique is not covered by the other groups of this subclass, unspecified or not relevant
- G01B21/02—Measuring arrangements or details thereof, where the measuring technique is not covered by the other groups of this subclass, unspecified or not relevant for measuring length, width, or thickness
- G01B21/08—Measuring arrangements or details thereof, where the measuring technique is not covered by the other groups of this subclass, unspecified or not relevant for measuring length, width, or thickness for measuring thickness
-
- 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
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/05—Geographic models
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/20—Image enhancement or restoration using local operators
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/50—Image enhancement or restoration using two or more images, e.g. averaging or subtraction
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/70—Denoising; Smoothing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
- G06T7/33—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
- G06T7/337—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods involving reference images or patches
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20212—Image combination
- G06T2207/20221—Image fusion; Image merging
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Remote Sensing (AREA)
- Radar, Positioning & Navigation (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Geometry (AREA)
- Computer Networks & Wireless Communication (AREA)
- Software Systems (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Electromagnetism (AREA)
- Computer Graphics (AREA)
- Image Analysis (AREA)
Abstract
本发明公开了一种基于多源遥感数据的高时间分辨率冰川厚度反演方法，如下：融合第三SAR影像对、光学影像对的各方向偏移量、与冰川三维流速矢量建立冰川表面三维流速方程；求解冰川表面三维流速方程，通过引入L1范数最小法来探测粗差并确定第一次平差的残差值；对残差值进行赋权得到权值；结合权值利用最小二乘准则得到三维流速的最佳估计值；三维流速的最佳估计值包括东西向、南北和垂直向；将由垂直向位移量求取的垂直速率分解为表面平行流和非表面平行流分量；结合东西向位移量、南北向位移量计算得到的水平速度矢量，与某一高度空间分辨率DEM数据进行高斯平滑滤波得到的表面坡度，计算表面平行流和非表面平行流分量。
Description
技术领域
本发明涉及SAR数据处理技术领域，更具体的，涉及一种基于多源遥感数据的高时间分辨率冰川厚度反演方法。
背景技术
随着当代卫星重返频率和空间分辨率的逐步提高，为冰川物质平衡研究提供了海量的观测数据，可以实现大范围、高空间分辨率、连续动态观测，同时也给高精度的冰川厚度反演带来了新的挑战。对于冰川厚度监测主要采用dem差分的方法，时间分辨率低，滞后性严重。这使得在实现近实时冰厚监测要求和时间效率上难以兼顾。目前，时序冰川厚度已进入基于多源平台的高时频监测阶段，其核心思想是对各个不同平台的数据进行偏移量估计，构建设计矩阵后先利用一范数最小法再通过最小二乘求解冰川的三维流速，最后将垂直流速进一步分解得到冰川厚度变化。
基于多源平台数据求解冰川三维流速以此探究冰川厚度变化仍是待攻克的难题，传统的基于小基线集技术求解方程组存在秩亏的问题。同时，不同观测几何或不同卫星传感器的观测结果精度不尽相同，不可避免有粗差存在于观测值中。而现有的三维流速求解中没有考虑到粗差对解算的影响，不具备抵抗粗差干扰的能力。当观测中存在粗差时，会将其无差别等概率分配到其他观测值中，导致无法从残差结果进行粗差观测值的定位，进一步损害三维流速的平差结果。因此，在进行三维流速反演中，必须同时考虑方程的设计和观测数据中存在异常值干扰的问题。
发明内容
本发明为了解决针对SAR时序解算方程组秩亏和粗差对三维流速估计造成不良影响的问题，提供了一种基于多源遥感数据的高时间分辨率冰川厚度反演方法，其采用L1范数最小法对含有异常误差的观测值进行探测，再利用最小二乘法适当降低在平差中的权重，使权重与其实际精度相适应，降低了三维求解结果受粗差干扰的影响，极大程度的保证了三维估计结果的可靠性。
为实现上述本发明目的，采用的技术方案如下：
一种基于多源遥感数据的高时间分辨率冰川厚度反演方法，所述的方法包括步骤如下：
融合第三SAR影像对、光学影像对的各方向偏移量、与冰川三维流速矢量建立冰川表面三维流速方程；所述的第三SAR影像与光学影像对的偏移量结果分辨率一致；
求解冰川表面三维流速方程，具体通过引入L1范数最小法来探测粗差并确定第一次平差的残差值V；对得到的残差值V进行赋权得到权值P；结合权值P并利用最小二乘准则得到三维地表形变的最佳估计值；
所述的三维地表形变的最佳估计值包括东西向位移量、南北向位移量、垂直向位移量；
将由垂直向位移量求取的垂直速率分解为由沿着冰川表面斜坡运动引起的表面平行流和由于冰川表面高程变化引起的非表面平行流分量；
结合东西向位移量、南北向位移量计算得到的水平速度矢量，与某一高度空间分辨率DEM数据进行高斯平滑滤波得到的表面坡度，计算表面平行流和非表面平行流分量。
优选地，所述的第三SAR影像对通过以下处理得到：
获取第一平台遥感数据和第二平台遥感数据组成第一SAR影像对，通过设置时间基线阈值、空间基线阈值选择得到最优的第一SAR影像对；
根据某一高度空间分辨率DEM数据对最优的第一SAR影像对进行精确配准，得到第二SAR影像对；
采用强度跟踪算法对第二SAR影像对进行偏移量的估算，每对SAR影像得到视线向、方位向的第一偏移量；对第一偏移量进行地理编码并依据最邻近法重采样至网格单元中，得到与光学影像对的偏移量结果分辨率一致的第三SAR影像对。
优选地，在将光学影像对与第三SAR影像对融合之前，先采用强度跟踪算法对光学影像对进行偏移量的估算，剔除误差，得到光学影像对的东西向、南北向的第二偏移量。
融合第三SAR影像对与光学影像对的各方向偏移量，与冰川三维流速矢量的东西方向、北南方向和垂直方向建立冰川表面三维流速方程如下：
其中：
式中，A和D分别代表第一平台遥感数据的升轨数据和降轨数据，代表第二平台遥感数据的升轨数据；/>代表光学影像；/>是SAR卫星的方位角，即卫星飞行方向与正北方向沿顺时针的夹角；/>是卫星的入射角；/>和/>分别表示由SAR影像对得到的地面点的视线向位移量和方位向的第一偏移量；/>和/>分别表示由光学影像对得到的地面点的东西向位移量和南北向的第二偏移量；/>是在/>时间间隔内计算的偏移量；B为系数矩阵，X为观测点三维流速，L为观测量；/>表示冰川三维流速矢量的东西方向，/>表示冰川三维流速矢量的北南方向，/>表示冰川三维流速矢量的垂直方向。
进一步地，通过引入L1范数最小法来探测粗差并确定第一次平差的残差值V，公式如下：
其中，V为8x1残差矩阵；B为已知的8x3系数矩阵；X为3x1未知参数向量；L为8x1已知观测矩阵，观测量的数量。
再进一步地，所述的权值P为权矩阵，采用表示权矩阵中的某个元素，/>为权因子，/>表示对角矩阵；
其中，通过Bisquare权函数构造得到：
。
再进一步地，结合权值P并利用最小二乘准则得到冰川三维流速矢量的最佳估计值，法方程公式如下：
其中：
冰川三维流速矢量的最佳估计值可由以下公式计算得到：
其中，、/>和/>分别代表东西方向、北南方向和垂直方向的流速最佳估计值。
再进一步地，计算表面平行流和非表面平行流分量，具体如下：
引入冰川运动学边界条件与冰川表面高程变化率，将计算公式如下：
等价于：
其中，为垂直冰粒速度；/>为水平速度矢量；/>表示某一高度空间分辨率DEM数据进行高斯平滑滤波得到的表面坡度，/>，/>表示高程，/>表示沿东西方向，/>表示沿北南方向；/>表示质量平衡率；
根据Sorge恒定密度定律，忽略速度的垂直剖面，将冰川厚度变化转换为等效质量，将由垂直向位移量求取的垂直速率分解为表面平行流和非表面平行流分量；
其中，表面平行流的计算公式为，非表面平行流分量的计算公式为。
一种计算机设备，包括存储器和处理器，所述存储器存储有可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现如上所述方法的步骤。
一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现如上所述方法的步骤。
本发明的有益效果如下：
本发明通过引入多源遥感数据的不同方向上偏移量，与冰川三维流速矢量建立冰川表面三维流速方程进行三维流速反演，使得冰川表面三维流速方程的自由度提高的同时，在冰川表面三维流速方程求解中，引入L1范数最小法探测粗差；首先采用L1范数最小法对含有异常误差的观测值进行探测，再利用最小二乘法适当降低在平差中的权重，使权重与其实际精度相适应，降低了三维求解结果受粗差干扰的影响，极大程度的保证了三维估计结果的可靠性，从而提高求解参数的精度，再以此将由垂直向位移量求取的垂直速率分解为表面平行流和非表面平行流分量，得到可靠的冰川厚度变化。
附图说明
图1是实施例1所述的基于多源遥感数据的高时间分辨率冰川厚度反演方法的步骤流程图。
图2是实施例2所述的基于多源遥感数据的高时间分辨率冰川厚度反演方法的步骤流程图。
图3是三维误差椭球图，描述了使用不同数据集进行三维速度反演的不确定性；其中(a)仅使用S1数据反演；(b)使用S1和S2数据联合反演；(c) S1、S2和GF-3联合反演。
具体实施方式
下面结合附图和具体实施方式对本发明做详细描述。
实施例1
如图1所示，一种基于多源遥感数据的高时间分辨率冰川厚度反演方法，所述的方法包括步骤如下：
融合第三SAR影像对、光学影像对的各方向偏移量、与冰川三维流速矢量建立冰川表面三维流速方程；所述的第三SAR影像与光学影像对的偏移量结果分辨率一致；
求解冰川表面三维流速方程，具体通过引入L1范数最小法来探测粗差并确定第一次平差的残差值V；对得到的残差值V进行赋权得到权值P；结合权值P并利用最小二乘准则得到三维地表形变的最佳估计值；
所述的三维地表形变的最佳估计值包括东西向位移量、南北向位移量、垂直向位移量；
将由垂直向位移量求取的垂直速率分解为由沿着冰川表面斜坡运动引起的表面平行流和由于冰川表面高程变化引起的非表面平行流分量；
结合东西向位移量、南北向位移量计算得到的水平速度矢量，与某一高度空间分辨率DEM数据进行高斯平滑滤波得到的表面坡度，计算表面平行流和非表面平行流分量。
本实施例通过引入多源遥感数据的不同方向上偏移量，与冰川三维流速矢量建立冰川表面三维流速方程进行三维流速反演，使得冰川表面三维流速方程的自由度提高的同时，在冰川表面三维流速方程求解中，引入L1范数最小法探测粗差；首先采用L1范数最小法对含有异常误差的观测值进行探测，再利用最小二乘法适当降低在平差中的权重，使权重与其实际精度相适应，降低了三维求解结果受粗差干扰的影响，极大程度的保证了三维估计结果的可靠性，从而提高求解参数的精度，再以此将由垂直向位移量求取的垂直速率分解为表面平行流和非表面平行流分量，得到可靠的冰川厚度变化。
实施例2
基于实施例1所述的基于多源遥感数据的高时间分辨率冰川厚度反演方法，更具体地，如图2所示，本实施例提供的基于多源遥感数据的高时间分辨率冰川厚度反演方法，具体如下：
基于多源遥感数据的高时间分辨率冰川厚度反演方法，所述的方法包括步骤如下：
融合第三SAR影像对、光学影像对的各方向偏移量，与冰川三维流速矢量建立冰川表面三维流速方程；所述的第三SAR影像与光学影像对的偏移量结果分辨率一致；
所述的第三SAR影像对通过以下处理得到：
获取第一平台遥感数据和第二平台遥感数据组成第一SAR影像对，通过设置时间基线阈值、空间基线阈值选择得到最优的第一SAR影像对；
根据某一高度空间分辨率DEM数据对最优的第一SAR影像对进行精确配准，得到第二SAR影像对；
采用强度跟踪算法对第二SAR影像对进行偏移量的估算，每对SAR影像得到视线向（也即距离向）、方位向的第一偏移量；对第一偏移量进行地理编码并依据最邻近法重采样至网格单元中，得到与光学影像对的偏移量结果分辨率一致的第三SAR影像对。
本实施例可以采用哨兵1号卫星（Sentinel-1）的升降轨数据作为第一平台遥感数据，和高分三号卫星（GF-3）的升轨数据作为第二平台遥感数据组成一个主辅的第一SAR影像对。
本实施例在采用强度跟踪算法对第一SAR影像对进行偏移量的估算时，相关系数阈值设置为0.2，以此掩膜在失相干严重区域获得的不可靠偏移量。
求解冰川表面三维流速方程，具体通过引入L1范数最小法来探测粗差并确定第一次平差的残差值V；对得到的残差值V进行赋权得到权值P；结合权值P并利用最小二乘准则得到三维地表形变的最佳估计值；
所述的三维地表形变的最佳估计值包括东西向位移量、南北向位移量、垂直向位移量；
将由垂直向位移量求取的垂直速率分解为由沿着冰川表面斜坡运动引起的表面平行流和由于冰川表面高程变化引起的非表面平行流分量；
结合东西向位移量、南北向位移量计算得到的水平速度矢量，与某一高度空间分辨率DEM数据进行高斯平滑滤波得到的表面坡度，计算表面平行流和非表面平行流分量。
在一个具体的实施例中，在将光学影像对与第三SAR影像对融合之前，先采用强度跟踪算法对光学影像对进行偏移量的估算，剔除误差，得到光学影像对的东西向、南北向的第二偏移量。
在一个具体的实施例中，基于傅里叶偏移算法对第三平台遥感数据的某一波段进行偏移量计算；其依据是傅里叶偏移理论，进行傅里叶逆变换可以得到一个冲击函数，该函数在其他位置都为零，只在偏移位置处不为零，该位置即为水平地表形变量。在本实施例所述的第三平台遥感数据可以采用哨兵2号卫星（Sentinel-2）的遥感数据，具体基于傅里叶偏移算法对哨兵2号卫星遥感数据的第八波段进行偏移量计算。
在一个具体的实施例中，设和/>是两幅图像，/>是两幅图像的偏移量，则有：
进行傅里叶变化之后为：
式中，分别代表行、列的频率变化；
设是/>和/>的互能量谱，则有：
式中，*代表复数共轭相乘；
基于上述公式进行傅里叶逆变换得到一个冲击函数，冲击函数在其他位置都为零，只在偏移位置处不为零，该位置即为水平地表形变量。
在一个具体的实施例中，设方位向（也即沿卫星轨道方向）为正，视线向（也即远离雷达卫星方向）为正，垂直向上、正东、正北方向为正；
本实施例还在建立冰川表面三维流速方程之前，分别对第三SAR影像对、光学影像对均进行正北方向校正、基准校正。
其中，真北方向校正具体如下：
GAMMA软件可以将从Sentinel-1和GF-3追踪到的距离向和方位向变形转到地理坐标系下。光学变形也需从格网坐标系转化为地理东西方向和地理南北方向下的值。对于每个像元，计算子午线收敛角，和子午线收敛，然后生成旋转矩阵，用其转换坐标变换过程基于下式：
其中，表示子午线收敛角，/>表示像元纬度，/>表示像元经度与中央子午线之差，和/>表示网格坐标系中东西和南北方向的速度，/>和/>表示地理坐标系中东西和南北方向的偏移量。
融合第三SAR影像对与光学影像对的各方向偏移量，与冰川三维流速矢量的东西方向、北南方向和垂直方向建立冰川表面三维流速方程如下：
其中：
式中，A和D分别代表第一平台遥感数据的升轨数据和降轨数据，代表第二平台遥感数据的升轨数据；/>代表光学影像；/>是SAR卫星的方位角，即卫星飞行方向与正北方向沿顺时针的夹角；/>是卫星的入射角；/>和/>分别表示由SAR影像对得到的地面点的视线向位移量和方位向的第一偏移量；/>和/>分别表示由光学影像对得到的地面点的东西向位移量和南北向的第二偏移量；/>是在/>时间间隔内计算的偏移量；B为系数矩阵，X为观测点三维流速，L为观测量；/>表示冰川三维流速矢量的东西方向，/>表示冰川三维流速矢量的北南方向，/>表示冰川三维流速矢量的垂直方向。
在一个具体的实施例中，考虑到最小二乘估计在线性规划问题中对粗差很敏感，通过引入L1范数最小法来探测粗差并确定第一次平差的残差值V，公式如下：
其中，V为8x1残差矩阵；B为已知的8x3系数矩阵；X为3x1未知参数向量；L为8x1已知观测矩阵；观测量的数量，取值为8。根据单纯形解法求出初始残差矩阵/>。
在本实施例中，所述的权值P为权矩阵，采用表示权矩阵中的某个元素，/>为权因子，/>表示对角矩阵；
其中，通过Bisquare权函数构造得到：
。
在本实施例中，结合权值P并利用最小二乘准则得到
冰川三维流速矢量的最佳估计值，法方程公式如下：
其中
冰川三维流速矢量的最佳估计值可由以下公式计算得到：
其中，、/>和/>分别代表东西方向、北南方向和垂直方向的流速最佳估计值。
在冰川表面三维流速方程求解中，引入L1范数最小法探测粗差，对粗差赋予更小的权重再采用加权最小二乘法进行解算，以抵制观测值中可能存在粗差的影响，保证解算参数的可靠性。
在一个具体的实施例中，为进一步探究冰川物质平衡，将由垂直位移量求取的垂直速率进一步分解为由沿着冰川表面斜坡运动引起的表面平行流（SPF）和由于冰川表面高程变化引起的非表面平行流分量（nSPF），而计算表面平行流和非表面平行流分量，具体如下：
引入冰川运动学边界条件与冰川表面高程变化率，表示为，将计算公式如下：
等价于：
其中，为垂直冰粒速度；/>为水平速度矢量；/>表示某一高度空间分辨率DEM数据进行高斯平滑滤波得到的表面坡度，/>，/>表示高程，/>表示沿东西方向，/>表示沿北南方向；/>表示质量平衡率；
根据Sorge恒定密度定律，忽略速度的垂直剖面，将冰川厚度变化转换为等效质量，将由垂直向位移量求取的垂直速率分解为表面平行流和非表面平行流分量；
其中，表面平行流的计算公式为，非表面平行流分量的计算公式为。
本发明效果可以通过以下实验结果进一步说明。
实验描述：
根据2018/8/11~2018/9/20时间段Sentinel-1的2景升轨数据和2景降轨数据、2景GF-3号升轨数据和2景Sentinel-2数据计算得到不同方向的流速，Sentinel-1的升轨数据的流速均值分别为-0.09（方位向）和-0.01m/d（视线向）；Sentinel-1的降轨数据的流速均值分别为-0.05（方位向）和-0.03m/d（视线向）；GF-3号升轨数据的流速均值分别为-0.066（方位向）和0.01m/d（视线向）；Sentinel-2数据的流速均值分别为-0.12（东西方向）和0.05m/d（北南方向）。
融合多源SARSAR影像对数据与光学影像对数据的各方向位移观测与地表三方向位移量建立冰川表面三维流速方程，求解得到东西、南北、垂直向流速，流速均值分别为0.04、-0.08和-0.02m/d。
根据本发明的方法将垂直方向流速分解为表面平行流和非表面平行流，流速均值分别为-0.01和-0.01m/d，其中冰川中下游消融速度相较于上游更快。
为了评估加入多源数据对解算的影响，分别求取仅用Sentinel-1、使用Sentinel-1和Sentinel-2、以及使用Sentinel-1、Sentinel-2和GF-3联合解算的三维误差椭球，如图3所示。融合多源数据比仅使用单一数据的解算精度，在东-西、北-南和垂直方向上分别提高了51%、30%和31%。
为评估基于L1范数和最小二乘解算对结果的影响，求取两种方法解算在三维方向上的标准偏差在东西方向分别为0.0254和0.0189，在南北方向为0.0553和0.0154，在垂直方向为0.0239和0.0197。本发明的方法具有较好的抗差性，可以抵制一定数量的观测误差的影响。
实施例3
本实施例还提供了一种计算机系统，包括存储器、处理器以及存储在存储器上并可在处理器上运行的计算机程序，所述的处理器执行所述的计算机程序时，实现如实施例1或实施例2所述的方法的步骤。
其中，存储器和处理器采用总线方式连接，总线可以包括任意数量的互联的总线和桥，总线将一个或多个处理器和存储器的各种电路连接在一起。总线还可以将诸如外围设备、稳压器和功率管理电路等之类的各种其他电路连接在一起，这些都是本领域所公知的，因此，本文不再对其进行进一步描述。总线接口在总线和收发机之间提供接口。收发机可以是一个元件，也可以是多个元件，比如多个接收器和发送器，提供用于在传输介质上与各种其他装置通信的单元。经处理器处理的数据通过天线在无线介质上进行传输，进一步，天线还接收数据并将数据传送给处理器。
还提供了一种计算机可读存储介质，其上存储有计算机程序，所述的计算机程序被处理器执行时，实现如实施例1或实施例2所述的方法的步骤。
即，本领域技术人员可以理解，实现上述实施例方法中的全部或部分步骤是可以通过程序来指令相关的硬件来完成，该程序存储在一个存储介质中，包括若干指令用以使得一个设备(可以是单片机，芯片等)或处理器(processor)执行本申请各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-OnlyMemory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。
显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。
Claims (10)
1.一种基于多源遥感数据的高时间分辨率冰川厚度反演方法，其特征在于：所述的方法包括步骤如下：
融合第三SAR影像对、光学影像对的各方向偏移量、与冰川三维流速矢量建立冰川表面三维流速方程；所述的第三SAR影像与光学影像对的偏移量结果分辨率一致；
求解冰川表面三维流速方程，具体通过引入L1范数最小法来探测粗差并确定第一次平差的残差值V；对得到的残差值V进行赋权得到权值P；结合权值P并利用最小二乘准则得到三维地表形变的最佳估计值；
所述的三维地表形变的最佳估计值包括东西向位移量、南北向位移量、垂直向位移量；
将由垂直向位移量求取的垂直速率分解为由沿着冰川表面斜坡运动引起的表面平行流和由于冰川表面高程变化引起的非表面平行流分量；
结合东西向位移量、南北向位移量计算得到的水平速度矢量，与某一高度空间分辨率DEM数据进行高斯平滑滤波得到的表面坡度，计算表面平行流和非表面平行流分量。
2.根据权利要求1所述的基于多源遥感数据的高时间分辨率冰川厚度反演方法，其特征在于：所述的第三SAR影像对通过以下处理得到：
获取第一平台遥感数据和第二平台遥感数据组成第一SAR影像对，通过设置时间基线阈值、空间基线阈值选择得到最优的第一SAR影像对；
根据某一高度空间分辨率DEM数据对最优的第一SAR影像对进行精确配准，得到第二SAR影像对；
采用强度跟踪算法对第二SAR影像对进行偏移量的估算，每对SAR影像得到视线向、方位向的第一偏移量；对第一偏移量进行地理编码并依据最邻近法重采样至网格单元中，得到与光学影像对的偏移量结果分辨率一致的第三SAR影像对。
3.根据权利要求1所述的基于多源遥感数据的高时间分辨率冰川厚度反演方法，其特征在于：在将光学影像对与第三SAR影像对融合之前，先采用强度跟踪算法对光学影像对进行偏移量的估算，剔除误差，得到光学影像对的东西向、南北向的第二偏移量。
4.根据权利要求1所述的基于多源遥感数据的高时间分辨率冰川厚度反演方法，其特征在于：设方位向为正，视线向为正，垂直向上、正东、正北方向为正；融合第三SAR影像对与光学影像对的各方向偏移量，与冰川三维流速矢量的东西方向、北南方向和垂直方向建立冰川表面三维流速方程如下：
其中：
式中，A和D分别代表第一平台遥感数据的升轨数据和降轨数据，代表第二平台遥感数据的升轨数据；/>代表光学影像；/>是SAR卫星的方位角，即卫星飞行方向与正北方向沿顺时针的夹角；/>是卫星的入射角；/>和/>分别表示由SAR影像对得到的地面点的视线向位移量和方位向的第一偏移量；/>和/>分别表示由光学影像对得到的地面点的东西向位移量和南北向的第二偏移量；/>是在/>时间间隔内计算的偏移量；B为系数矩阵，X为观测点三维流速，L为观测量；/>表示冰川三维流速矢量的东西方向，/>表示冰川三维流速矢量的北南方向，/>表示冰川三维流速矢量的垂直方向。
5.根据权利要求4所述的基于多源遥感数据的高时间分辨率冰川厚度反演方法，其特征在于：通过引入L1范数最小法来探测粗差并确定第一次平差的残差值V，公式如下：
其中，V为8x1残差矩阵；B为已知的8x3系数矩阵；X为3x1未知参数向量；L为8x1已知观测矩阵； 观测量的数量。
6.根据权利要求5所述的基于多源遥感数据的高时间分辨率冰川厚度反演方法，其特征在于：所述的权值P为权矩阵，采用表示权矩阵中的某个元素，/>为权因子，/>表示对角矩阵；
其中，通过Bisquare权函数构造得到：
。
7.根据权利要求5所述的基于多源遥感数据的高时间分辨率冰川厚度反演方法，其特征在于：结合权值P并利用最小二乘准则得到冰川三维流速矢量/>的最佳估计值，法方程公式如下：
其中：
冰川三维流速矢量的最佳估计值可由以下公式计算得到：
其中，、/>和/>分别代表东西方向、北南方向和垂直方向的流速最佳估计值。
8.根据权利要求7所述的基于多源遥感数据的高时间分辨率冰川厚度反演方法，其特征在于：计算表面平行流和非表面平行流分量，具体如下：
引入冰川运动学边界条件与冰川表面高程变化率，将计算公式如下：
等价于：
其中，为垂直冰粒速度；/>为水平速度矢量；/>表示某一高度空间分辨率DEM数据进行高斯平滑滤波得到的表面坡度，/>，/>表示高程，/>表示沿东西方向，/>表示沿北南方向；/>表示质量平衡率；
根据Sorge恒定密度定律，忽略速度的垂直剖面，将冰川厚度变化转换为等效质量，将由垂直向位移量求取的垂直速率分解为表面平行流和非表面平行流分量；
其中，表面平行流的计算公式为，非表面平行流分量的计算公式为。
9.一种计算机设备，包括存储器和处理器，所述存储器存储有可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现如权利要求1至8中任一项所述方法的步骤。
10.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1至8中任一项所述方法的步骤。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310497009.7A CN116485857A (zh) | 2023-05-05 | 2023-05-05 | 一种基于多源遥感数据的高时间分辨率冰川厚度反演方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310497009.7A CN116485857A (zh) | 2023-05-05 | 2023-05-05 | 一种基于多源遥感数据的高时间分辨率冰川厚度反演方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116485857A true CN116485857A (zh) | 2023-07-25 |
Family
ID=87217649
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310497009.7A Pending CN116485857A (zh) | 2023-05-05 | 2023-05-05 | 一种基于多源遥感数据的高时间分辨率冰川厚度反演方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116485857A (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117236668A (zh) * | 2023-11-15 | 2023-12-15 | 山东锋士信息技术有限公司 | 基于供给与消耗的区域水资源配置分析决策方法及系统 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108919262A (zh) * | 2018-04-27 | 2018-11-30 | 中国国土资源航空物探遥感中心 | Dem辅助强度相关的冰川表面运动三维矢量反演方法 |
| CN109472810A (zh) * | 2018-07-10 | 2019-03-15 | 湖南科技大学 | 一种基于遥感图像的冰川流速可视化提取方法 |
| CN114924268A (zh) * | 2022-04-01 | 2022-08-19 | 长安大学 | 一种大量级形变的三维位移反演方法 |
| WO2022213673A1 (zh) * | 2021-04-06 | 2022-10-13 | 中国矿业大学 | 融合无人机dom和星载sar影像的地表三维形变提取方法 |
-
2023
- 2023-05-05 CN CN202310497009.7A patent/CN116485857A/zh active Pending
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108919262A (zh) * | 2018-04-27 | 2018-11-30 | 中国国土资源航空物探遥感中心 | Dem辅助强度相关的冰川表面运动三维矢量反演方法 |
| CN109472810A (zh) * | 2018-07-10 | 2019-03-15 | 湖南科技大学 | 一种基于遥感图像的冰川流速可视化提取方法 |
| WO2022213673A1 (zh) * | 2021-04-06 | 2022-10-13 | 中国矿业大学 | 融合无人机dom和星载sar影像的地表三维形变提取方法 |
| CN114924268A (zh) * | 2022-04-01 | 2022-08-19 | 长安大学 | 一种大量级形变的三维位移反演方法 |
Non-Patent Citations (5)
| Title |
|---|
| JIALIANG LIU ET AL: "Three-Dimensional Flow Velocity Estimation of Mountain Glacier Based on SAR Interferometry and Offset-Tracking Technology: A Case of the Urumqi Glacier No.1", 《WATER》, pages 1 - 25 * |
| XUEMIN XING ET AL: "Slow Deformation Time-Series Monitoring for Urban Areas Based on the AWHPSPO Algorithm and TELM: A Case Study of Changsha, China", 《REMOTE SENSING》, pages 1 - 20 * |
| YANJUN LI ET AL: "Unmanned Aerial Vehicle Remote Sensing for Antarctic Research", 《IEEE》, pages 73 - 93 * |
| YUELING SHI ET AL: "Investigating the Intra-Annual Dynamics of Kunlun Glacier in the West Kunlun Mountains, China, From Ascending and Descending Sentinel-1 SAR", 《IEEE》, pages 1272 - 1282 * |
| 李刚 等: "基于多源 SAR 数据的 2022 年门源 Ms 6.9 级地震同震破裂模型反演研究", 《武汉大学学报(信息科学版)》, pages 1 - 19 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117236668A (zh) * | 2023-11-15 | 2023-12-15 | 山东锋士信息技术有限公司 | 基于供给与消耗的区域水资源配置分析决策方法及系统 |
| CN117236668B (zh) * | 2023-11-15 | 2024-03-08 | 山东锋士信息技术有限公司 | 基于供给与消耗的区域水资源配置分析决策方法及系统 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110058236B (zh) | 一种面向三维地表形变估计的InSAR和GNSS定权方法 | |
| CN109212522B (zh) | 一种基于双基星载sar的高精度dem反演方法及装置 | |
| CN109100723B (zh) | 基于多普勒天气雷达数据的高空风反演方法 | |
| CN104597471A (zh) | 面向时钟同步多天线gnss接收机的定向测姿方法 | |
| CN109031301A (zh) | 基于PSInSAR技术的山区地形形变提取方法 | |
| CN105301601A (zh) | 一种适用于全球区域的gnss电离层延迟三维建模方法 | |
| CN114488228B (zh) | 一种适用于动态载体平台的gnss多路径误差削弱方法 | |
| CN111856459B (zh) | 一种改进的DEM最大似然约束多基线InSAR相位解缠方法 | |
| CN116485857A (zh) | 一种基于多源遥感数据的高时间分辨率冰川厚度反演方法 | |
| CN110780327A (zh) | 基于星载ais与红外相机的海上目标协同定位方法 | |
| CN115201825B (zh) | 一种InSAR震间形变监测中的大气延迟校正方法 | |
| CN113238228B (zh) | 基于水准约束的三维地表形变获取方法、系统及装置 | |
| CN110146904B (zh) | 一种适用于区域电离层tec的精确建模方法 | |
| CN106157258B (zh) | 一种星载sar图像几何校正方法 | |
| CN111722250B (zh) | 基于gnss时间序列的地壳形变影像共模误差提取方法 | |
| CN107909606A (zh) | 一种sar图像配准联系点粗差剔除方法 | |
| CN109783846B (zh) | 基于gnss海洋浮标的潮位测量不确定度评定方法 | |
| CN114488227B (zh) | 一种基于空间相关性的多路径误差改正方法 | |
| CN116052009A (zh) | 基于InSAR形变的GNSS测站布局方法及其应用 | |
| CN115097450A (zh) | 跨轨道高分三号sar偏移量大梯度滑坡形变估计方法 | |
| Luo et al. | Monitoring Surface Deformation of Transmission Corridors in Mountain Areas Based on SBAS-INSAR | |
| CN115201823B (zh) | 一种利用BDS-InSAR数据融合的地表形变监测方法 | |
| CN111443366A (zh) | 一种gnss区域网中异常点探测方法及系统 | |
| CN111308570A (zh) | 基于载波相位微分速度构建全球重力场的方法 | |
| Yin et al. | Constrained robust unscented Kalman filter for BDS navigation in dense urban areas |
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 |