CN115201879A - 基于北斗geo卫星反射信号的洪涝监测方法 - Google Patents

基于北斗geo卫星反射信号的洪涝监测方法 Download PDF

Info

Publication number
CN115201879A
CN115201879A CN202210887952.4A CN202210887952A CN115201879A CN 115201879 A CN115201879 A CN 115201879A CN 202210887952 A CN202210887952 A CN 202210887952A CN 115201879 A CN115201879 A CN 115201879A
Authority
CN
China
Prior art keywords
signal
monitoring
noise ratio
formula
flood
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
Application number
CN202210887952.4A
Other languages
English (en)
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.)
Hangzhou Dianzi University
Original Assignee
Hangzhou Dianzi University
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 Hangzhou Dianzi University filed Critical Hangzhou Dianzi University
Priority to CN202210887952.4A priority Critical patent/CN115201879A/zh
Publication of CN115201879A publication Critical patent/CN115201879A/zh
Pending legal-status Critical Current

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
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/42—Determining position
    • 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/24—Acquisition or tracking or demodulation of signals transmitted by the system
    • G01S19/25—Acquisition or tracking or demodulation of signals transmitted by the system involving aiding data received from a cooperating element, e.g. assisted GPS
    • G01S19/258—Acquisition or tracking or demodulation of signals transmitted by the system involving aiding data received from a cooperating element, e.g. assisted GPS relating to the satellite constellation, e.g. almanac, ephemeris data, lists of satellites in view
    • 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/24—Acquisition or tracking or demodulation of signals transmitted by the system
    • G01S19/29—Acquisition or tracking or demodulation of signals transmitted by the system carrier including Doppler, related
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A10/00—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE at coastal zones; at river basins
    • Y02A10/40—Controlling or monitoring, e.g. of flood or hurricane; Forecasting, e.g. risk assessment or mapping

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明公开了一种基于北斗GEO卫星反射信号的洪涝监测方法，包括如下步骤：S1、读取监测日北斗GEO卫星原始观测值；S2、对原始信噪比进行数据处理得到反射信号信噪比；S3、通过反射信号信噪比获取反射信号信噪比幅度值；S4、采用分卫星分频率的方法建立北斗GEO卫星基于第一监测量和第二监测量的参考监测模型；S5、联合北斗GEO三频信号的第一监测量和第二监测量通过参考监测模型对洪涝进行监测，并根据监测阈值将结果进行标记。该方法不仅消除了直射信号的影响，同时反射信号信噪比幅度值受环境影响小，且北斗GEO卫星属于静止卫星，相应的反射信号信噪比时间序列特征更加明显，因此可以有效提高洪涝监测的精确性和实时性。

Description

基于北斗GEO卫星反射信号的洪涝监测方法
技术领域
本发明涉及导航卫星遥感反演技术领域，具体指一种基于北斗GEO（Geostationary Orbit，地球静止轨道卫星）卫星反射信号信噪比和幅度值的洪涝监测方法。
背景技术
洪涝灾害有着发生速度快、影响范围广和灾害量极大等显著特点，特别是城乡环境下的短时洪涝灾害位于人口集中地区，严重影响着人民的生命和财产安全，是威胁人类生存和发展的主要自然灾害之一。据统计显示，近二十年来，全球范围受洪涝灾害影响的总人口数超过10亿人，经济损失高达3万亿美元。在我国，由于水文和气候分布具有显著特点，南方地区雨季会导致洪涝灾害频繁发生，而北方地区的极端气候也极易引起洪涝灾害。比如，2021年在北方就发生过三次极端强降雨导致的城市洪涝灾害，已造成300多万人次受灾，经济损失高达2000多亿元。因此，如何建立科学、有效和高精度的洪涝监测方法，提高洪涝灾害监测的实时性，保障人民的生命和财产安全，是当前洪涝灾害监测管理的迫切需求，也是洪涝灾害预警和防范领域亟待解决的问题之一。
当前，洪涝监测方法主要可以分为三大类，第一类是基于传统雨量测量计的监测方法，其测量精度较高，单个测站成本低。然而，该方法测量范围小，导致监测的时间和空间分辨率较低，且难以大范围布站。第二类是基于遥感成像卫星的洪涝监测，该方法测量精度较高，监测范围广。然而，遥感成像卫星受天气环境影响较大，恶劣天气环境会导致监测方法失效或精度降低，且监测数据需要回传地面处理，导致实时性较低。第三类是基于导航系统卫星的洪涝监测。该方法通过大地测量型接收机接收导航卫星信号进行反演测量，测量精度高，且测量范围广。然而，现有基于导航卫星的洪涝监测是利用信号原始信噪比直接建模，而原始信噪比包含的直射信号会导致建模精度降低，且采用的相位值监测法受环境影响大，对测量地的土壤类型更加敏感，导致监测精度低等问题。此外，现有的监测方法往往基于IGSO（Inclined Geosynchronous Satellite Orbit）倾斜轨道卫星或MEO（MediumEarth Orbit）中轨道卫星进行监测，而这两类卫星轨道始终处于运动状态，卫星的信噪比误差较大，导致洪涝监测精度降低。
发明内容
本发明提出一种基于北斗GEO卫星反射信号的洪涝监测方法,不仅消除了直射信号的影响，同时反射信号信噪比幅度值受环境影响小，且北斗GEO卫星属于静止卫星，相应的反射信号信噪比时间序列特征更加明显，因此可以有效提高洪涝监测的精确性和实时性，也可以保证洪涝监测的成功率和稳定性，为洪涝灾害预警提供强有力的支撑。
为了解决上述技术问题，本发明的技术方案为：
一种基于北斗GEO卫星反射信号的洪涝监测方法，包括如下步骤：
S1、读取监测日北斗GEO卫星原始观测值，并计算相应历元时刻的高度角信息和时间信息，所述原始观测值包括原始信噪比、原始伪距和原始载波；
S2、对原始信噪比进行数据处理得到反射信号信噪比，并作为第一监测量；
S3、通过反射信号信噪比获取反射信号信噪比幅度值，并作为第二监测量；
S4、采用分卫星分频率的方法建立北斗GEO卫星基于第一监测量和第二监测量的参考监测模型,所述参考监测模型是以数据库文本的形式将历年相应历元时刻的第一监测量和第二监测量存储在本地电脑中，另外还存储相应历元时刻的高度角和时间信息，同时也将相应的监测阈值进行存储；
S5、联合北斗GEO三频信号的第一监测量和第二监测量通过参考监测模型对洪涝进行监测，并根据监测阈值将结果进行标记。
作为优选，所述步骤S2中，对原始信噪比进行数据处理的方法为：
S2-1、采用基于高斯拟合方法消除直射信号信噪比的影响；
S2-2、采用基于经验模态分解去噪算法滤除高频随机噪声，得到仅保留含有反射信号信噪比的信号。
作为优选，所述步骤S2-1的具体方法如下：
S2-1-1、利用高斯拟合算法对原始信噪比的观测值进行拟合，得到高斯拟合函数，
其中高斯拟合函数求解过程表示如下：假设原始信噪比观测数据集为（x i,S i）(i=1,2,3,....,N)，N表示信号长度，x表示信号的时间序列，S表示该时间下的原始信噪比，该数据集可以用高斯函数进行描述，描述公式为：
Figure 759350DEST_PATH_IMAGE001
（1）
式中，S i 和x i是已知的，exp表示指数运算，数学领域常识计算，a，b和c为未知参数，分别表示高斯函数曲线的峰高、峰的位置和半宽度信息，
a，b和c这三个参数可以通过下列过程求解：将式（1）两边取自然对数，可以将式（1）转换为：
Figure 87563DEST_PATH_IMAGE002
（2）
令：
Figure 961978DEST_PATH_IMAGE003
，
Figure 401050DEST_PATH_IMAGE004
,
Figure 677310DEST_PATH_IMAGE005
,
Figure 179354DEST_PATH_IMAGE006
，则式（2）可以化为二次多项式拟合函数，表示如下：
Figure 541066DEST_PATH_IMAGE007
（3）
考虑全部数据和测量误差，并以矩阵形式表示如下：
Figure 518249DEST_PATH_IMAGE008
（4）
上式可以化简为：
Figure 180174DEST_PATH_IMAGE009
（5）
式中，Z为已知量，由
Figure 584611DEST_PATH_IMAGE010
计算得到，X为已知量，由x i得到，因此，在不考虑总体测量误差E影响的情况下，根据最小二乘原理可以求得拟合常数b 0， b 1， b 2构成的矩阵B的广义最小二乘解，表示为：
Figure 433618DEST_PATH_IMAGE011
（6）
式中，X和Z都是已知量，T表示矩阵转置运算，通过式（6）即可求出拟合常数b 0， b 1，b 2，进一步可以通过下列式子求出未知参数a, b, c：
Figure 480072DEST_PATH_IMAGE012
（7）
式（7）中b 0， b 1， b 2已通过式（6）求出，e为自然常数，因此，通过式（7）即可求出未知参数a， b， c，进而可以得到由原始信噪比数据拟合的高斯函数，也即式（1）；
S2-1-2、利用上步求到的高斯拟合函数计算每一个历元时刻的拟合值，表示为S G ，计算过程可以表示为：
Figure 262083DEST_PATH_IMAGE013
（8）
式（8）中a， b， c已经通过式（7）求出，而x i 代表时间序列，因此可以求出对应时刻的拟合值S G ；
S2-1-3、利用原始信噪比S d+r+n 减去上步计算出来的拟合值S G ，得到仅包含反射信号信噪比和随机噪声的信号，公式表示为：
S r+n =S d+r+n -S G （9）
通过上述过程可以从原始信噪比观测值中消除直射信号信噪比影响，得到仅包含反射信号信噪比和随机噪声的信号S r+n 。
作为优选，所述步骤S2-2的具体方法如下：
S2-2-1、信号分解
对给定信号S r+n (x)进行处理，找到该信号上所有局部极大值和极小值所处的位置，然后将所有的极大值点和极小值点通过插值函数生成给定信号S r+n (x)的上下两条包络线，上包络线由极大值插值生成，表示为a 1，下包络线由极小值插值生成,表示为b 1，
将得到的两条包络线求均值，表示为：
Figure 837421DEST_PATH_IMAGE014
（10）
式中a 1和b 1即通过插值函数生成的上下包络线，然后用给定信号S r+n (x)减去式（10）中得到的
Figure 173724DEST_PATH_IMAGE015
，得到第一组残差
Figure 492710DEST_PATH_IMAGE016
，表示为：
Figure 863648DEST_PATH_IMAGE017
（11）
判断式（11）中计算得到的
Figure 875467DEST_PATH_IMAGE018
是否符合本征模态函数的两个假设条件，若不符合，将
Figure 699066DEST_PATH_IMAGE019
作为新的输入信号重复上述分解过程，直到分解出来的
Figure 556164DEST_PATH_IMAGE020
符合本征模态函数的两个假设条件为止，
假设分解到第k次后
Figure 47188DEST_PATH_IMAGE021
符合本征模态函数（IMF）的两个假设条件，第k次分解表示为：
Figure 964328DEST_PATH_IMAGE022
（12）
则此时的
Figure 275224DEST_PATH_IMAGE023
便是S r+n (x)分解出来的第一个本征模态函数，表示为IMF1，它的周期是后续所有分解出来的本征模态函数中最短的，频率最高的，
用原始信号减去第一个本征模态函数IMF1，得到去除IMF1的新的数据序列r 1，对这个新序列r 1继续重复进行IMF1的求取过程，得到第二个本征模态函数IMF2，接着用r 1减去第二个本征模态函数IMF2，得到新的数据序列r 2；如此重复计算，直到最后的r m 无法再继续分解为止，上述过程表示为：
Figure 936013DEST_PATH_IMAGE024
（13）
此时，r m 便是原始信号S r+n (x)经过经验模态分解后剩下的残差信号，表示为：Res，
S2-2-2、信号重构
经过上述信号经验模态分解后，每一个本征模态函数IMF都表示原始信号中的一个独立频率成分，频率由高到低依次排列，因此，可以用分离出来的所有的IMF与残差Res对信号进行重构，得到去噪后的信号，重构公式表示如下：
Figure 281543DEST_PATH_IMAGE025
（14）
通过上述过程处理后，可以有效滤除高频随机噪声，得到仅保留反射信号信噪比的信号,即
Figure 369585DEST_PATH_IMAGE026
，表示为
Figure 170706DEST_PATH_IMAGE027
，并以此信号作为第一监测量。
作为优选，所述步骤S3中，反射信号信噪比幅度值的获取方法：
S3-1、利用信噪比线性-周期转换模型将以dB-Hz为单位的线性反射信号信噪比转换为以volts/volts为单位的周期反射信号信噪比，线性-周期转换模型的转换公式表示如下：
Figure 635186DEST_PATH_IMAGE028
（15）
式中，
Figure 835223DEST_PATH_IMAGE029
为仅包含反射信号信噪比的信号，该信号是以dB-Hz为单位的线性信号，S r_v 为经过线性-周期转换模型后得到的以volts/volts为单位的周期反射信号信噪比；
S3-2、采用Lomb-Scargle谱分析法求卫星有效反射高
Figure 94166DEST_PATH_IMAGE030
；
S3-3、采用高斯粒子滤波算法求出反射信号信噪比幅度值，具体公式表示如下：
Figure 114075DEST_PATH_IMAGE031
式中，A和
Figure 647824DEST_PATH_IMAGE032
表示该信噪比的幅度和相位，是公式中的未知量，S r_v 是周期反射信号信噪比，
Figure 436788DEST_PATH_IMAGE033
是卫星有效反射高，λ是信号的波长，属于已知量，EL是对应的高度角。
作为优选，所述步骤S3-2中卫星有效反射高
Figure 132212DEST_PATH_IMAGE034
的求解方法：
S3-2-1、通过Lomb-Scargle谱分析法求周期反射信号信噪比S r_v 的功率最大值时的频率
Figure 904996DEST_PATH_IMAGE035
，具体过程表示如下：
采用Lomb-Sacragle谱分析法求信号的功率谱，公式表示如下：
Figure 711278DEST_PATH_IMAGE036
(17)
式中，P x ( f )是频率为 f 的周期信号的功率；S r_v (x)为反射信号信噪比，N为序列长度，τ为时间平移不变量，可通过下式计算得到：
Figure 354749DEST_PATH_IMAGE037
通过式（17-18），可求出P x ( f )达到峰值时对应的
Figure 486653DEST_PATH_IMAGE038
；
S3-2-2、利用
Figure 215575DEST_PATH_IMAGE039
和卫星有效反射高
Figure 825547DEST_PATH_IMAGE040
的关系求出卫星有效反射高，反射信号信噪比频率
Figure 589104DEST_PATH_IMAGE041
和卫星有效反射高的关系表示如下：
Figure 626330DEST_PATH_IMAGE042
式中，λ是信号的波长，属于已知量，通过上述过程可求出卫星的有效反射高
Figure 108127DEST_PATH_IMAGE043
。
作为优选，所述步骤S4中监测阈值为第一监测量的波动不超过2dB-Hz，第二监测量的波动不超过0.5。
作为优选，所述步骤S5包括如下步骤：
S5-1、采用基于北斗GEO卫星三频反射信号信噪比对洪涝进行监测，并根据监测结果进行标记，所述标记的级别包括一级标记和二级标记；
S5-2、采用基于北斗GEO卫星三频反射信号信噪比的幅度值对洪涝进行监测，并根据监测结果进行标记，所述标记的级别包括一级标记和二级标记；
S5-3、将经过第一监测量和第二监测量的洪涝监测结果进行联合处理，若第一监测量和第二监测量存在一级标记，则认为存在洪涝危险，需进行一级洪涝预警处理，若第一监测量和第二监测量无一级标记但有二级标记，则进行二级洪涝预警，若两者都没有标记，不进行预警。
作为优选，所述步骤S5-1中对洪涝进行监测的方法为：
首先，获取当前历元时刻的原始信噪比通过步骤S2计算得到当前历元时刻的第一监测量，并计算该历元时刻的高度角信息和时间信息；
其次，将计算得到的北斗GEO卫星三频，即B1频段、B2频段和B3频段的高度角信息和时间信息，并输入参考监测模型中搜索相应参考值中的第一监测量，并将计算得到的当前历元时刻的第一监测量与参考值中的第一监测量作对比，若当前历元时刻的第一监测量波动大于2dB-Hz时存在洪涝，则进行标记；
最后，根据历元标记结果进行统计分析，若标记历元占总历元数的80%以上，则认为存在强洪涝相关，做一级标记，若标记历元占总历元数的60%-80%内，则认为存在较强洪涝相关，做二级标记，若标记历元占比低于60%，则不进行标记，最后，将监测结果进行保存并输出。
作为优选，所述步骤S5-2中对洪涝进行监测的方法为：
首先，根据步骤S5-1计算得到的当前历元的反射信号信噪比通过步骤S3计算得到反射信号信噪比的幅度值，即当前历元时刻的第二检测量，并计算该历元时刻的高度角和时间信息；
其次，将计算得到的北斗GEO卫星三频即B1频段、B2频段和B3频段的高度角信息和时间信息，并输入参考监测模型中搜索相应参考值中的第二监测量，并将计算得到的当前历元时刻的第二监测量与参考值中的第二监测量作对比，若当前历元时刻的反射信号信噪比幅度值波动大于0.5时存在洪涝，则进行标记；
最后，根据历元标记结果进行统计分析，若标记历元占总历元数的80%以上，则认为存在强洪涝相关，做一级标记，若标记历元占总历元数的60%-80%内，则认为存在较强洪涝相关，做二级标记，若标记历元占比低于60%，则不进行标记，最后，将监测结果进行保存并输出。
本发明具有以下的特点和有益效果：
采用上述技术方案，可以利用现有北斗导航GEO卫星对洪涝进行监测，与现有基于传统雨量计测量方法相比，可以解决监测范围小等问题。与现在基于遥感成像卫星测量方法相比，可以解决受环境影响大、实时性低等问题。与现有基于导航系统卫星测量方法相比，既可以解决基于IGSO和MEO卫星因轨道运动导致的监测精度低等问题，又可以解决因忽略直射信号信噪比影响和忽略反射信号信噪比幅度值信息而导致监测精度低和性能不稳定等问题。同时，联合北斗GEO三频信号联合监测可以保证监测成功率和稳定性等问题，从而满足高精度、高实时性和高稳定性的洪涝监测。考虑到直射信号信噪比导致的监测精度低等问题，采用基于高斯拟合方法对原始信噪比观测值中的直射信号信噪比进行剔除。利用基于经验模态分解去噪算法滤除高频随机噪声，提取高精度的反射信号信噪比，作为第一检测量。利用信噪比线性-周期转换模型将以dB-Hz为单位的线性反射信号信噪比转换为以volts/volts为单位的周期反射信号信噪比信号。采用基于Lomb-Scargle谱分析法计算卫星有效反射高，消除高度角非均匀特征影响。利用高斯粒子滤波算法计算反射信号信噪比幅度值，以此作为洪涝监测的第二检测量，提高监测精度。联合北斗GEO卫星B1频段、B2频段和B3频段进行联合监测，进一步保障洪涝监测的成功率和准确率。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1 基于北斗GEO卫星反射信号的洪涝监测方法流程图。
图2 基于北斗GEO卫星反射信号信噪比的洪涝监测方法流程图。
图3 基于北斗GEO卫星反射信号信噪比幅度值的洪涝监测方法流程图。
具体实施方式
需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。
在本发明的描述中，需要理解的是，术语“中心”、“纵向”、“横向”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”等仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”等的特征可以明示或者隐含地包括一个或者更多个该特征。在本发明的描述中，除非另有说明，“多个”的含义是两个或两个以上。
在本发明的描述中，需要说明的是，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以通过具体情况理解上述术语在本发明中的具体含义。
本发明提供了一种基于北斗GEO卫星反射信号的洪涝监测方法，如图1所示，包括如下步骤：
S1、读取监测日北斗GEO卫星原始观测值，并计算相应历元时刻的高度角信息和时间信息，所述原始观测值包括原始信噪比、原始伪距和原始载波。
具体的，原始观测值的获取是采用普通大地测量型接收机连续收集非洪涝日的北斗GEO卫星原始数据，包括信噪比、伪距和载波相位观测值，卫星星历等观测数据。同时计算相应历元时刻的高度角，方位角以及历元时间等信息，为后续建立洪涝监测模型提供辅助参数信息。
需要说明的是，信噪比可以直接从观测文件中提取，相应历元时刻的高度角和时间信息可以通过伪距、载波相位和卫星星历计算得到。该过程属于定位领域常识，这里不再详细列举计算过程。为方便后续陈述，这里提取的原始信噪比在后续步骤中用S d+r+n 表示，计算出来的高度角用EL表示。
S2、对原始信噪比进行数据处理得到反射信号信噪比，并作为第一监测量。
可以理解的，提取的原始信噪比S d+r+n 中包括直射信号信噪比（表示为S d ）、反射信号信噪比（表示为S r ）、以及随机噪声（表示为S n ），因此需要对原始信噪比进行数据处理，具体的处理方法如下：
S2-1、采用基于高斯拟合方法消除直射信号信噪比的影响。
S2-1-1、利用高斯拟合算法对原始信噪比的观测值进行拟合，得到高斯拟合函数。
其中，高斯拟合函数求解过程表示如下：假设原始信噪比观测数据集为（x i,S i）(i=1,2,3,....,N)，N表示信号长度，x表示信号的时间序列，S表示该时间下的原始信噪比，该数据集可以用高斯函数进行描述，描述公式为：
Figure 256212DEST_PATH_IMAGE044
（1）
式中，S i 和x i是已知的，exp表示指数运算，数学领域常识计算，a，b和c为未知参数，分别表示高斯函数曲线的峰高、峰的位置和半宽度信息，
a，b和c这三个参数可以通过下列过程求解：将式（1）两边取自然对数，可以将式（1）转换为：
Figure 139854DEST_PATH_IMAGE045
（2）
令：
Figure 347982DEST_PATH_IMAGE046
，
Figure 317075DEST_PATH_IMAGE047
,
Figure 268850DEST_PATH_IMAGE048
,
Figure 9929DEST_PATH_IMAGE049
，则式（2）可以化为二次多项式拟合函数，表示如下：
Figure 654537DEST_PATH_IMAGE050
（3）
考虑全部数据和测量误差，并以矩阵形式表示如下：
Figure 110926DEST_PATH_IMAGE051
（4）
上式可以化简为：
Figure 600813DEST_PATH_IMAGE052
（5）
式中，Z为已知量，由
Figure 459047DEST_PATH_IMAGE053
计算得到，X为已知量，由x i得到，因此，在不考虑总体测量误差E影响的情况下，根据最小二乘原理可以求得拟合常数b 0，b 1，b 2构成的矩阵B的广义最小二乘解，表示为：
Figure 743398DEST_PATH_IMAGE054
（6）
式中，X和Z都是已知量，T表示矩阵转置运算，通过式（6）即可求出拟合常数b 0，b 1，b 2，进一步可以通过下列式子求出未知参数a, b, c：
Figure 952663DEST_PATH_IMAGE055
（7）
式（7）中b 0， b 1， b 2已通过式（6）求出，e为自然常数，属于数学常识，这里不在解释。因此，通过式（7）即可求出未知参数a， b， c，进而可以得到由原始信噪比数据拟合的高斯函数，也即式（1）；
S2-1-2、利用上步求到的高斯拟合函数计算每一个历元时刻的拟合值，表示为S G ，计算过程可以表示为：
Figure 246241DEST_PATH_IMAGE056
（8）
式（8）中a， b， c已经通过式（7）求出，而x i 代表时间序列，因此可以求出对应时刻的拟合值S G ；
S2-1-3、利用原始信噪比S d+r+n 减去上步计算出来的拟合值S G ，得到仅包含反射信号信噪比和随机噪声的信号，公式表示为：
S r+n =S d+r+n -S G （9）
通过上述过程可以从原始信噪比观测值中消除直射信号信噪比影响，得到仅包含反射信号信噪比和随机噪声的信号S r+n 。
S2-2、采用基于经验模态分解去噪算法滤除高频随机噪声，得到仅保留含有反射信号信噪比的信号。
可以理解的，经验模态分解方法本质上是对一个含噪信号进行平稳化处理，将信号循环分解为本征模态函数(Intrinsic Model Function，表示为：IMF)和残差，直到分解出来的本征模态函数符合以下两个假设为止：1）本征模态函数绘制的曲线内过零点的个数与曲线上所有极值点个数必须满足差值不大于一的条件。2）本征模态函数绘制的曲线中，同一点上下包络线的值的和为零。达到上述两个条件后，则停止分解。并利用上述分解所得的本征模态函数对信号进行重构，从而达到去噪的效果。该方法能够有效滤除随机噪声，最大程度保留有用信号。基于经验模态分解的去噪算法主要有两部分组成：
S2-2-1、信号分解
对给定信号S r+n (x)进行处理，找到该信号上所有局部极大值和极小值所处的位置，然后将所有的极大值点和极小值点通过插值函数生成给定信号S r+n(x)的上下两条包络线，上包络线由极大值插值生成，表示为a 1，下包络线由极小值插值生成,表示为b 1，需要说明的是，插值函数可以使用拉格朗日插值法或三次样条插值法，这两种方法都是数学领域常识方法，这里不再描述。
进一步的，将得到的两条包络线求均值，表示为：
Figure 958982DEST_PATH_IMAGE057
（10）
式中a 1和b 1即通过插值函数生成的上下包络线，然后用给定信号S r+n (x)减去式（10）中得到的
Figure 679813DEST_PATH_IMAGE058
，得到第一组残差
Figure 845215DEST_PATH_IMAGE059
，表示为：
Figure 942484DEST_PATH_IMAGE060
（11）
判断式（11）中计算得到的
Figure 775311DEST_PATH_IMAGE061
是否符合本征模态函数的两个假设条件，若不符合，将
Figure 667044DEST_PATH_IMAGE062
作为新的输入信号重复上述分解过程，直到分解出来的
Figure 319742DEST_PATH_IMAGE063
符合本征模态函数的两个假设条件为止，
假设分解到第k次后
Figure 955123DEST_PATH_IMAGE064
符合本征模态函数（IMF）的两个假设条件，第k次分解表示为：
Figure 642456DEST_PATH_IMAGE065
（12）
则此时的
Figure 236248DEST_PATH_IMAGE066
便是S r+n (x)分解出来的第一个本征模态函数，表示为IMF1，它的周期是后续所有分解出来的本征模态函数中最短的，频率最高的，
用原始信号减去第一个本征模态函数IMF1，得到去除IMF1的新的数据序列r 1，对这个新序列r 1继续重复进行IMF1的求取过程，得到第二个本征模态函数IMF2，接着用r 1减去第二个本征模态函数IMF2，得到新的数据序列r 2；如此重复计算，直到最后的r m 无法再继续分解为止，上述过程表示为：
Figure 203391DEST_PATH_IMAGE067
（13）
此时，r m 便是原始信号S r+n(x)经过经验模态分解后剩下的残差信号，表示为：Res。
S2-2-2、信号重构
经过上述信号经验模态分解后，每一个本征模态函数IMF都表示原始信号中的一个独立频率成分，频率由高到低依次排列，因此，可以用分离出来的所有的IMF与残差Res对信号进行重构，得到去噪后的信号，重构公式表示如下：
Figure 642463DEST_PATH_IMAGE068
（14）
通过上述过程处理后，可以有效滤除高频随机噪声，得到仅保留反射信号信噪比的信号,即
Figure 184303DEST_PATH_IMAGE069
，表示为
Figure 417838DEST_PATH_IMAGE070
，并以此信号作为第一监测量。
S3、通过反射信号信噪比获取反射信号信噪比幅度值，并作为第二监测量。
对得到的高精度反射信号信噪比继续进行处理，利用信噪比线性-周期转换模型将以dB-Hz为单位的线性反射信号信噪比转换为以volts/volts为单位的周期反射信号信噪比信号。采用基于Lomb-Scargle谱分析法计算卫星有效反射高，消除高度角非均匀特征影响。利用高斯粒子滤波算法计算反射信号信噪比幅度值，以此作为洪涝监测的第二监测量。
具体的，反射信号信噪比幅度值的获取方法：
S3-1、利用信噪比线性-周期转换模型将以dB-Hz为单位的线性反射信号信噪比转换为以volts/volts为单位的周期反射信号信噪比，线性-周期转换模型的转换公式表示如下：
Figure 45128DEST_PATH_IMAGE071
（15）
式中，
Figure 22311DEST_PATH_IMAGE072
为仅包含反射信号信噪比的信号，该信号是以dB-Hz为单位的线性信号，
Figure 684237DEST_PATH_IMAGE073
为经过线性-周期转换模型后得到的以volts/volts为单位的周期反射信号信噪比；
S3-2、采用Lomb-Scargle谱分析法求卫星有效反射高
Figure 88673DEST_PATH_IMAGE074
。
该过程主要分为两部分，第一部分是通过Lomb-Scargle谱分析法求周期反射信号信噪比
Figure 203260DEST_PATH_IMAGE075
的功率最大值时的频率(表示为：
Figure 984134DEST_PATH_IMAGE076
)。第二部分是利用
Figure 234987DEST_PATH_IMAGE039
和卫星有效反射高(表示为：
Figure 810325DEST_PATH_IMAGE077
) 的关系求出卫星有效反射高。具体过程表示如下：
S3-2-1、通过Lomb-Scargle谱分析法求周期反射信号信噪比S r_v 的功率最大值时的频率
Figure 677787DEST_PATH_IMAGE038
，具体过程表示如下：
采用Lomb-Sacragle谱分析法求信号的功率谱，公式表示如下：
Figure 996773DEST_PATH_IMAGE078
(16)
式中，P x ( f )是频率为 f 的周期信号的功率；S r_v (x)为反射信号信噪比，N为序列长度，τ为时间平移不变量，可通过下式计算得到：
Figure 367711DEST_PATH_IMAGE079
通过式（16-17），可求出P x ( f )达到峰值时对应
Figure 379529DEST_PATH_IMAGE080
；
S3-2-2、利用
Figure 203129DEST_PATH_IMAGE080
和卫星有效反射高
Figure 325806DEST_PATH_IMAGE081
的关系求出卫星有效反射高，反射信号信噪比频率
Figure 554180DEST_PATH_IMAGE076
和卫星有效反射高的关系表示如下：
Figure 736900DEST_PATH_IMAGE082
式中，
Figure 47796DEST_PATH_IMAGE080
已通过公式（16-17）求出，λ是信号的波长，属于已知量，比如北斗B1频段的波长为：19.2cm，属于领域内的公知常识。通过上述过程可求出卫星的有效反射高。
S3-3、采用高斯粒子滤波算法求出反射信号信噪比幅度值，具体公式表示如下：
Figure 974163DEST_PATH_IMAGE083
式中，A和
Figure 319694DEST_PATH_IMAGE084
表示该信噪比的幅度和相位，是公式中的未知量，S r_v 是周期反射信号信噪比，
Figure 673315DEST_PATH_IMAGE085
是卫星有效反射高，λ是信号的波长，属于已知量，EL是对应的高度角。
可以理解的，上式可以通过最小二乘算法、卡尔曼滤波或高斯粒子滤波算法求解。
考虑到此时高度角呈非均匀分布，因此为了提高参数解算的精度，本实施例中，采用高斯粒子滤波算法进行求解。高斯粒子滤波算法结合了粒子滤波和高斯滤波的优点，是在粒子滤波的基本框架下用单个高斯分布取近似系统状态分布，从而得到最佳的参数解算。该算法是一种很成熟的算法，有很多开源代码可以在Matlab平台中实现。因此，这里仅作简单的算法流程说明，如下所示：
假设一个长度为N的m维高斯随机变量X的概率密度可以写作：
Figure 205927DEST_PATH_IMAGE086
式中，
Figure 935986DEST_PATH_IMAGE087
是信号均值，∑是信号的方差阵。T表示转置运算。exp为指数运算，这两种运算属于数学常识。假设在k-1时刻可以得到条件概率分布
Figure 401602DEST_PATH_IMAGE088
的近似如下：
Figure 926125DEST_PATH_IMAGE089
当输入新的观测值时，可以通过粒子滤波的预测步骤和更新步骤得到条件概率分布
Figure 946033DEST_PATH_IMAGE090
的近似如下：
Figure 214204DEST_PATH_IMAGE091
大多数情况下，上式中的均值
Figure 534326DEST_PATH_IMAGE092
和方差
Figure 229750DEST_PATH_IMAGE093
的解析表达式得不到，需要利用粒子滤波中的粒子及其权值来求
Figure 2534DEST_PATH_IMAGE094
和
Figure 543237DEST_PATH_IMAGE095
的估计值。因此，假设信号长度为N时，高斯粒子滤波算法的执行过程可以总结如下：
1）从
Figure 717866DEST_PATH_IMAGE096
中抽样得到样本集合
Figure 584191DEST_PATH_IMAGE097
；
2）分别从
Figure 847201DEST_PATH_IMAGE098
(i=1,2,...,N)中抽样，得到样本集合
Figure 191594DEST_PATH_IMAGE099
；
3）计算样本
Figure 220730DEST_PATH_IMAGE100
的均值
Figure 257956DEST_PATH_IMAGE101
和方差
Figure 739753DEST_PATH_IMAGE102
，表示如下：
Figure 153417DEST_PATH_IMAGE103
式中，T为转置运算，属于线性代数常识运算。
4）从重要性概率密度
Figure 771480DEST_PATH_IMAGE104
中抽取样本
Figure 245187DEST_PATH_IMAGE105
。通常
Figure 214280DEST_PATH_IMAGE106
可以采用
Figure 900476DEST_PATH_IMAGE107
或
Figure 904204DEST_PATH_IMAGE108
近似。
5）利用得到的观测值y k 计算各x k(i)的权值
Figure 283233DEST_PATH_IMAGE109
，计算公式表示如下：
Figure 739622DEST_PATH_IMAGE110
6）归一化权值
Figure 495089DEST_PATH_IMAGE111
得到
Figure 822165DEST_PATH_IMAGE112
，计算公式表示如下：
Figure 637674DEST_PATH_IMAGE113
7）利用
Figure 581359DEST_PATH_IMAGE114
计算均值
Figure 140517DEST_PATH_IMAGE115
和方差
Figure 587678DEST_PATH_IMAGE116
，计算公式表示如下：
Figure 311439DEST_PATH_IMAGE117
(26)
从而得到k时刻
Figure 476842DEST_PATH_IMAGE118
的近似如下：
Figure 574111DEST_PATH_IMAGE119
式中，均值
Figure 406937DEST_PATH_IMAGE120
可以看作为k时刻
Figure 33091DEST_PATH_IMAGE121
的估计值。
具体到本案例中，x就是包含未知变量的时间序列值，y就是对应的已知量。多维时x和y以矩阵的形式存在。因此，通过上述高斯粒子滤波算法，可以精确的求出公式（19）中的反射信号信噪比的幅度值A和相位值。本案例实施例中只用幅度值A。以此反射信号信噪比的幅度值A作为第二监测量。
S4、采用分卫星分频率的方法建立北斗GEO卫星基于第一监测量和第二监测量的参考监测模型,所述参考监测模型是以数据库文本的形式将历年相应历元时刻的第一监测量和第二监测量存储在本地电脑中，另外还存储相应历元时刻的高度角和时间信息，同时也将相应的监测阈值进行存储；监测阈值为第一监测量的波动不超过2dB-Hz，第二监测量的波动不超过0.5。
可以理解的，建立参考监测模型（模型以数据库的形式存放在本地电脑上）。利用计算得到的第一监测量和第二监测量，以及相应的辅助信息建立基于北斗GEO卫星反射信号信噪比和反射信号信噪比幅度值的洪涝监测模型。通过对北斗GEO卫星进行分卫星分频率分别建立第一监测量和第二监测量的参考监测模型，并按照相应历元的时间、高度角等信息进行存储，为后续洪涝监测奠定基础。具体的，本实施例中存储方式如下：北斗B01卫星B1频段12点12分30秒（按数据采样率为30秒进行举例）时刻的第一检测量为SB01_B1，第二检测量为SAB01_B1，对应的高度角为：ELB01_B1。时间信息如该例子中的12点12分30记录符号表示为：TB01_B1。
S5、联合北斗GEO三频信号的第一监测量和第二监测量通过参考监测模型对洪涝进行监测，并根据监测阈值将结果进行标记。具体的，如图2所示，
S5-1、采用基于北斗GEO卫星三频反射信号信噪比对洪涝进行监测，并根据监测结果进行标记，所述标记的级别包括一级标记和二级标记。
首先，获取当前历元时刻的原始信噪比通过步骤S2计算得到当前历元时刻的第一监测量，并计算该历元时刻的高度角信息和时间信息，为后续搜索参考模型提供基准。高度角计算过程较简单，属于该领域常规方法，此处不再描述。
其次，将计算得到的北斗GEO卫星三频，即B1频段、B2频段和B3频段的高度角信息和时间信息，并输入参考监测模型中搜索相应参考值中的第一监测量，并将计算得到的当前历元时刻的第一监测量与参考值中的第一监测量作对比，若当前历元时刻的第一监测量波动大于2dB-Hz时存在洪涝，则进行标记；
最后，根据历元标记结果进行统计分析，若标记历元占总历元数的80%以上，则认为存在强洪涝相关，做一级标记，若标记历元占总历元数的60%-80%内，则认为存在较强洪涝相关，做二级标记，若标记历元占比低于60%，则不进行标记，最后，将监测结果进行保存并输出。
S5-2、采用基于北斗GEO卫星三频反射信号信噪比的幅度值对洪涝进行监测，并根据监测结果进行标记，所述标记的级别包括一级标记和二级标记；具体的，如图3所示，
首先，根据步骤S5-1计算得到的当前历元的反射信号信噪比通过步骤S3计算得到反射信号信噪比的幅度值，即当前历元时刻的第二检测量，并计算该历元时刻的高度角和时间信息；
其次，将计算得到的北斗GEO卫星三频即B1频段、B2频段和B3频段的高度角信息和时间信息，并输入参考监测模型中搜索相应参考值中的第二监测量，并将计算得到的当前历元时刻的第二监测量与参考值中的第二监测量作对比，若当前历元时刻的反射信号信噪比幅度值波动大于0.5时存在洪涝，则进行标记；
最后，根据历元标记结果进行统计分析，若标记历元占总历元数的80%以上，则认为存在强洪涝相关，做一级标记，若标记历元占总历元数的60%-80%内，则认为存在较强洪涝相关，做二级标记，若标记历元占比低于60%，则不进行标记，最后，将监测结果进行保存并输出。
S5-3、将经过第一监测量和第二监测量的洪涝监测结果进行联合处理，若第一监测量和第二监测量存在一级标记，则认为存在洪涝危险，需进行一级洪涝预警处理，若第一监测量和第二监测量无一级标记但有二级标记，则进行二级洪涝预警，若两者都没有标记，不进行预警。
经过上述洪涝监测算法处理，可以有效解决现有基于IGSO和MEO卫星因轨道波动导致监测精度低等问题，也可以解决因忽略直射信号影响和反射信号信噪比幅度值作用等影响而导致监测精度低和监测性能不稳定等问题，为及时准确的洪涝监测提供保障。
以上结合附图对本发明的实施方式作了详细说明，但本发明不限于所描述的实施方式。对于本领域的技术人员而言，在不脱离本发明原理和精神的情况下，对这些实施方式包括部件进行多种变化、修改、替换和变型，仍落入本发明的保护范围内。

Claims (10)

1.一种基于北斗GEO卫星反射信号的洪涝监测方法，其特征在于，包括如下步骤：
S1、读取监测日北斗GEO卫星原始观测值，并计算相应历元时刻的高度角信息和时间信息，所述原始观测值包括原始信噪比、原始伪距和原始载波；
S2、对原始信噪比进行数据处理得到反射信号信噪比，并作为第一监测量；
S3、通过反射信号信噪比获取反射信号信噪比幅度值，并作为第二监测量；
S4、采用分卫星分频率的方法建立北斗GEO卫星基于第一监测量和第二监测量的参考监测模型,所述参考监测模型是以数据库文本的形式将历年相应历元时刻的第一监测量和第二监测量存储在本地电脑中，另外还存储相应历元时刻的高度角和时间信息，同时也将相应的监测阈值进行存储；
S5、联合北斗GEO三频信号的第一监测量和第二监测量通过参考监测模型对洪涝进行监测，并根据监测阈值将结果进行标记。
2.根据权利要求1所述的基于北斗GEO卫星反射信号的洪涝监测方法，其特征在于，所述步骤S2中，对原始信噪比进行数据处理的方法为：
S2-1、采用基于高斯拟合方法消除直射信号信噪比的影响；
S2-2、采用基于经验模态分解去噪算法滤除高频随机噪声，得到仅保留含有反射信号信噪比的信号。
3.根据权利要求2所述的基于北斗GEO卫星反射信号的洪涝监测方法，其特征在于，所述步骤S2-1的具体方法如下：
S2-1-1、利用高斯拟合算法对原始信噪比的观测值进行拟合，得到高斯拟合函数，
其中高斯拟合函数求解过程表示如下：假设原始信噪比观测数据集为（x i,S i）(i=1,2,3,....,N)，N表示信号长度，x表示信号的时间序列，S表示该时间下的原始信噪比，该数据集可以用高斯函数进行描述，描述公式为：
Figure 207264DEST_PATH_IMAGE001
（1）
式中，S i 和x i是已知的，exp表示指数运算，a，b和c为未知参数，分别表示高斯函数曲线的峰高、峰的位置和半宽度信息，
a，b和c这三个参数可以通过下列过程求解：将式（1）两边取自然对数，可以将式（1）转换为：
Figure 593246DEST_PATH_IMAGE002
（2）
令：
Figure 841825DEST_PATH_IMAGE003
，
Figure 971455DEST_PATH_IMAGE004
,
Figure 989089DEST_PATH_IMAGE005
,
Figure 545973DEST_PATH_IMAGE006
，则式（2）可以化为二次多项式拟合函数，表示如下：
Figure 281848DEST_PATH_IMAGE007
（3）
考虑全部数据和测量误差，并以矩阵形式表示如下：
Figure 684010DEST_PATH_IMAGE008
（4）
上式可以化简为：
Figure 87310DEST_PATH_IMAGE009
（5）
式中，Z为已知量，由
Figure 549515DEST_PATH_IMAGE010
计算得到，X为已知量，由x i得到，因此，在不考虑总体测量误差E影响的情况下，根据最小二乘原理可以求得拟合常数b 0， b 1， b 2构成的矩阵B的广义最小二乘解，表示为：
Figure 507107DEST_PATH_IMAGE011
（6）
式中，X和Z都是已知量，T表示矩阵转置运算，通过式（6）即可求出拟合常数b 0， b 1， b 2，进一步可以通过下列式子求出未知参数a, b, c：
Figure 244118DEST_PATH_IMAGE012
（7）
式（7）中b 0， b 1， b 2已通过式（6）求出，e为自然常数，因此，通过式（7）即可求出未知参数a， b， c，进而可以得到由原始信噪比数据拟合的高斯函数，也即式（1）；
S2-1-2、利用上步求到的高斯拟合函数计算每一个历元时刻的拟合值，表示为S G ，计算过程可以表示为：
Figure 236345DEST_PATH_IMAGE013
（8）
式（8）中a， b， c已经通过式（7）求出，而x i 代表时间序列，因此可以求出对应时刻的拟合值S G ；
S2-1-3、利用原始信噪比S d+r+n 减去上步计算出来的拟合值S G ，得到仅包含反射信号信噪比和随机噪声的信号，公式表示为：
S r+n =S d+r+n -S G （9）
通过上述过程可以从原始信噪比观测值中消除直射信号信噪比影响，得到仅包含反射信号信噪比和随机噪声的信号S r+n 。
4.根据权利要求3所述的基于北斗GEO卫星反射信号的洪涝监测方法，其特征在于，所述步骤S2-2的具体方法如下：
S2-2-1、信号分解
对给定信号S r+n (x)进行处理，找到该信号上所有局部极大值和极小值所处的位置，然后将所有的极大值点和极小值点通过插值函数生成给定信号S r+n (x)的上下两条包络线，上包络线由极大值插值生成，表示为a 1，下包络线由极小值插值生成,表示为b 1，
将得到的两条包络线求均值，表示为：
Figure 869452DEST_PATH_IMAGE014
（10）
式中a 1和b 1即通过插值函数生成的上下包络线，然后用给定信号S r+n (x)减去式（10）中得到的
Figure 845498DEST_PATH_IMAGE015
，得到第一组残差
Figure 852113DEST_PATH_IMAGE016
，表示为：
Figure 698846DEST_PATH_IMAGE017
（11）
判断式（11）中计算得到的
Figure 34012DEST_PATH_IMAGE018
是否符合本征模态函数的两个假设条件，若不符合，将
Figure 700617DEST_PATH_IMAGE019
作为新的输入信号重复上述分解过程，直到分解出来的
Figure 779432DEST_PATH_IMAGE020
符合本征模态函数的两个假设条件为止，
假设分解到第k次后
Figure 746251DEST_PATH_IMAGE021
符合本征模态函数（IMF）的两个假设条件，第k次分解表示为：
Figure 721160DEST_PATH_IMAGE022
（12）
则此时的
Figure 406219DEST_PATH_IMAGE023
便是S r+n (x)分解出来的第一个本征模态函数，表示为IMF1，它的周期是后续所有分解出来的本征模态函数中最短的，频率最高的，
用原始信号减去第一个本征模态函数IMF1，得到去除IMF1的新的数据序列r 1，对这个新序列r 1继续重复进行IMF1的求取过程，得到第二个本征模态函数IMF2，接着用r 1减去第二个本征模态函数IMF2，得到新的数据序列r 2；如此重复计算，直到最后的r m 无法再继续分解为止，上述过程表示为：
Figure 757566DEST_PATH_IMAGE024
（13）
此时，r m 便是原始信号S r+n (x)经过经验模态分解后剩下的残差信号，表示为：Res；
S2-2-2、信号重构
经过上述信号经验模态分解后，每一个本征模态函数IMF都表示原始信号中的一个独立频率成分，频率由高到低依次排列，因此，可以用分离出来的所有的IMF与残差Res对信号进行重构，得到去噪后的信号，重构公式表示如下：
Figure 844471DEST_PATH_IMAGE025
（14）
通过上述过程处理后，可以有效滤除高频随机噪声，得到仅保留反射信号信噪比的信号,即
Figure 990281DEST_PATH_IMAGE026
，表示为
Figure 897057DEST_PATH_IMAGE027
，并以此信号作为第一监测量。
5.根据权利要求4所述的基于北斗GEO卫星反射信号的洪涝监测方法，其特征在于，所述步骤S3中，反射信号信噪比幅度值的获取方法：
S3-1、利用信噪比线性-周期转换模型将以dB-Hz为单位的线性反射信号信噪比转换为以volts/volts为单位的周期反射信号信噪比，线性-周期转换模型的转换公式表示如下：
Figure 52095DEST_PATH_IMAGE028
（15）
式中，
Figure 727927DEST_PATH_IMAGE029
为仅包含反射信号信噪比的信号，该信号是以dB-Hz为单位的线性信号，
Figure 575797DEST_PATH_IMAGE030
为经过线性-周期转换模型后得到的以volts/volts为单位的周期反射信号信噪比；
S3-2、采用Lomb-Scargle谱分析法求卫星有效反射高
Figure 969870DEST_PATH_IMAGE031
；
S3-3、采用高斯粒子滤波算法求出反射信号信噪比幅度值，具体公式表示如下：
Figure 663019DEST_PATH_IMAGE032
式中，A和
Figure 196287DEST_PATH_IMAGE033
表示该信噪比的幅度和相位，是公式中的未知量，
Figure 949480DEST_PATH_IMAGE030
是周期反射信号信噪比，
Figure 830848DEST_PATH_IMAGE034
是卫星有效反射高，λ是信号的波长，属于已知量，EL是对应的高度角。
6.根据权利要求5所述的基于北斗GEO卫星反射信号的洪涝监测方法，其特征在于，所述步骤S3-2中卫星有效反射高
Figure 327688DEST_PATH_IMAGE035
的求解方法：
S3-2-1、通过Lomb-Scargle谱分析法求周期反射信号信噪比S r_v 的功率最大值时的频率
Figure 978113DEST_PATH_IMAGE036
，具体过程表示如下：
采用Lomb-Sacragle谱分析法求信号的功率谱，公式表示如下：
Figure 902206DEST_PATH_IMAGE037
(17)
式中，P x ( f )是频率为 f 的周期信号的功率；S r_v (x)为反射信号信噪比，N为序列长度，τ为时间平移不变量，可通过下式计算得到：
Figure 536450DEST_PATH_IMAGE038
通过式（17-18），可求出P x ( f )达到峰值时对应
Figure 571402DEST_PATH_IMAGE039
；
S3-2-2、利用
Figure 76333DEST_PATH_IMAGE040
和卫星有效反射高
Figure 171328DEST_PATH_IMAGE041
的关系求出卫星有效反射高，反射信号信噪比频率
Figure 761709DEST_PATH_IMAGE042
和卫星有效反射高的关系表示如下：
Figure 865931DEST_PATH_IMAGE043
式中，λ是信号的波长，属于已知量，通过上述过程可求出卫星的有效反射高。
7.根据权利要求1所述的基于北斗GEO卫星反射信号的洪涝监测方法，其特征在于，所述步骤S4中监测阈值为第一监测量的波动不超过2dB-Hz，第二监测量的波动不超过0.5。
8.根据权利要求6所述的基于北斗GEO卫星反射信号的洪涝监测方法，其特征在于，所述步骤S5包括如下步骤：
S5-1、采用基于北斗GEO卫星三频反射信号信噪比对洪涝进行监测，并根据监测结果进行标记，所述标记的级别包括一级标记和二级标记；
S5-2、采用基于北斗GEO卫星三频反射信号信噪比的幅度值对洪涝进行监测，并根据监测结果进行标记，所述标记的级别包括一级标记和二级标记；
S5-3、将经过第一监测量和第二监测量的洪涝监测结果进行联合处理，若第一监测量和第二监测量存在一级标记，则认为存在洪涝危险，需进行一级洪涝预警处理，若第一监测量和第二监测量无一级标记但有二级标记，则进行二级洪涝预警，若两者都没有标记，不进行预警。
9.根据权利要求8所述的基于北斗GEO卫星反射信号的洪涝监测方法，其特征在于，所述步骤S5-1中对洪涝进行监测的方法为：
首先，获取当前历元时刻的原始信噪比通过步骤S2计算得到当前历元时刻的第一监测量，并计算该历元时刻的高度角信息和时间信息；
其次，将计算得到的北斗GEO卫星三频，即B1频段、B2频段和B3频段的高度角信息和时间信息，并输入参考监测模型中搜索相应参考值中的第一监测量，并将计算得到的当前历元时刻的第一监测量与参考值中的第一监测量作对比，若当前历元时刻的第一监测量波动大于2dB-Hz时存在洪涝，则进行标记；
最后，根据历元标记结果进行统计分析，若标记历元占总历元数的80%以上，则认为存在强洪涝相关，做一级标记，若标记历元占总历元数的60%-80%以内，则认为存在较强洪涝相关，做二级标记，若标记历元占比低于60%，则不进行标记，最后，将监测结果进行保存并输出。
10.根据权利要求9所述的基于北斗GEO卫星反射信号的洪涝监测方法，其特征在于，所述步骤S5-2中对洪涝进行监测的方法为：
首先，根据步骤S5-1计算得到的当前历元的反射信号信噪比通过步骤S3计算得到反射信号信噪比的幅度值，即当前历元时刻的第二检测量，并计算该历元时刻的高度角和时间信息；
其次，将计算得到的北斗GEO卫星三频即B1频段、B2频段和B3频段的高度角信息和时间信息，并输入参考监测模型中搜索相应参考值中的第二监测量，并将计算得到的当前历元时刻的第二监测量与参考值中的第二监测量作对比，若当前历元时刻的反射信号信噪比幅度值波动大于0.5时存在洪涝，则进行标记；
最后，根据历元标记结果进行统计分析，若标记历元占总历元数的80%以上，则认为存在强洪涝相关，做一级标记，若标记历元占总历元数的60%-80%以内，则认为存在较强洪涝相关，做二级标记，若标记历元占比低于60%，则不进行标记，最后，将监测结果进行保存并输出。
CN202210887952.4A 2022-07-27 2022-07-27 基于北斗geo卫星反射信号的洪涝监测方法 Pending CN115201879A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210887952.4A CN115201879A (zh) 2022-07-27 2022-07-27 基于北斗geo卫星反射信号的洪涝监测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210887952.4A CN115201879A (zh) 2022-07-27 2022-07-27 基于北斗geo卫星反射信号的洪涝监测方法

Publications (1)

Publication Number Publication Date
CN115201879A true CN115201879A (zh) 2022-10-18

Family

ID=83584697

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210887952.4A Pending CN115201879A (zh) 2022-07-27 2022-07-27 基于北斗geo卫星反射信号的洪涝监测方法

Country Status (1)

Country Link
CN (1) CN115201879A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116400385A (zh) * 2023-03-21 2023-07-07 湖北珞珈实验室 一种底层大气与电离层耦合异常探测系统及方法

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116400385A (zh) * 2023-03-21 2023-07-07 湖北珞珈实验室 一种底层大气与电离层耦合异常探测系统及方法
CN116400385B (zh) * 2023-03-21 2024-01-12 湖北珞珈实验室 一种底层大气与电离层耦合异常探测系统及方法

Similar Documents

Publication Publication Date Title
Labroue et al. First quality assessment of the Cryosat-2 altimetric system over ocean
Wen et al. How well were the early 2017 California Atmospheric River precipitation events captured by satellite products and ground‐based radars?
Laabs et al. Timing of the last glaciation and subsequent deglaciation in the Ruby Mountains, Great Basin, USA
CN113075706A (zh) 一种基于gnss-r的雪深反演方法及其应用
Wang et al. Can the GPM IMERG hourly products replicate the variation in precipitation during the wet season over the Sichuan Basin, China?
CN115201879A (zh) 基于北斗geo卫星反射信号的洪涝监测方法
CN117516636B (zh) 一种海岸堤坝安全监测预警方法及系统
Nair et al. Exploring the potential of SWOT mission for reservoir monitoring in Mahanadi basin
Reagan et al. World Ocean Atlas 2023, Volume 2: Salinity
Mohamad et al. Monitoring groundwater depletion due to drought using satellite gravimetry: A review
Sun et al. Construction of the Mean Sea Surface model combined HY-2A with DTU18 MSS in the Antarctic ocean
Chen et al. Multisensor remote sensing and the multidimensional modeling of extreme flood events: a case study of hurricane harvey–triggered floods in Houston, Texas, USA
Li et al. GNSS-IR dual-frequency data fusion for soil moisture inversion based on Helmert variance component estimation
Wang et al. A study of the relationship between surface roughness and GNSS-R coherent returns over land
Gong et al. Methods of INSAR atmosphere correction for volcano activity monitoring
Durmaz Non-parametric and semi-parametric regional modeling of the ionospheric vertical total electron content using ground-based gps observations
Shakya Estimation of Stage-Area-Storage Relationships in Reservoirs and Stage-Discharge Relationships in Rivers Using Remotely Sensed Data
Ma et al. Using CYGNSS and L-band Radiometer Observations to Retrieve Surface Water Fraction: A Case Study of the Catastrophic Flood of 2022 in Pakistan
CN114254568B (zh) 基于人工智能决策树的gps遥感洪水预警方法
CN114137575B (zh) 顾及卫星偏差和载噪比弧段影响的洪水探测方法
Yang Satellite Altimetry Applications on Lake Ice Thickness and Land Subsidence
Khajeh et al. Interferometric path models for GNSS ground-based phase altimetry
Donor River Discharge Estimation Using Multisource Remote Sensing: Application in the Lancang Mekong River Basin
Zheng et al. Research on GNSS-IR soil moisture retrieval based on random forest algorithm
Karimi Sea level variability and mean sea level determination around Australia from satellite altimetry

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