CN101916439B - 基于希尔伯特-黄变换的空间碎片高光谱序列探测方法 - Google Patents

基于希尔伯特-黄变换的空间碎片高光谱序列探测方法 Download PDF

Info

Publication number
CN101916439B
CN101916439B CN2010102369741A CN201010236974A CN101916439B CN 101916439 B CN101916439 B CN 101916439B CN 2010102369741 A CN2010102369741 A CN 2010102369741A CN 201010236974 A CN201010236974 A CN 201010236974A CN 101916439 B CN101916439 B CN 101916439B
Authority
CN
China
Prior art keywords
curve
spectrum
imf component
amplitude
hyperspectral
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.)
Expired - Fee Related
Application number
CN2010102369741A
Other languages
English (en)
Other versions
CN101916439A (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.)
Harbin Institute of Technology
Original Assignee
Harbin Institute of Technology
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 Harbin Institute of Technology filed Critical Harbin Institute of Technology
Priority to CN2010102369741A priority Critical patent/CN101916439B/zh
Publication of CN101916439A publication Critical patent/CN101916439A/zh
Application granted granted Critical
Publication of CN101916439B publication Critical patent/CN101916439B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Image Analysis (AREA)

Abstract

基于希尔伯特-黄变换的空间碎片高光谱序列探测方法，属于图像处理领域，本发明为解决采用传统的危险空间碎片探测方法时，必须依赖有关空间碎片图像的样本信息，算法适应性差的问题，本发明方法包括：一、连续T次采样疑似目标中心位置的高光谱曲线，并将获取的T段高光谱曲线处理合成为待处理曲线；二、一维经验模态分解，获取两个本征模态函数分量；三、对二阶IMF分量进行希尔伯特变换，得到其幅值和瞬时频率；四、保留幅值高于均值一半的部分，保留瞬时频率高于均值的部分，并分割形成特征波段集合；五、循环搜索特征波段集合，判断目标是否旋转，如旋转，则确定疑似目标为空间碎片，如不旋转，则确定疑似目标不是空间碎片。

Description

基于希尔伯特-黄变换的空间碎片高光谱序列探测方法
技术领域
本发明涉及基于希尔伯特-黄变换的空间碎片高光谱序列探测方法，属于图像处理领域。
背景技术
航天技术的日益发展有利地促进了世界各国的经济建设和社会进步，但与此同时，由于人类对太空开发的狂热，致使大量的人造物体被发射到太空之中。随着航天器发射数量的增加，失效的航天器，以及脱落、碰撞、爆炸等形成的空间碎片（space debris）的数量也在逐年增加。大、中型空间碎片主要包括废弃的卫星、火箭及其脱落物等；小型空间碎片主要包括卫星表面脱落物、发动机喷出的碎屑、碰撞或爆炸产生的碎片等等。
探测是掌握空间碎片分布情况的前提和基础，也是空间碎片观测研究中至关重要的工作。质量很大的大型碎片主要在地面进行观测，也较为容易识别，测量方法主要是可见光观测、雷达观测等。数量众多的中型空间碎片（尺寸在1mm-10cm之间），由于尺度较小，地基望远镜和雷达一般无法观测或测定其轨道，而航天器表面采样分析方法对其撞击后再进行分析又可能为时已晚，因此他们常被称为危险空间碎片。目前探测危险空间碎片的研究方法集中在天基可见光遥感探测，无论是探测方法还是识别方法都处于发展的初级阶段，且可见光相机较高光谱成像仪在视域信息提供量上存在较大差距。在传统的危险空间碎片（中型空间碎片）探测手段中，主要依靠天基可见光相机的成像序列并利用空间碎片的运动特性进行探测，必须依赖有关空间碎片图像的样本信息，对疑似目标只能获取其位置及灰度信息，天基可见光相机无法利用空间碎片围绕它的主轴进行自旋运动，算法适应性差。
发明内容
本发明目的是为了解决采用传统的危险空间碎片探测方法时，天基可见光相机无法利用空间碎片围绕它的主轴进行自旋运动，必须依赖有关空间碎片图像的样本信息，算法适应性差的问题，提供了基于希尔伯特-黄变换的空间碎片高光谱序列探测方法。
本发明包括如下步骤：
步骤一、连续T次采样疑似目标中心位置的高光谱曲线，并将获取的T段高光谱曲线处理合成为待处理曲线；
步骤二、对步骤一获取的待处理曲线数据进行一维经验模态分解，获取两个本征模态函数分量和残差，两个本征模态函数分量为一阶IMF分量和二阶IMF分量；
步骤三、对所述二阶IMF分量进行希尔伯特变换，得到二阶IMF分量的幅值和瞬时频率；
步骤四、保留同时满足幅值高于均值一半和瞬时频率高于均值两个条件的二阶IMF分量，并分割形成特征波段集合；
步骤五、循环搜索特征波段集合，判断目标是否旋转，进而确定疑似目标是否为空间碎片，如旋转，则确定疑似目标为空间碎片，如不旋转，则确定疑似目标不是空间碎片。
本发明的优点：
1）本发明利用高光谱成像仪进行天基采样，能够对目标因旋转所产生的光谱微弱变化进行检测，克服了天基可见光相机无法利用空间碎片围绕它的主轴进行自旋运动这一事实来进行探测，使得空间碎片的判定能够基于对象的特殊运动信息，而无需任何有关空间碎片图像的样本信息，因而算法适应性强。
2）本发明所提出的探测方法，在确定空间碎片的同时，不但可以通过高光谱图像信息定位对象的位置，而且可以估计对象的自旋周期，可为空间碎片的天基探测搜集更多的目标信息。
附图说明
图1为本发明方法流程图；
图2为一维经验模态分解流程图；
图3为实施例中观测像素的20段原始采样数据合成后的高光谱曲线；
图4为合成高光谱曲线的一阶本征模态函数（IMF1）曲线；
图5为合成高光谱曲线的二阶本征模态函数（IMF2）曲线；
图6为合成高光谱曲线的残差函数（RES）曲线；
图7为IMF2经希尔伯特变换后的幅值曲线；
图8为IMF2经希尔伯特变换后的瞬时频率曲线；
图9为本发明方法筛选后的特征波段定位图。
具体实施方式
具体实施方式一：下面结合图1和图2说明本实施方式，
1998年NASA的黄锷（Norden E Huang）提出了一种依据数据本身的时间尺度特征来对信号进行分解的新算法——希尔伯特-黄变换（Hilbert-Huang Transform，HHT）。其中经验模态分解法（Empirical Mode Decomposition，EMD）是HHT算法的核心步骤，利用信号内部时间尺度的变化做能量与频率的解析，将信号展开成若干个本征模态函数（Intrinsic Mode Function，IMF），再利用希尔伯特变换（Hilbert Transform，HT）获得IMF的瞬时频率及振幅，最后求得时间-频率-振幅的三维谱分布，即Hilbert谱。
其中，IMF必须满足下列条件：
1）在整个函数中，极值点的数目与穿越零点的数目相等或者相差1；
2）在任何时刻，由局部极值包络线所定义的包络线局部均值为零。
依托这两个条件构建起来的EMD及HHT被认为是对以傅立叶变换为基础的线性及稳态谱分析的重大突破，是近年来强有力地求解非线性、非平稳信号的自适应方法。
在传统的危险空间碎片（中型空间碎片）探测手段中，主要依靠天基可见光相机的成像序列并利用空间碎片的运动特性进行探测，对疑似目标只能获取其位置及灰度信息。而遥感探测器若改用高光谱成像仪，则可以增加捕获对象的光谱信息，进而可利用空间碎片通常会围绕它的主轴进行简单的自旋运动这一事实来进行探测。由于目标像素的光谱序列属于非线性、非平稳信号，因此采用HHT进行自旋特征提取进而检测出空间碎片是非常适合的。
本实施方式方法通过以下技术方案实现的：选定观察对象中心位置处目标像素的                                                
Figure 553860DEST_PATH_IMAGE001
段连续高光谱数据序列，合成1段数据后进行一维经验模态分解（Empirical Mode Decomposition，EMD），对得到的二阶本征模态函数（Intrinsic Mode Function，IMF）进行希尔伯特变换（Hilbert Transform，HT）获得其瞬时频率及振幅，进而综合判断特征位置并分割为
Figure 471000DEST_PATH_IMAGE001
段，各段反映的是每个采样时刻得到的高光谱数据的高频特征信息所在波段，通过循环搜索判断目标是否旋转从而检测出空间碎片。具体方法包括以下步骤：
步骤一、连续T次采样疑似目标中心位置的高光谱曲线，并将获取的T段高光谱曲线处理合成为待处理曲线；
步骤二、对步骤一获取的待处理曲线数据进行一维经验模态分解，获取两个本征模态函数分量和残差，两个本征模态函数分量为一阶IMF分量和二阶IMF分量；
步骤三、对所述二阶IMF分量进行希尔伯特变换，得到二阶IMF分量的幅值和瞬时频率；
步骤四、保留同时满足幅值高于均值一半和瞬时频率高于均值两个条件的二阶IMF分量，并分割形成特征波段集合；
步骤五、循环搜索特征波段集合，判断目标是否旋转，进而确定疑似目标是否为空间碎片，如旋转，则确定疑似目标为空间碎片，如不旋转，则确定疑似目标不是空间碎片。
高光谱成像仪每次采样会得到观测区域的若干幅连续光谱图像，选择其中疑似目标的中心位置像素，可以得到一条光谱曲线。连续
Figure 562322DEST_PATH_IMAGE001
次观测采样，则可得到该位置像素的
Figure 19848DEST_PATH_IMAGE001
段光谱曲线。如果疑似目标为空间碎片，大多数情况下都会发生自旋，表现在同一位置像素的光谱曲线上就是会有周期性的变化。下面在步骤一中详细说明。
步骤一中的连续T次采样疑似目标中心位置的高光谱曲线，并将获取的T段高光谱曲线处理合成为待处理曲线的过程为：
步骤11、采用高光谱成像仪采集观测区域的多幅连续高光谱图像，选择其中疑似目标的中心位置像素，获取一条疑似目标中心位置的高光谱曲线，所述疑似目标中心位置的高光谱曲线的长度为m，
步骤12、按照步骤11的方法连续T次采样疑似目标中心位置的高光谱曲线，获取具有连续数据的T段高光谱曲线，
步骤13、对每段高光谱曲线进行左右翻转，并连接在原始高光谱曲线之后，形成T段长度为2m的新的高光谱曲线，
步骤14、将所述T段长度为2m的新的高光谱曲线依次拼接，合成为待处理曲线。待处理曲线的长度为
Figure 834220DEST_PATH_IMAGE002
。
步骤二中的获取两个本征模态函数分量和残差的过程为：
设定输入的待处理曲线为
Figure 469732DEST_PATH_IMAGE003
，
Figure 799082DEST_PATH_IMAGE004
，
步骤21、IMF分解过程初始化：
Figure 309567DEST_PATH_IMAGE005
，且满足关系式成立，其中
Figure 34126DEST_PATH_IMAGE007
为第
Figure 601505DEST_PATH_IMAGE008
次分解后剩余的残差函数；
步骤22、筛选过程初始化，
Figure 338517DEST_PATH_IMAGE009
，且满足关系式
Figure 924219DEST_PATH_IMAGE010
成立，其中
Figure 354063DEST_PATH_IMAGE011
为第次本征模态函数分解中经过第
Figure 244714DEST_PATH_IMAGE013
次筛选后的剩余函数；
步骤23、根据筛选程序获取输入的待处理曲线
Figure 888185DEST_PATH_IMAGE003
经过第次本征模态函数分解的剩余的残差函数中经过第
Figure 827639DEST_PATH_IMAGE015
次筛选后的剩余函数
Figure 906453DEST_PATH_IMAGE016
；
步骤24、采用标准偏差准则判断输入的待处理曲线
Figure 716015DEST_PATH_IMAGE003
经过第
Figure 549979DEST_PATH_IMAGE014
次本征模态函数分解的剩余的残差函数中经过第
Figure 500618DEST_PATH_IMAGE015
次筛选后的剩余函数
Figure 196172DEST_PATH_IMAGE017
是否满足本征模态函数的条件，即
Figure 610973DEST_PATH_IMAGE018
是否小于阈值
Figure 553521DEST_PATH_IMAGE019
，
Figure 303040DEST_PATH_IMAGE020
；
判断结果为是，执行步骤25，判断结果为否，则
Figure 785974DEST_PATH_IMAGE021
，然后执行步骤23，
步骤25、提取一个本征模态函数分量
Figure 258544DEST_PATH_IMAGE022
；
步骤26、获取输入的待处理曲线
Figure 185043DEST_PATH_IMAGE003
经过第
Figure 172590DEST_PATH_IMAGE014
次本征模态函数分解的剩余的残差函数
Figure 396898DEST_PATH_IMAGE023
；
步骤27、判断是否满足下述关系式成立：
Figure 769980DEST_PATH_IMAGE024
；
判断结果为否，则
Figure 116647DEST_PATH_IMAGE025
，然后执行步骤22，判断结果为是，完成提取过程，获得两个本征模态函数分量：一阶IMF分量、二阶IMF分量
Figure 529174DEST_PATH_IMAGE026
；和
Figure 370222DEST_PATH_IMAGE027
个经过第2次本征模态函数分解的剩余的残差RES：
Figure 614122DEST_PATH_IMAGE028
。即
Figure 803795DEST_PATH_IMAGE029
。
步骤三所述的获取二阶IMF分量的幅值和瞬时频率的过程为：
步骤31、对二阶IMF分量
Figure 15202DEST_PATH_IMAGE030
进行离散褶积，获得其希尔伯特变换
Figure 643630DEST_PATH_IMAGE031
；
Figure 945298DEST_PATH_IMAGE032
，
步骤32、获取二阶IMF分量的幅值：
Figure 118921DEST_PATH_IMAGE030
的解析信号为
Figure 568357DEST_PATH_IMAGE033
，将所述解析信号的包络振幅
Figure 938159DEST_PATH_IMAGE034
作为二阶IMF分量的幅值，解析信号的包络振幅
Figure 140339DEST_PATH_IMAGE034
按如下公式计算：
Figure 999710DEST_PATH_IMAGE035
，
步骤33、获取二阶IMF分量的瞬时频率
Figure 608546DEST_PATH_IMAGE036
：先求取所述解析信号的相角
Figure 595088DEST_PATH_IMAGE037
：
，
然后，根据下式获取二阶IMF分量的瞬时频率
Figure 636042DEST_PATH_IMAGE036
：
Figure 43759DEST_PATH_IMAGE039
。
合成高光谱曲线记录了同一目标位置像素的连续采样光谱序列，如果目标发生自旋，则该序列会存在循环的光谱。但是同样环境下采集的光谱曲线是具有高度相关性的，只在一部分高频信息所在波段存在差异，因此自适应地提取高频信息所在波段是非常有意义的，可以用这些波段的集合作为各采样时刻的特征位置，下面在步骤四中详细说明。
步骤四所述的特征波段集合的获取过程为：
步骤41、计算二阶IMF分量的幅值平均值和瞬时频率的平均值
Figure 682867DEST_PATH_IMAGE041
：
Figure 634774DEST_PATH_IMAGE042
，其中
Figure 280519DEST_PATH_IMAGE043
，
Figure 530235DEST_PATH_IMAGE044
，
步骤42、舍弃二阶IMF分量中幅值小于幅值平均值一半的波段，同时舍弃瞬时频率小于瞬时频率平均值，获取特征时刻集合
Figure 269893DEST_PATH_IMAGE045
：
Figure 136348DEST_PATH_IMAGE046
，
步骤43、在范围内依次截取长度为的集合，获得
Figure 976129DEST_PATH_IMAGE001
段采样数据各自对应的特征波段集合。
如果对象为空间碎片，在自旋一周之后，相同位置的观察像素会具有一定相关性，即便观察像素在碎片上的位置稍有偏移也是如此，只是相关性有所减弱。下面在步骤五进行详细说明。
步骤五所述的循环搜索特征波段集合，判断目标是否旋转，进而确定疑似目标是否为空间碎片的过程为：
步骤51、计算特征波段集合
Figure 123131DEST_PATH_IMAGE048
两两之间特征集合的相异波段数目，
步骤52、计算所述相异波段数目的均值
Figure 736832DEST_PATH_IMAGE050
；
步骤53、计算等间隔采样数据之间的相异波段数目采用下式计算均值
Figure 201442DEST_PATH_IMAGE051
：
Figure 246759DEST_PATH_IMAGE052
，
由于采样数据量为，循环周期的有效数值经验证应在
Figure 588100DEST_PATH_IMAGE053
到之间。
步骤54、获取
Figure 5492DEST_PATH_IMAGE055
中最小值所对应的间隔
Figure 17441DEST_PATH_IMAGE056
：
Figure 420741DEST_PATH_IMAGE057
，
在到
Figure 807915DEST_PATH_IMAGE054
之间选取多个k值，计算出每个k值对应的均值
Figure 544926DEST_PATH_IMAGE051
，并将其中最小均值找出来，这个最小均值对应的间隔是用于衡量相关性的一个周期参数。
步骤55、判断
Figure 560473DEST_PATH_IMAGE058
是否成立，
判断结果为是，断定目标旋转导致观测像素在间隔
Figure 349568DEST_PATH_IMAGE059
的周期具有较高的光谱相关性，使得特征集合相异波段数目减少，进而判断疑似目标为空间碎片；判断结果为否，判断疑似目标不是空间碎片。
具体实施方式二、下面结图1至图9，结合某空间碎片模型的高光谱旋转成像序列给出一个具体实施例：
执行步骤一：选定观察对象中心位置处目标像素的多个连续高光谱数据序列进行数据合成，形成待处理曲线。
对连续采集的20个高光谱数据，调用目标中心位置像素的高光谱曲线（共计20段，每段长度m为220），接着对每段高光谱曲线进行左右翻转，并连接在原始光谱曲线之后，形成新的高光谱曲线（共计20段，每段长度为2m=440），最后将20段新的光谱曲线依次连接得到一条合成曲线为待处理曲线（如图3所示，长度=8800）作为后续步骤的一维输入光谱数据——待处理曲线。
执行步骤二：对合成的一维输入光谱数据进行一维经验模态分解，分解至一阶、二阶本征模态函数以及残差函数即停止分解，具体描述如下：
设定输入的待处理曲线为，
Figure 242810DEST_PATH_IMAGE060
，
步骤21、IMF分解过程初始化：
Figure 34048DEST_PATH_IMAGE005
，且满足关系式
Figure 112863DEST_PATH_IMAGE006
成立，其中
Figure 423890DEST_PATH_IMAGE007
为第次分解后剩余的残差函数；
步骤22、筛选过程初始化，
Figure 208492DEST_PATH_IMAGE009
，且满足关系式成立，其中为第
Figure 759931DEST_PATH_IMAGE012
次本征模态函数分解中经过第
Figure 10915DEST_PATH_IMAGE013
次筛选后的剩余函数；
步骤23、根据筛选程序获取输入的待处理曲线
Figure 493849DEST_PATH_IMAGE003
经过第
Figure 966418DEST_PATH_IMAGE014
次本征模态函数分解的剩余的残差函数中经过第次筛选后的剩余函数
Figure 113421DEST_PATH_IMAGE016
；
步骤24、采用标准偏差准则判断输入的待处理曲线经过第
Figure 477854DEST_PATH_IMAGE014
次本征模态函数分解的剩余的残差函数中经过第
Figure 824522DEST_PATH_IMAGE015
次筛选后的剩余函数是否满足本征模态函数的条件，即是否小于阈值，
Figure 10204DEST_PATH_IMAGE061
；
判断结果为是，执行步骤25，判断结果为否，则
Figure 723076DEST_PATH_IMAGE021
，然后执行步骤23，
步骤25、提取一个本征模态函数分量
Figure 351504DEST_PATH_IMAGE022
；
步骤26、获取输入的待处理曲线
Figure 387593DEST_PATH_IMAGE003
经过第
Figure 325331DEST_PATH_IMAGE014
次本征模态函数分解的剩余的残差函数；
步骤27、判断是否满足下述关系式成立：
Figure 878989DEST_PATH_IMAGE024
；
判断结果为否，则
Figure 848213DEST_PATH_IMAGE025
，然后执行步骤22，判断结果为是，完成提取过程，获得两个本征模态函数分量：一阶IMF分量（见图4）、二阶IMF分量（见图5）：
Figure 707585DEST_PATH_IMAGE026
；和
Figure 316421DEST_PATH_IMAGE027
个经过第2次本征模态函数分解的剩余的残差RES（见图6）：
Figure 801497DEST_PATH_IMAGE028
。即
Figure 874496DEST_PATH_IMAGE029
。
执行步骤三：对IMF2进行希尔伯特变换，得到二阶IMF分量的幅值和瞬时频率。
首先，按如下算式对二阶IMF分量：
Figure 655501DEST_PATH_IMAGE062
进行离散褶积得到其希尔伯特变换
Figure 751633DEST_PATH_IMAGE031
：
Figure 259975DEST_PATH_IMAGE063
；
接着，按如下算式计算的解析信号
Figure 841184DEST_PATH_IMAGE065
的包络振幅
Figure 486929DEST_PATH_IMAGE066
：
Figure 549694DEST_PATH_IMAGE067
；
然后，按如下算式计算解析信号的相角
Figure 597284DEST_PATH_IMAGE068
和瞬时频率
Figure 907043DEST_PATH_IMAGE069
：
Figure 289351DEST_PATH_IMAGE070
，
Figure 405075DEST_PATH_IMAGE071
；
该步骤完毕，得到二阶IMF分量的幅值
Figure 244855DEST_PATH_IMAGE066
和瞬时频率
Figure 538564DEST_PATH_IMAGE072
，分别如图7和图8所示。
执行步骤四：综合判断高光谱数据的高频信息所在波段，得到特征波段集合。
对于二阶IMF分量的瞬时频率中小于其平均值（
Figure 893322DEST_PATH_IMAGE073
）的波段予以排除以保留具有特征意义的高频信息所在波段；同时由于二阶IMF分量的幅值和瞬时频率是依波段对应的，而过小的幅值其相应波段以噪声成分居多，因此对于幅值小于其平均值（
Figure 750420DEST_PATH_IMAGE074
）二分之一的波段也要予以排除，最后得到特征时刻集合（见图9）：
，
然后对
Figure 515485DEST_PATH_IMAGE076
依次截取长度为440的集合，即得到20段采样数据各自对应的特征波段集合
Figure 723744DEST_PATH_IMAGE077
。
执行步骤五：循环搜索判断目标是否旋转，进而确定是否为空间碎片。
首先，计算两两之间特征集合的相异波段数目
Figure 422895DEST_PATH_IMAGE078
；
接着，计算相异波段数目的均值
Figure 1513DEST_PATH_IMAGE079
；
由于采样数据量为
Figure 997151DEST_PATH_IMAGE080
，循环周期的有效数值经验证应在到
Figure 738022DEST_PATH_IMAGE082
之间，对等间隔采样数据之间的相异波段数目采用下式计算均值：
Figure 289089DEST_PATH_IMAGE083
；
找出
Figure 291680DEST_PATH_IMAGE084
（见表1）中最小值所对应的间隔
Figure 126650DEST_PATH_IMAGE085
；
由于
Figure 290915DEST_PATH_IMAGE086
大于阈值2，因此可判断目标旋转导致观测像素在间隔10的周期具有较高的光谱相关性，使得特征集合相异波段数目减少，进而判断目标为空间碎片。
表1  等间隔采样数据之间的相异波段数目均值
间隔k 相异波段数目均值Mk
2 46.23
3 47.86
4 50.05
5 36.17
6 48.17
7 51.95
8 48.56
9 48.69
10 16.30
11 52.33
12 55.13
13 49.43
14 56.83
15 44.20
16 39.50
17 57.33
18 52.00

Claims (6)

1.基于希尔伯特-黄变换的空间碎片高光谱序列探测方法，其特征在于，它包括如下步骤：
步骤一、连续T次采样疑似目标中心位置的高光谱曲线，并将获取的T段高光谱曲线处理合成为待处理曲线；
步骤二、对步骤一获取的待处理曲线数据进行一维经验模态分解，获取两个本征模态函数分量和残差，两个本征模态函数分量为一阶IMF分量和二阶IMF分量；
步骤三、对所述二阶IMF分量进行希尔伯特变换，得到二阶IMF分量的幅值和瞬时频率；
步骤四、保留同时满足幅值高于均值一半和瞬时频率高于均值两个条件的二阶IMF分量，并分割形成特征波段集合；
步骤五、循环搜索特征波段集合，判断目标是否旋转，进而确定疑似目标是否为空间碎片，如旋转，则确定疑似目标为空间碎片，如不旋转，则确定疑似目标不是空间碎片，
其过程为：
步骤51、计算特征波段集合Ψ1，Ψ2，...，ΨT两两之间特征集合的相异波段数目Dij，1≤i＜j≤T，
步骤52、计算所述相异波段数目的均值 M &OverBar; = mean { D ij | 1 &le; i < j &le; T } ;
步骤53、计算等间隔采样数据之间的相异波段数目采用下式计算均值Mk：
Mk＝mean{Dij|i，j间隔为k，1≤i＜j≤T}，0.1T≤k≤0.9T，
步骤54、获取Mk中最小值所对应的间隔k*：
k*＝arg min{Mk|0.1T≤k≤0.9T}，
步骤55、判断是否成立，
判断结果为是，断定目标旋转导致观测像素在间隔k*的周期具有较高的光谱相关性，使得特征集合相异波段数目减少，进而判断疑似目标为空间碎片；判断结果为否，判断疑似目标不是空间碎片。
2.根据权利要求1所述的基于希尔伯特-黄变换的空间碎片高光谱序列探测方法，其特征在于，步骤一中的连续T次采样疑似目标中心位置的高光谱曲线，并将获取的T段高光谱曲线处理合成为待处理曲线的过程为：
步骤11、采用高光谱成像仪采集观测区域的多幅连续高光谱图像，选择其中疑似目标的中心位置像素，获取一条疑似目标中心位置的高光谱曲线，所述疑似目标中心位置的高光谱曲线的长度为m，
步骤12、按照步骤11的方法连续T次采样疑似目标中心位置的高光谱曲线，获取具有连续数据的T段高光谱曲线，
步骤13、对每段高光谱曲线进行左右翻转，并连接在原始高光谱曲线之后，形成T段长度为2m的新的高光谱曲线，
步骤14、将所述T段长度为2m的新的高光谱曲线依次拼接，合成为待处理曲线。
3.根据权利要求1所述的基于希尔伯特-黄变换的空间碎片高光谱序列探测方法，其特征在于，步骤二中的获取两个本征模态函数分量和残差的过程为：
设定输入的待处理曲线为x(t)，t＝1，2，...，N，
步骤21、经验模态分解过程初始化：n＝1，且满足关系式rn-1(t)＝x(t)成立，其中rn-1(t)为第(n-1)次经验模态分解后剩余的残差函数；
步骤22、筛选过程初始化，k＝1，且满足关系式hn(k-1)(t)＝rn-1(t)成立，其中hn(k-1)(t)为第n次本征模态函数分解中经过第(k-1)次筛选后的剩余函数；
步骤23、根据筛选程序获取输入的待处理曲线x(t)经过第n次本征模态函数分解中经过第k次筛选后的剩余函数hnk(t)；
步骤24、采用标准偏差准则判断输入的待处理曲线x(t)经过第n次本征模态函数分解的剩余的残差函数中经过第k次筛选后的剩余函数hnk(t)是否满足本征模态函数的条件，即 ( h n ( k - 1 ) ( t ) - h nk ( t ) ) 2 / h n ( k - 1 ) 2 ( t ) 是否小于阈值HSD，0.2≤HSD≤0.3；
判断结果为是，执行步骤25，判断结果为否，则k＝k+1，然后执行步骤23，
步骤25、提取一个本征模态函数分量cn(t)＝hnk(t)；
步骤26、获取输入的待处理曲线x(t)经过第n次本征模态函数分解的剩余的残差函数rn(t)＝rn-1(t)-cn(t)；
步骤27、判断是否满足下述关系式成立：n≥2；
判断结果为否，则n＝n+1，然后执行步骤22，判断结果为是，完成提取过程，获得两个本征模态函数分量：一阶IMF分量、二阶IMF分量{c1(t)，c2(t)}；和1个经过第2次本征模态函数分解的剩余的残差RES：r2(t)。
4.根据权利要求3所述的基于希尔伯特-黄变换的空间碎片高光谱序列探测方法，其特征在于，步骤24的中HSD＝0.25。
5.根据权利要求1所述的基于希尔伯特-黄变换的空间碎片高光谱序列探测方法，其特征在于，步骤三所述的获取二阶IMF分量的幅值和瞬时频率的过程为：
步骤31、对二阶IMF分量c2(t)进行离散褶积，获得其希尔伯特变换y(t)；
y ( t ) = c 2 ( t ) &times; 1 &pi;t ,
步骤32、获取二阶IMF分量的幅值：c2(t)的解析信号为c2(t)+jy(t)，将所述解析信号的包络振幅a(t)作为二阶IMF分量的幅值，解析信号的包络振幅a(t)按如下公式计算：
a ( t ) = c 2 2 ( t ) + y 2 ( t ) ,
步骤33、获取二阶IMF分量的瞬时频率f(t)：先求取所述解析信号的相角θ(t)：
&theta; ( t ) = arctan ( y ( t ) c 2 ( t ) ) ,
然后，根据下式获取二阶IMF分量的瞬时频率f(t)：
f ( t ) = 1 2 &pi; d dt &theta; ( t ) .
6.根据权利要求1所述的基于希尔伯特-黄变换的空间碎片高光谱序列探测方法，其特征在于，步骤四所述的特征波段集合的获取过程为：
步骤41、计算二阶IMF分量的幅值平均值
Figure FDA0000115533990000035
和瞬时频率的平均值
Figure FDA0000115533990000036
a &OverBar; = mean { a ( t ) | 1 &le; t &le; N } , 其中N＝2m×T，
f &OverBar; = mean { f ( t ) | 1 &le; t &le; N } ,
步骤42、舍弃二阶IMF分量中幅值小于幅值平均值
Figure FDA0000115533990000041
一半的波段，同时舍弃瞬时频率小于瞬时频率平均值
Figure FDA0000115533990000042
获取特征时刻集合Ψ：
&Psi; = { t f ( t ) &GreaterEqual; f &OverBar; a ( t ) &GreaterEqual; a &OverBar; 2 , 1 &le; t &le; N } ,
步骤43、在1≤t≤N范围内依次截取长度为2m的集合，获得T段采样数据各自对应的特征波段集合Ψ1，Ψ2，...，ΨT。
CN2010102369741A 2010-07-27 2010-07-27 基于希尔伯特-黄变换的空间碎片高光谱序列探测方法 Expired - Fee Related CN101916439B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN2010102369741A CN101916439B (zh) 2010-07-27 2010-07-27 基于希尔伯特-黄变换的空间碎片高光谱序列探测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN2010102369741A CN101916439B (zh) 2010-07-27 2010-07-27 基于希尔伯特-黄变换的空间碎片高光谱序列探测方法

Publications (2)

Publication Number Publication Date
CN101916439A CN101916439A (zh) 2010-12-15
CN101916439B true CN101916439B (zh) 2012-08-08

Family

ID=43323943

Family Applications (1)

Application Number Title Priority Date Filing Date
CN2010102369741A Expired - Fee Related CN101916439B (zh) 2010-07-27 2010-07-27 基于希尔伯特-黄变换的空间碎片高光谱序列探测方法

Country Status (1)

Country Link
CN (1) CN101916439B (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106373149B (zh) * 2015-07-23 2019-03-22 中国人民解放军海军大连舰艇学院 基于辅助图像分量和Hilbert变换的图像分解方法
CN105571412B (zh) * 2015-12-11 2017-04-12 中国人民解放军63850部队 基于Hilbert变换的弹丸进动周期提取方法
FR3071929A1 (fr) 2017-10-02 2019-04-05 Ecole Polytechnique Radar laser pour la detection des debris spatiaux
CN111198385A (zh) * 2019-12-26 2020-05-26 北京旷视机器人技术有限公司 障碍物检测方法、装置、计算机设备和存储介质
CN114298106B (zh) * 2021-12-29 2022-07-22 中国铁道科学研究院集团有限公司铁道建筑研究所 路基碾压中特征波识别方法
CN116992293B (zh) * 2023-09-26 2023-12-08 北京豪迈生物工程股份有限公司 一种用于化学发光仪器的智能化数据处理方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101685435A (zh) * 2008-09-26 2010-03-31 财团法人工业技术研究院 用于图像纹理分析的多维度经验模态分析方法
CN101763336A (zh) * 2010-01-18 2010-06-30 哈尔滨工业大学 一种具有激励噪声添加参数选取功能的集合经验模态分解方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
FR2831302A1 (fr) * 2001-10-19 2003-04-25 St Microelectronics Sa Codage d'informations concentriques

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101685435A (zh) * 2008-09-26 2010-03-31 财团法人工业技术研究院 用于图像纹理分析的多维度经验模态分析方法
CN101763336A (zh) * 2010-01-18 2010-06-30 哈尔滨工业大学 一种具有激励噪声添加参数选取功能的集合经验模态分解方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
JP特开2003-228707A 2003.08.15
王坚 等.基于经验模态分解的高分辨率影像融合.《遥感学报》.2007,第11卷(第1期), *
陈志刚 等.高光谱图像光谱域噪声去除的经验模态分解方法.《红外与毫米波学报》.2008,第27卷(第5期), *

Also Published As

Publication number Publication date
CN101916439A (zh) 2010-12-15

Similar Documents

Publication Publication Date Title
Dong et al. Early-season mapping of winter wheat in China based on Landsat and Sentinel images
CN101916439B (zh) 基于希尔伯特-黄变换的空间碎片高光谱序列探测方法
Wang et al. A combined GLAS and MODIS estimation of the global distribution of mean forest canopy height
CN102959354B (zh) 用于利用LiDAR数据来分析树冠层的方法和装置
CN109031344B (zh) 一种全波形激光雷达和高光谱数据联合反演森林结构参数的方法
Einzmann et al. Windthrow detection in European forests with very high-resolution optical data
CN108492260B (zh) 基于张量投票耦合霍夫变换的地质线性体提取方法
CN103400364B (zh) 一种森林资源变化监测方法
CN108766203B (zh) 一种紧致极化水稻制图的方法及系统
Zhang et al. Ground-penetrating radar railroad ballast inspection with an unsupervised algorithm to boost the region of interest detection efficiency
CN106778640B (zh) 一种三维可视化环境下地表植被覆盖模型的生成方法
CN103712560A (zh) 基于多个传感器信息融合的零件检测方法、系统及装置
CN103926203B (zh) 一种针对地物光谱不确定性的光谱角度制图方法
CN107764797A (zh) 一种基于低秩张量算法的拉曼光谱图像数据预处理方法
KR20170063182A (ko) 표적신호 분리를 이용하여 탐지성능이 향상된 초분광 영상의 표적물질 탐지방법
Lantini et al. A deep learning approach for tree root detection using GPR spectrogram imagery
Zhou et al. ICESat waveform-based land-cover classification using a curve matching approach
Ovakoglou et al. Spatial enhancement of Modis leaf area index using regression analysis with Landsat vegetation index
Xu et al. On-the-fly extraction of polyhedral buildings from airborne LiDAR data
CN102156872B (zh) 一种基于多光谱数据的物体识别方法和装置
Fang et al. Retrieval and mapping of heavy metal concentration in soil using time series landsat 8 imagery
Liu et al. The impacts of smoothing methods for time-series remote sensing data on crop phenology extraction
Yao et al. Corn area extraction by the integration of MODIS-EVI time series data and China’s environment satellite (HJ-1) data
CN105510900A (zh) 基于小波包分析神经网络的全波形激光雷达数据去噪方法
Olfindo Jr et al. Sugarcane plantation mapping using dynamic time warping from multi-temporal Sentinel-1A radar images

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20120808

CF01 Termination of patent right due to non-payment of annual fee