CN107421496A - 一种高精度的水面高程提取算法 - Google Patents

一种高精度的水面高程提取算法 Download PDF

Info

Publication number
CN107421496A
CN107421496A CN201710617728.2A CN201710617728A CN107421496A CN 107421496 A CN107421496 A CN 107421496A CN 201710617728 A CN201710617728 A CN 201710617728A CN 107421496 A CN107421496 A CN 107421496A
Authority
CN
China
Prior art keywords
waveform
mrow
correction
level
water level
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
CN201710617728.2A
Other languages
English (en)
Other versions
CN107421496B (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.)
Qingdao Meiji Marine Geographic Information Technology Co Ltd
Original Assignee
Qingdao Meiji Marine Geographic Information Technology 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 Qingdao Meiji Marine Geographic Information Technology Co Ltd filed Critical Qingdao Meiji Marine Geographic Information Technology Co Ltd
Priority to CN201710617728.2A priority Critical patent/CN107421496B/zh
Publication of CN107421496A publication Critical patent/CN107421496A/zh
Application granted granted Critical
Publication of CN107421496B publication Critical patent/CN107421496B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/14—Receivers specially adapted for specific applications
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C5/00—Measuring height; Measuring distances transverse to line of sight; Levelling between separated points; Surveyors' levels
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01F—MEASURING VOLUME, VOLUME FLOW, MASS FLOW OR LIQUID LEVEL; METERING BY VOLUME
    • G01F23/00—Indicating or measuring liquid level or level of fluent solid material, e.g. indicating in terms of volume or indicating by means of an alarm
    • G01F23/22—Indicating or measuring liquid level or level of fluent solid material, e.g. indicating in terms of volume or indicating by means of an alarm by measuring physical variables, other than linear dimensions, pressure or weight, dependent on the level to be measured, e.g. by difference of heat transfer of steam or water
    • G01F23/28—Indicating or measuring liquid level or level of fluent solid material, e.g. indicating in terms of volume or indicating by means of an alarm by measuring physical variables, other than linear dimensions, pressure or weight, dependent on the level to be measured, e.g. by difference of heat transfer of steam or water by measuring the variations of parameters of electromagnetic or acoustic waves applied directly to the liquid or fluent solid material
    • G01F23/284—Electromagnetic waves

Abstract

本发明公开了一种高精度的水面高程提取算法，该算法包括通用类波形的模拟；通用类波形与Level 1b干涉波形互相关重跟踪算法；大气及地球物理修正；高程计算与异常值排除；基于SSA的水位滤波与噪声分析。由于卫星高度计本身的限制和复杂的地理环境，Cryosat‑2 Level 2高度计数据产品在内陆或近岸区域存在大量的数据缺失和精度损失，本发明通过改进高度计回波波形重跟踪算法，提取高精度的水位序列产品，弥补当前内陆水体或近岸测高数据数量少、精度不高的缺陷。

Description

一种高精度的水面高程提取算法
技术领域
本发明属于卫星遥感与地理信息技术领域，具体涉及一种基于卫星测高的水面高程提取方法。
背景技术
卫星测高被用作提取水位(水面高程)的一种有效工具，其测距精度从分米级到厘米级。水位通常通过卫星高度计理论回波模型与接收波形拟合来计算，这个过程称之为波形重跟踪(Waveform retracking)。卫星高度计的二级GDR数据在发布前已经经过波形重跟踪处理，对受到陆地地形影响较大的波形进行了可能的数据剔除，这导致有些原本过境湖面的测高卫星在此出现数据缺失，使卫星高度计数据观测的湖泊数量十分有限。由于通常雷达高度计(T/P、Jason、HY-2)照射的地面足印直径约为15-20km，在这么大的范围内，周围的地形将有可能污染观测到的近岸波形。如果波形被严重污染，波形重跟踪器(Waveformretracker)可能会给出不精确的或无效的水位值。
由于卫星雷达高度计自身的局限和复杂的地理环境，测高卫星在海洋近岸的回波波形也受到陆地地形的影响，导致观测数据质量较差，为了增加和提高ERS、Topex/Poseidon、Jason和Envisat等传统脉冲式高度计数据的数量和质量，研究人员针对受陆地影响的海洋近岸回波提出了多种波形重跟踪方法，这些重跟踪算法也被应用到湖泊水位的提取上。应用较广泛的有Beta-5参数法(Martin et al.,1983)、重心偏差OCOG法(Winghamet al.,1986)、海洋算法(Callahan et al.,2004)和阈值法(Davis et al.,1997)等，基于这些方法又提出了许多改进的算法，比如许厚泽等(2007)结合Beta-5算法和OCOG算法对中国南海区域四轨Topex/Poseidon和ERS-1的波形进行处理，并将结果与验潮站的实测数据进行了比较，用于重构近海平均海平面；Yang et al.,(2008)利用Jason-1高度计在中国近海一年的测量波形，采用升降轨交叉点的海面高度异常值及验潮站的实测水位比较了海洋算法、OCOG算法、Ice-2算法和阈值算法四种波形重构算法的结果，指出在近海的波形处理方面OCOG算法优于其他算法；Guo et al.,2006和Hwang et al.,2006提出了一种新的算法，通过从有多个前沿的回波中提取出子波形，再判断并确定与大地水准面模型结果最合适的海面高度结果。该方法的本质还是阈值算法，是不基于模型的算法，缺乏明确的物理概念，而且只能由波形得到卫星到海面的距离这一个量，没有有效波高和后向散射系数；Gomez-Enri et al.,2009针对近海Brown和尖峰混合的回波，提出了实验性的混合波形重跟踪算法；Yang et al.(2012)基于波形分类和子波形的波形重跟踪算法(OceanCS)，该算法结合了海洋算法和重心偏差算法的优势。
虽然以上波形重跟踪技术已被证明能够有效改善测高数据的精度，但是没有一种算法能够使环境极为复杂的内陆地区各种测高波形的重跟踪效果均达到最优。这是因为各种波形重跟踪方法都是针对特定地区的应用提出的，只对某些特定的波形有较好的重跟踪效果，例如Beta-5参数法适用于海洋波形，阈值法是为了处理冰面波形提出的。地理环境非常复杂的内陆湖泊存在多种多样的波形，仅用一种算法难以满足各种波形的重跟踪，也就不能最大限度地提高测高数据的数量和质量。若能根据波形形状自动选取最佳重跟踪算法，则有可能极大提高可用的测高数据。Freeman et al.,2006提出了一个处理近海回波的专家系统，通过把波形分为十类，再分别采用四种不同的波形重构方法(海洋算法、OCOG算法、Ice-2算法和阈值算法)进行计算。但这仍有不足，由于是专家系统，该算法也需要一定的先验知识及人为干预。而且，Berry et al.(2005)也曾依据观测波形的类别，使用几种不同的重跟踪算法去改进观测精度，但与位于玛瑙的验潮站实测数据比较发现，其与Frappart et al.(2006)仅使用Ice重跟踪算法并无明显优势。这些算法都是针对ERS、Envisat、Jason-1/2等传统脉冲有限式高度计的回波波形特征提出的，对于Cryosat-2干涉模式的波形特征并不完全适用。
在2010年，搭载着合成孔径干涉雷达高度计SIRAL的Cryosat-2卫星发射，它是一种新型的Ku波段雷达高度计，也是目前国际上第一颗应用干涉技术的测高卫星，标志着雷达高度计技术的一大突破。这种高度计能够采用三种不同的模式运行：低分辨率(LR)模式、合成孔径雷达(SAR)模式和SAR interferometry(SARin)模式。对于具有复杂地形的地区(比如青藏高原)，SIRAL使用SARin模式运行，沿轨分辨率大约是300m，但是由于另外一个干涉天线的使用，它能进行垂直轨迹(Cross-track)的坡向改正(Wingham et al.,2006)。在SARin模式中，Cryosat-2数据不仅具有更高的沿轨分辨率，而且相较ERS、Envisat、Jason等传统高度计也具有不同的波形形状，它们的波形后缘具有更快的下降速率，并且具有多个可分辨的波峰，可以预期Cryosat-2能够找到更加可靠的水位采样。
目前存在的众多重跟踪算法(比如：ocean、Ice-1、Ice-2、sea-ice和OceanCS等)几乎都是针对不受污染的传统脉冲有限式高度计的波形特征提出的，而对Cryosat-2卫星受污染的干涉回波波形尚无法很好地解决。
发明内容
本发明的目的是提供一种高精度的海面高程提取与验证算法，能够适用Cryosat-2干涉模式的波形特征，以弥补现有技术的不足。
在Cryosat-2的标准Level 2处理器中，对获取的Level 1b波形进行了标准化的波形重跟踪处理，每个回波波形仅仅能够获取一个高程。如果波形中存在两个以上的、可分辨的不同 波峰，标准的重跟踪处理器可能重跟踪到波形中那个错误的波峰，导致测量值无效，基于获取的Level 1b观测波形与模拟的Cryosat-2 SARin通用类波形的互相关关系的波形重跟踪处理算法，能够计算每个波形存在的多个波峰对应的高程，例如湖泊水位一个高程，湖泊周围的地形一个高程。由于同一观测周期内湖泊水位沿轨保持常数，而周围地形是随机变化的，因此通过多个沿轨测量值可以有效移除地形的高程而找到最优的湖泊水位估计值。
基于上述基本思路，本发明将针对Cryosat-2雷达高度计Level 1b干涉模式的波形数据开展波形模拟、重跟踪与湖泊水位的精提取研究。最后，利用奇异谱分析技术和外部其他实测数据对水位进行噪声评估。
因此，基于上述目的和思路，本发明采取的具体技术方案为：
一种高精度的水面高程提取算法，包括以下步骤：
(1)通用类波形的模拟；
(2)通用类波形与Level 1b干涉波形互相关重跟踪算法；
(3)大气及地球物理效应修正；
(4)高程计算与异常值排除；
(5)基于SSA的水位滤波与噪声分析。
进一步的，上述步骤(1)包括不同散射体去斜回波的模拟、去斜干涉回波波形的模拟和最优波形后向散射系数的设定。
进一步的，上述步骤(2)具体为，基于步骤(1)模拟的通用类波形与观测的Level1b干涉波形进行互相关关系，找到观测波形中每个波峰位置与参考波形的偏离程度，最终重跟踪获取距离改正值rrc；基于L1b波形与模拟的类波形的互相关关系，对L1b的距离改正值rrc进行求解；
观测的Level 1b干涉波形与模拟的通用类波形的互相关关系如公式1所示：
其中，ρ(n)是L1b波形与模拟波形的互相关函数，wL1b和wsim分别代表L1b波形和模拟波形的能量，通过互相关函数可以得到多个波形中波峰的位置偏移np，然后将np转化成距离改正值，如公式(2)
进一步的，上述步骤(4)中，基于公式(3)，对重跟踪及大气、物理校正后湖面高程进 行计算：
h＝halt-(rwd+rrc+rgc)-hgeoid (3)
其中，h代表计算的湖面水位，halt代表卫星的高度，hgeoid代表大地水准面高程，rwd代表雷达信号的参考距离reference range，直接通过观测的L1b波形的窗口延迟时间计算，rrc代表上述重跟踪获取的距离改正值，rgc代表大气及地球物理改正。
进一步的，上述大气及地球物理效应的改正计算如公式(4)：
rgc＝Δrdry+Δrwet+Δriono+ΔrSSB+Δrset+Δrpol (4)
其中，rgc代表总的大气及地球物理效应改正；Δrdry是干对流层改正；Δrwet是湿对流层改正；Δriono是电离层改正；ΔrSSB是海况(电磁)偏差改正；Δrset是固体潮改正；Δrpol是极潮改正。
进一步的，上述提取算法得到相应的高程后，再进行验证，即验证了本算法能够得到高精度的水面高程。
本发明的优点和有益效果：
(1)由于卫星高度计本身的限制和复杂的地理环境，Cryosat-2 Level 2高度计数据产品在内陆或近岸区域存在大量的数据缺失和精度损失，本发明通过改进高度计回波波形重跟踪算法，提取高精度的水位序列产品，弥补当前内陆水体(比如湖泊)或近岸测高数据数量少、精度不高的缺陷。
(2)通用类波形的模拟是波形重跟踪成功实现的关键，它涉及星上、星下的一系列信号调制过程，我们对Cryosat-2干涉模式下的数据处理技术进行详细总结与改进，实现对观测波形的理论模拟，为实际观测回波中的多个波峰的探测提供了可能。
(3)GIM模型提供了从GPS地面跟踪站到GPS卫星(20000km)天顶方向的TEC，而Cryosat-2卫星轨道高度仅为717.2km，利用IRI2016模型计算Cryosat-2卫星轨道以上与高度20000km以下区间内的TEC，由IRI2016模型在Cryosat-2卫星轨道以下和以上的TEC之比，计算GIM模型在Cryosat-2卫星轨道以下的TEC，以此实现现有电离层GIM模型的改进。经过与Jason-2 GDR校正值比较发现，改进的GIM校正值与Jason-2双频校正结果差值的绝对值在2cm之内。改进的GIM方法可以作为一种普适的、精度较好的电离层延迟误差校正方法。
(4)由于Cryosat-2回波波形多个波峰的存在，导致重跟踪的湖面高程中存在的异常值比例很高(有时超过70％)，传统的RANSAC算法和3倍中误差算法无法满足要求，本发明根据沿轨轨迹上的水位保持常数这一标准，总结了一种新的异常值移除算法，实现湖泊水位的精确提取。另外，湖泊水位除了长期线性趋势外，还表现为明显的周期波动性，通过对水位序列的SSA分析和验潮站比较，实现对了Cryosat-2高度计重跟踪数据噪声评估和精度验证。
附图说明
图1是本发明基于卫星测高的水面高程提取流程图。
图2是GPS和Cryosat-2高度计轨道高度对比示意图。
图3是本发明采用的对应多个波峰的水面高程精提取的方法实例结果图。
图4是基于奇异谱分析(SSA)的青海湖测高水位滤波图。
图5是使用实例数据对于高度计提取水位的精度验证评估结果图。
图6是水文站获取的湖泊水位与SSA提取的湖泊水位的比较结果图。
具体实施方式
以下通过具体实施例并结合附图对本发明进一步解释和说明。
实施例：
本实施例的具体技术方案流程如图1所示，主要包括以下五个步骤：
(1)通用类波形的模拟；
(2)通用类波形与Level 1b干涉波形互相关重跟踪算法；
(3)大气及地球物理效应修正；
(4)高程计算与异常值排除；
(5)基于SSA的水位滤波与噪声分析。
以下是对上述流程的具体分析：
(1)通用类波形的模拟
通过对Cryosat-2星上、星下回波信号干涉处理过程的研究，我们在模拟通用类波形的过程中主要解决以下三个方面的信号调制问题：不同散射体去斜回波(derampedecho)的模拟、去斜回波干涉波形(waveform from deramped echoes)的模拟、最优波形后向散射系数的设定。最终经过一系列信号模拟与调制，获取能够反映真实湖面波形的通用类模拟波形，并以此重跟踪实际观测的L1b波形。
通过对Cryosat-2SARin模式星上、星下的信号处理过程初步的理论推导，我们试验通用类波形的模拟，主要包括：
a)不同散射体去斜回波的模拟
主要参考星上去斜信号的处理，它将接收的信号由瞬时可变的频率调制为固定频率，这里主要考虑多普勒效应、天线增益、及后向散射系数。最终去斜回波(deramped echo)可按照公式5进行计算：
其中，固定去斜频率fd＝Q(tr-tw)，多普勒频移fd,Doppler＝(D-1)fc，相位fc为SIRAL的中心频率13.575Ghz，为多普勒效应，V为水平卫星速度(约7500m/s)，nv和ns分别为归一化的卫星速度矢量和散 射体方向矢量，c为光速，频率坡度τt为传播的脉冲长度49μs，tr为回波到达时间，B为信号带宽，振幅Am与后向散射系数σ0和天线增益G有关，θ和ν分别为相对天线孔径的极角和方位角。
b)去斜干涉回波波形的模拟——多视处理
主要参考星下SAR干涉的多视处理(Multilooking processing)。过程包括：
第一步，将上面模拟的去斜回波echos控制到星下点的地面位置Ground location(称为多视配准)。
第二步，将去斜回波进行方位向快速傅里叶变换形成波束。
第三步，对去斜波束进行多普勒改正移除多普勒频移。
第四步，对去斜波束进行倾斜距离改正，移除参考距离与观测距离的差异。
第五步，为避免后续处理中方位向分辨率的损失，对去斜波束信号通过零补齐进行2倍分辨率的重采样。
第六步，对去斜波束信号进行距离向快速傅里叶变换形成单视波形。
第七步，通过平均波数簇中照射到相同条带的多个单视波形的功率形成平均波形，即获取得单个天线获取的去斜多视回波波形。
第八步，由于Cryosat-2具有两个干涉天线，通过两个天线获取的去斜多视回波进行干涉处理(共轭相乘)，获取干涉回波波形，它将被用于与实测L1b波形的互相关重跟踪。
c)最优波形后向散射系数的设定
实际上，经过上面星上、星下的信号调制，理论上应该已经获得了模拟的通用干涉类波形，但是为了获得最优的模拟波形，应该选择最优的后向散射系数，如公式6：
其中，αi为相对于湖面的入射角，σ为水面坡度的RMS。
通过一系列的σ设置(0.005,0.02)，当使所有轨迹的重跟踪后高程偏差最小时的σ为设定的最优后向散射系数。
(2)波形重跟踪
基于模拟的通用类波形与观测的L1b波形进行互相关关系，找到观测波形中每个波峰位置与参考波形的偏离程度，最终重跟踪获取距离改正值rrc。
为了得到比Cryosat-2卫星原始L1b数据的range分辨率0.47m更加精确的湖泊水位估计，首先考虑对观测的L1b数据和模拟的通用类波形以原采样10倍的分辨率进行重采样。理论上，这样可使最终获取的湖泊水位精度达到0.047m。下面基于L1b波形与模拟的类波形的互相关关系，对L1b的range改正值rrc进行求解。
观测的L1b波形与模拟的通用类波形的互相关关系如公式7所示：
其中，ρ(n)是L1b波形与模拟波形的互相关函数，wL1b和wsim分别代表L1b波形和模拟波形的能量。通过互相关函数可以搜索到波形中波峰的位置偏移np，他们要么代表湖泊的水位，要么代表陆地的高程，可能多个波峰同时找到。然后将np转化成range的改正值，如公式(8)
(3)大气及地球物理效应的改正
假设星地之间不存在大气(即真空环境)，那么水面高度即为卫星高度(alt)和卫星测距(Retracked range)之差。但有研究指出，Retracked range受到大气环境效应(对流层、电离层)及湖面地球物理效应(海况偏差、固体潮、极潮)的影响，通过干涉技术，在短基线情况下脉冲往返并无法完全抵消大气及地球物理效应的影响(Joana et al.,2015)。由于干涉天线的入射角很小，这些误差的性质与它们在传统脉冲有限式星下点测高是是一样的，总的大气及地球物理效应改正如公式(9)：
rgc＝Δrdry+Δrwet+Δriono+ΔrSSB+Δrset+Δrpol (9)
其中，rgc代表总的大气及地球物理效应改正；Δrdry是干对流层改正；Δrwet是湿对流层改正；Δriono是电离层改正；ΔrSSB是海况(电磁)偏差改正；Δrset是固体潮改正；Δrpol是极潮改正。
公式(9)中的各变量可以从Cryosat-2L1b数据文件中直接读取，并且具有较高的精度，比如其中干对流层改正是利用湖泊表面的气压值来计算，目前用模式再分析数据计算的Δrdry具有毫米级的精度，但由于Cryosat-2L1b数据中对大部分湖泊都没有提供相应的湿对流层、电离层以及海况偏差的改正，而它们却是影响测距精度的最主要误差源，比如:湿对流层改正与传播路径上的水汽压力有关，它对测高的影响在-0.9～35.7cm之间，电离层延迟随昼夜季节和纬度有显著的变化，对测高距离的测量影响为0.2～5cm，而由海面(湖面)存在波浪而引起的计算海面高度的不确定性为2cm，因此有必要对这三项测距误差进行单独计算。另外，Jason-2卫星由NASA与CNES共同设计发射，目前仍在轨运行，我从Jason-2二级产品数据读取相应湿对流层、电离层、海况偏差校正的辐射计反演值和模式计算值，为了保证计算精度，将我们反演的大气与地球物理改正值验证后对Retracked range进行改正。
a)湿对流层改正
大气湿对流层延迟校正主要有两种方法：第一种是采用校正微波辐射计测得的亮温，根据亮温与路径延迟的转换关系进行校正；第二种方法是基于当地气压，利用校正模型计算得到。但不管用哪种方法，湿对流层校正算法是相同的。由于Cryosat-2高度计数据中在没有微波辐射计观测值，而美国国家环境预报发布的NCEP数据，提供了高精度的相对湿度、液态水含量、温度、水蒸气分压等大气剖面数据，它为本发明通过第二种方法计算湿对流层改正值提供了可能。
湿对流层路径延迟包括两部分：水蒸气导致的路径延迟和云中液态水导致的路径延迟。水蒸气导致的路径延迟为：
ρv＝1.739×109×RH×θ5×exp(-22.64θ) (10)
其中，T为大气温度，RH为相对湿度，θ＝300/T；
云中液态水导致的路径延迟为：
其中，ρa(h)和ρa(h0)分别表示h高度与相对湿度为94％处的绝对湿度，以上参数都可以直接从NCEP中读取，经过必要的地理配准后计算相关改正值。
总的湿对流层改正为：
Δrwet＝PDV+PDL (12)
b)电离层改正-改进的IRI-GIM方法
当测高卫星信号穿过电离层时，会产生各种物理效应，其中最主要的是折射效应，其结果对传播信号产生延迟。雷达信号传播速度的变化程度与电子密度和信号频率相关，而折射弯曲效应对测距结果影响不大，一般可不予考虑。电子密度随太阳活动强度时间地理位置和信号传播路径的改变而变化，在使用中比较复杂。总电子含量(total electroncontent,TEC)是电子密度沿信号传播路径的积分，电离层距离延迟改正(单位：mm)与TEC含量的关系为：
其中，f为为测高卫星的信号频率(单位：Hz)，TEC的变化范围为0～100TECU(1TECU＝1016/m2)。
迄今尚无法建立TEC的严格计算公式，实用中主要根据全球各电离层观测站积累的大量观测资料建立拟合改正模型。这里我们按照Cryosat-2产品用户手册建议的全球电离层分布模型(GIM，Global Ionosphere Map)对电离层改正量进行计算。我们选取由国际上4个机构NASA/JPL、CODE、ESA和UPC分别建立的GIM模型的加权平均值，其时间分辨率为2h，空间分辨率为纬向2.5°，经向5°，采用450km高度的单层模型数据。GIM模型提供了从地面接收站到GPS卫星天顶方向大约20000km高度范围内的TEC，而高度计卫星大都位于地球表面1000km左右，Cryosat-2卫星仅有717.2km。直接使用GIM模型所提供的电离层TEC值校正雷达高度计会高估电离层路径延迟。图2是GPS和高度计轨道高度对比示意图。
利用IRI-2016模型分别计算Cryosat-2高度计卫星轨道以下(h1，717.2km)和GPS卫星轨道以下(h2，20000km)的TEC，并计算二者之比，以此从GIM值中获取实际用于校正高度计电离层路径延迟的TEC分量TECalt，计算公式如下：
将按比例获得的有效电子含量TECalt带入公式(13)，取Ku波段频率13.58GHz，即可获得电离层路径延迟值。
经过与Jason-2GDR校正值比较发现，改进的GIM校正值与双频校正结果差值的绝对值在2cm之内。改进的GIM方法可以作为一种普适的、精度较好的电离层延迟误差校正方法。
c)海况偏差改正
海况偏差(SSB)是由于海表面存在波浪所导致实际的高度计测距大于真实距离的一种误差，湖泊同样受到这种因素的影响。海况偏差修正是高度计测距误差中最主要的误差源，目前有参数模型和非参数模型，当中等风速和浪高时，参数和非参数模型的修正效果均较好，但是在内陆湖泊地区会出现高风速低海况状况，带来较大误差。由于内陆湖泊一般面积较小，利用参数模型通过多项式的拟合一般会在局部有较好的反演效果，且实现相对容易，因此我们采用戈达德空间飞行中心GSFC推荐的四参数BM4模型计算海况偏差改正值，如公式(15):
ΔrSSB＝Hswh(a+bu+cu2+dHswh) (15)
其中，a、b、c、d四个参数的推荐值分别为-0.046200、-0.002350、0.000103和0.002540；u为通过后巷散射系数反演的风速；Hswh为有效波高。
BM4模型已成为广泛采用的海况偏差经验校正模型，此模型成功地应用在T/P，Jason等高度计GDR数据中，并且之后的众多的非参数模型也是根据与此模型的对比来分析模型的性能的。风场和波场数据是影响有效波高的两个重要因素，实验证明ECMWF的与Jason-2等高度计数据的一致性非常好，因此在构建模型时，我们以ECMWF再分析的有效波高数据作为模型中Hswh的输入变量。
由于湖面回波的特性与海面回波的特性是有差别的，而且雷达传感器在卫星经过水陆交接面时可能工作异常，导致有效波高异常。理想的海况偏差系数应该是非参数的，而非预先假定一个经验模型再利用测高数据去拟合。实际处理中，当观测有效波高大于5m时，本发明在现有非参数模型基础上，构建适用于Cryosat-2高度计的双要素(有效波高、海面风速)无参数模型。目前精度最高的海况校正模型是核平滑非参数(Non-Parametric，NPSSB)模型，最终的修正数值来自于根据风速u和有效波高Hswh建立的SSB查找表来进行双线性插值计算。
与参数模型不同，非参数校正模型没有具体的函数形式。双要素NPSSB模型是基于数据的统计结果构建的，在进行统计分析时，海面的高度差可以通过SSB的差异来表示:
HSSH2-HSSH1＝ΔrSSB2(u,Hswh)-ΔrSSB1(u,Hswh)+ε (16)
上式中，HSSH表示基于大地水准面的湖面高度，其值除了SSB以外其他的误差源均已校正；下标2和1对应不同的观测时间t2和t1，ΔrSSB(u,Hswh)表示包含了风速u和有效波高Hswh的海况偏差；ε表示误差项，其均值为0，包含了时间t2和t1测高数据与地理位置相关的误差、测距过程中其他校正项的误差等。湖面高度的差异可以用交叉点或共线轨道数据来获取。
(4)湖泊水位精提取——高程计算与异常值移除
基于公式(17)，对重跟踪及大气、物理校正后湖面高程进行计算。
h＝halt-(rwd+rrc+rgc)-hgeoid (17)
其中，h代表计算的湖面水位，halt代表卫星的高度，hgeoid代表大地水准面高程，rwd代 表雷达信号的参考距离reference range，它可以直接通过观测的L1b波形的窗口延迟时间计算，rrc代表上面重跟踪获取的range改正，rgc代表大气及地球物理改正。
一个L1b湖泊波形内可能有多个波峰，经过互相关重跟踪后，可能获取多个重跟踪改正值(公式2或8)，通过公式17，最终获得对应的多个湖面高程，如图3所示。但是，仅有一个代表真正的湖面水位，而这可以根据沿轨轨迹上的湖面水位应保持常数这一标准来进行选择。
卫星某个时刻过境湖面时，包含了多个高程的波形沿轨测量值会形成一个湖泊高程集合(图3)，但有时候异常值占总数的比例过大(有时超过70％)，标准的异常值移除算法(eg.RANSAC异常值移除算法)并不适用。为了移除异常值，我们尝试使用下面的异常值移除算法：
第一步，如果大部分沿轨高程值在某一高程间隔(elevation bin，可设置为1m)范围内(如图3横框内)，超出部分被剔除。
第二步，对剩余保留在红框内的高程值进行趋势拟合，有些沿轨L1b波形对应的多个高程值仍然保留在红框内，这时距离高程值拟合趋势最近的高程值代表湖泊水位，其余高程值进行剔除。
第三步，对于剩余高程值再次进行趋势拟合，剔除3倍中误差以外的高程值。最后保留的高程值数量应超过10个以上，这些剩余的高程值的均值作为最后沿轨的高程值。
图3中，竖框内的点代表具有多个波峰的波形对应的不同高程值，同一条轨迹上的湖面水位应当是一常数，因此对于明显偏离平均的点被移除，而横框内的点是异常值移除后保留的高程值(1m的设置范围框内)。
(5)基于SSA的水位滤波与噪声分析
奇异谱分析(Singular Spectrum Analysis,SSA)是一种特别适合于研究周期振荡行为的分析方法。以内陆湖泊为例，湖泊水位变化具有明显的周年/半年的季节振荡，因此借助SSA来对高度计获取的湖泊水位波动信息进行滤波，可有效降低水位序列的噪声水平。
通过精提取的Cryosat-2高度计水位序列、湖泊水文站实测序列与SSA滤波的水位序列的相互比较可以对精提取的高度计水位序列的噪声进行估计及分析。但基于SSA的水位时序滤波过程在实际应用中，它还存在着一系列问题，主要有三个方面的难点：a)窗宽M和前P项主成分的选择；b)缺值点的插补；c)粗差点的剔除。这三个难点的处理方法如下：
a).窗宽M和前P项主成分的选择。在进行奇异谱分析时，窗宽M的选择很关键，窗宽越大谱分辨率越灵敏，能分辨互相靠近的振荡周期，但同时将造成重建序列中震荡周期之间间歇性时段定位粗糙；前几项主成分对应着原序列信号的主要成分，后几项主成分主要是观测中的误差成分，因此，在实际重建成分过程中一般选择前面P项来进行重建，以达到消除噪声的目的。对于M和P值，本发明将采用交叉验证的方法完成上述两个参数的确定，我们以特定的M与P值进行循环试验，在利用数据建立最优模型的过程中，将数据m等分，选取第一等分数据作为验证值，利用剩余数据建模，然后利用模型值与验证数据的差值计算nrms值，以此类推验证计算m次，最后求解m次试验的平均nrms值以此作为该次建模中参数选择优劣程度评判标准。当平均nrms最小时，其对应的M和P就是所选定的最优参数值。该方法的优点在于在进行参数好坏评价的过程中，所有的数据都参与了模型建立，同时所有的数据又都参与了模型的验证。
b).缺值点的插补。以“天”为时间单位形成的长水位时序数据中，不可避免的会出现连续或不连续的缺值点，因此有必要研究对原始数据的插值方法。利用传统插值方法进行序列插值时存在一定的复杂性和不可靠性，一方面，不同序列的数据缺失位置及缺失部分时间跨度不同，另一方面，不同序列的波动特性各异，因此，要得到最优插值结果就要基于数据缺失特性差异及不同序列在不同方向上的内在数据结构差异基础上，选择不同模型及不同模型参数进行复杂的组合插值，这必然导致插值过程变的过于的繁琐，同时对于模型及其参数人 为选择的依赖性也必然使得插值结果会引入不可估计的人为噪声。基于奇异谱插值方法SSA-M取得了很好的效果，本发明将使用该方法完成水位序列缺值点的插补。
c).粗差点的剔除。以“天”为单位形成的水位长时间序列中也不可避免地存在粗差，准确定位并剔除时间序列中的粗差是序列分析的基础。实验表明少量粗差的存在不会影响SSA对整体数据结构的分析，可认为奇异谱分析具有一定的抗粗差能力。奇异谱分析的这些特性保证了可以在不了解序列特性的前提下正确分离有效信息和噪声的能力，在SSA分析时可以对时序中的水位点进行IQR判断，确定该点是否为粗差。
上述算法的结果验证策略
将SSA提取的波动信息(图4)与青海湖水文站历史水位数据进行比较，它们的高程差均值和均方根(root-mean-square-difference,RMSD)分别为0.2202m和0.252m(如图5)，表明高度计数据监测内陆湖泊的水位变化精度可以达到分米量级，并且通过SSA算法提取的水位与水文站的水位显示了很好一致性，其相关系数(R2)在99％置信水平上为0.708(如图6)，这表征了高度计提取的水面高精度和可靠性。

Claims (6)

1.一种高精度的水面高程提取算法，其特征在于，包括以下步骤：
(1)通用类波形的模拟；
(2)通用类波形与Level 1b干涉波形互相关重跟踪算法；
(3)大气及地球物理效应修正；
(4)高程计算与异常值排除；
(5)基于SSA的水位滤波与噪声分析。
2.如权利要求1所述的提取算法，其特征在于，所述步骤(1)包括不同散射体去斜回波的模拟、去斜干涉回波波形的模拟和最优波形后向散射系数的设定。
3.如权利要求1所述的提取算法，其特征在于，所述步骤(2)具体为，基于步骤(1)模拟的通用类波形与观测的Level 1b干涉波形进行互相关关系，找到观测波形中每个波峰位置与参考波形的偏离程度，最终重跟踪获取距离改正值rrc；基于L1b波形与模拟的类波形的互相关关系，对Level 1b的距离改正值rrc进行求解；
观测的Level 1b干涉波形与模拟的通用类波形的互相关关系如公式1所示：
<mrow> <mi>&amp;rho;</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mo>-</mo> <mi>&amp;infin;</mi> </mrow> <mi>&amp;infin;</mi> </munderover> <msub> <mi>w</mi> <mrow> <mi>L</mi> <mn>1</mn> <mi>b</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>m</mi> <mo>)</mo> </mrow> <msub> <mi>w</mi> <mrow> <mi>s</mi> <mi>i</mi> <mi>m</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>n</mi> <mo>+</mo> <mi>m</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>
其中，ρ(n)是L1b波形与模拟波形的互相关函数，wL1b和wsim分别代表Level 1b波形和模拟波形的能量，通过互相关函数可以得到到波形中波峰的位置偏移，他们代表湖泊的水位或者是代表陆地的高程，然后将np转化成距离改正值，如公式(2)
<mrow> <msub> <mi>r</mi> <mrow> <mi>r</mi> <mi>c</mi> </mrow> </msub> <mo>=</mo> <mfrac> <mrow> <msub> <mi>n</mi> <mi>p</mi> </msub> <mo>-</mo> <mn>2550</mn> </mrow> <mn>10</mn> </mfrac> <mo>&amp;CenterDot;</mo> <mn>0.47</mn> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> <mo>.</mo> </mrow>
4.如权利要求1所述的提取算法，其特征在于，所述步骤(4)中，基于公式(3)，对重跟踪及大气、物理校正后湖面高程进行计算：
h＝halt-(rwd+rrc+rgc)-hgeoid (3)
其中，h代表计算的湖面水位，halt代表卫星的高度，hgeoid代表大地水准面高程，rwd代表雷达信号的参考距离reference range，直接通过观测的Level 1b波形的窗口延迟时间计算，rrc代表上述重跟踪获取的距离改正值，rgc代表大气及地球物理改正。
5.如权利要求4所述的提取算法，其特征在于，所述大气及地球物理效应的改正计算如公式(4)：
rgc＝Δrdry+Δrwet+Δriono+ΔrSSB+Δrset+Δrpol (4)
其中，rgc代表总的大气及地球物理效应改正；Δrdry是干对流层改正；Δrwet是湿对流层改正；Δriono是电离层改正；ΔrSSB是海况(电磁)偏差改正；Δrset是固体潮改正；Δrpol是极潮改正。
6.如权利要求1所述的提取算法，其特征在于，所述提取算法还包括验证步骤。
CN201710617728.2A 2017-07-26 2017-07-26 一种高精度的水面高程提取方法 Active CN107421496B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710617728.2A CN107421496B (zh) 2017-07-26 2017-07-26 一种高精度的水面高程提取方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710617728.2A CN107421496B (zh) 2017-07-26 2017-07-26 一种高精度的水面高程提取方法

Publications (2)

Publication Number Publication Date
CN107421496A true CN107421496A (zh) 2017-12-01
CN107421496B CN107421496B (zh) 2018-09-04

Family

ID=60431367

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710617728.2A Active CN107421496B (zh) 2017-07-26 2017-07-26 一种高精度的水面高程提取方法

Country Status (1)

Country Link
CN (1) CN107421496B (zh)

Cited By (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108981658A (zh) * 2018-07-09 2018-12-11 中国科学院国家空间科学中心 一种基于星载干涉成像高度计的河流水面高程提取方法
CN109186561A (zh) * 2018-09-19 2019-01-11 南京大学 一种冰下湖体积变化的估算方法
CN109451866A (zh) * 2018-02-13 2019-03-08 北京小米移动软件有限公司 信息配置方法及装置、基站和用户设备
CN109556676A (zh) * 2018-03-02 2019-04-02 清华大学 河流水位的确定方法、装置、计算机设备及可读存储介质
CN109839623A (zh) * 2019-02-14 2019-06-04 北京遥感设备研究所 一种地外天体着陆测量雷达面目标回波信号测距处理方法
CN109839635A (zh) * 2019-03-13 2019-06-04 武汉大学 一种通过CryoSat-2 SARIn模式L1b级波形数据提取测高脚点高程的方法
CN110132373A (zh) * 2019-05-10 2019-08-16 宁波市测绘设计研究院 一种基于卫星测高数据的大型内陆湖泊和湿地水位反演方法
CN110378540A (zh) * 2019-08-02 2019-10-25 桂林理工大学 一种适用于广西北部湾地区的大气加权平均温度计算方法
CN110763302A (zh) * 2019-11-20 2020-02-07 北京航空航天大学 一种基于迭代频率估计的fmcw高精度液位测量方法
CN111812641A (zh) * 2020-07-22 2020-10-23 中国科学院空天信息创新研究院 一种基于多尺度峰值监测的雷达高度计波形重跟踪方法
CN112161670A (zh) * 2020-09-22 2021-01-01 长江水利委员会水文局长江上游水文水资源勘测局 一种顾及时间序列和空间结构河道水位改正方法
CN112414510A (zh) * 2020-10-10 2021-02-26 武汉大学 卫星测高河流水位监测方法及系统
CN114756640A (zh) * 2022-04-27 2022-07-15 国家卫星海洋应用中心 一种海面高度数据的评价方法及装置
CN115422981A (zh) * 2022-11-04 2022-12-02 自然资源部第一海洋研究所 面向单频机载激光测深数据的水陆分类方法、系统及应用
CN116719025A (zh) * 2023-08-09 2023-09-08 中国科学院地理科学与资源研究所 一种内陆水面高程测量方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2006258436A (ja) * 2005-03-15 2006-09-28 Japan Radio Co Ltd 衛星航法装置
CN103323840A (zh) * 2012-03-22 2013-09-25 中国科学院电子学研究所 干涉sar回波数据与平台运动及姿态数据的时间对准方法
CN105738881A (zh) * 2016-03-09 2016-07-06 中国石油大学(华东) 基于波形分类的近海Envisat卫星测高波形重定方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2006258436A (ja) * 2005-03-15 2006-09-28 Japan Radio Co Ltd 衛星航法装置
CN103323840A (zh) * 2012-03-22 2013-09-25 中国科学院电子学研究所 干涉sar回波数据与平台运动及姿态数据的时间对准方法
CN105738881A (zh) * 2016-03-09 2016-07-06 中国石油大学(华东) 基于波形分类的近海Envisat卫星测高波形重定方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
JINGJUAN LIAO; LE GAO; XIAOMING WANG: "Numerical Simulation and Forecasting of Water Level for Qinghai Lake Using Multi-Altimeter Data Between 2002 and 2012", 《JOURNAL OF SELECTED TOPICS IN APPLIED EARTH OBSERVATIONS AND REMOTE SENSING》 *
R. A. CULLEN; D. J. WINGHAM: "CryoSat Level lb Processing Algorithms and Simulation Results", 《IEEE INTERNATIONAL GEOSCIENCE & REMOTE SENSING SYMPOSIUM》 *
汪海洪等: "基于波形分类的近海卫星测高数据自适应重跟踪方法", 《测绘学报》 *

Cited By (24)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109451866A (zh) * 2018-02-13 2019-03-08 北京小米移动软件有限公司 信息配置方法及装置、基站和用户设备
CN109451866B (zh) * 2018-02-13 2023-06-16 北京小米移动软件有限公司 信息配置方法及装置、基站和用户设备
CN109556676A (zh) * 2018-03-02 2019-04-02 清华大学 河流水位的确定方法、装置、计算机设备及可读存储介质
CN108981658A (zh) * 2018-07-09 2018-12-11 中国科学院国家空间科学中心 一种基于星载干涉成像高度计的河流水面高程提取方法
CN108981658B (zh) * 2018-07-09 2019-11-19 中国科学院国家空间科学中心 一种基于星载干涉成像高度计的河流水面高程提取方法
CN109186561B (zh) * 2018-09-19 2020-10-02 南京大学 一种冰下湖体积变化的估算方法
CN109186561A (zh) * 2018-09-19 2019-01-11 南京大学 一种冰下湖体积变化的估算方法
CN109839623A (zh) * 2019-02-14 2019-06-04 北京遥感设备研究所 一种地外天体着陆测量雷达面目标回波信号测距处理方法
CN109839623B (zh) * 2019-02-14 2020-09-11 北京遥感设备研究所 一种地外天体着陆测量雷达面目标回波信号测距处理方法
CN109839635A (zh) * 2019-03-13 2019-06-04 武汉大学 一种通过CryoSat-2 SARIn模式L1b级波形数据提取测高脚点高程的方法
CN109839635B (zh) * 2019-03-13 2022-12-27 武汉大学 一种通过CryoSat-2 SARIn模式L1b级波形数据提取测高脚点高程的方法
CN110132373A (zh) * 2019-05-10 2019-08-16 宁波市测绘设计研究院 一种基于卫星测高数据的大型内陆湖泊和湿地水位反演方法
CN110378540B (zh) * 2019-08-02 2023-05-09 桂林理工大学 一种适用于广西北部湾地区的大气加权平均温度计算方法
CN110378540A (zh) * 2019-08-02 2019-10-25 桂林理工大学 一种适用于广西北部湾地区的大气加权平均温度计算方法
CN110763302A (zh) * 2019-11-20 2020-02-07 北京航空航天大学 一种基于迭代频率估计的fmcw高精度液位测量方法
CN111812641A (zh) * 2020-07-22 2020-10-23 中国科学院空天信息创新研究院 一种基于多尺度峰值监测的雷达高度计波形重跟踪方法
CN111812641B (zh) * 2020-07-22 2023-03-14 中国科学院空天信息创新研究院 一种基于多尺度峰值监测的雷达高度计波形重跟踪方法
CN112161670A (zh) * 2020-09-22 2021-01-01 长江水利委员会水文局长江上游水文水资源勘测局 一种顾及时间序列和空间结构河道水位改正方法
CN112414510A (zh) * 2020-10-10 2021-02-26 武汉大学 卫星测高河流水位监测方法及系统
CN114756640A (zh) * 2022-04-27 2022-07-15 国家卫星海洋应用中心 一种海面高度数据的评价方法及装置
CN114756640B (zh) * 2022-04-27 2022-10-21 国家卫星海洋应用中心 一种海面高度数据的评价方法及装置
CN115422981A (zh) * 2022-11-04 2022-12-02 自然资源部第一海洋研究所 面向单频机载激光测深数据的水陆分类方法、系统及应用
CN116719025A (zh) * 2023-08-09 2023-09-08 中国科学院地理科学与资源研究所 一种内陆水面高程测量方法
CN116719025B (zh) * 2023-08-09 2023-10-13 中国科学院地理科学与资源研究所 一种内陆水面高程测量方法

Also Published As

Publication number Publication date
CN107421496B (zh) 2018-09-04

Similar Documents

Publication Publication Date Title
CN107421496B (zh) 一种高精度的水面高程提取方法
Biancamaria et al. Satellite radar altimetry water elevations performance over a 200 m wide river: Evaluation over the Garonne River
Birkinshaw et al. Using satellite altimetry data to augment flow estimation techniques on the Mekong River
Nilsson et al. Improved retrieval of land ice topography from CryoSat-2 data and its impact for volume-change estimation of the Greenland Ice Sheet
Haase et al. Accuracy and variability of GPS tropospheric delay measurements of water vapor in the western Mediterranean
Cretaux et al. Hydrological applications of satellite altimetryrivers, lakes, man-made reservoirs, inundated areas
Crétaux et al. An absolute calibration site for radar altimeters in the continental domain: Lake Issykkul in Central Asia
Stanislav et al. Comparing wave heights simulated in the Black Sea by the SWAN model with satellite data and direct wave measurements
Frappart et al. Interannual variations of the terrestrial water storage in the Lower Ob'Basin from a multisatellite approach
Semmling et al. On the retrieval of the specular reflection in GNSS carrier observations for ocean altimetry
Den Ouden et al. Stand-alone single-frequency GPS ice velocity observations on Nordenskiöldbreen, Svalbard
Rypina et al. Eulerian and Lagrangian correspondence of high-frequency radar and surface drifter data: Effects of radar resolution and flow components
CN106768179A (zh) 基于连续运行gnss站信噪比数据的潮位的测量方法
CN106767383A (zh) 基于连续运行gnss站信噪比数据的积雪深度的测量方法
Hussain et al. Spatial and temporal variations of terrestrial water storage in upper Indus basin using GRACE and altimetry data
Guo et al. Sea level changes of China seas and neighboring ocean based on satellite altimetry missions from 1993 to 2012
Powell et al. An optimal filter for geostrophic mesoscale currents from along-track satellite altimetry
Din et al. The impact of sea level rise on geodetic vertical datum of Peninsular Malaysia
Xu et al. Lake level changes determined by Cryosat-2 altimetry data and water-induced loading deformation around Lake Qinghai
Lee Radar altimetry methods for solid earth geodynamics studies
Kang et al. Use of GNSS-derived PWV for predicting the path of typhoon: case studies of Soulik and Kongrey in 2018
Wang Ocean tide modeling in the Southern Ocean
Roohi Capability of pulse-limited satellite radar altimetry to monitor inland water bodies
Knudsen et al. A global mean ocean circulation estimation using goce gravity models-the DTU12MDT mean dynamic topography model
Zulkifli et al. Vertical land motion quantification using space-based geodetic methods: A review

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant