专利名称：基于希尔伯特-黄变换的空间碎片高光谱序列探测方法
技术领域：
本发明涉及基于希尔伯特-黄变换的空间碎片高光谱序列探测方法，属于图像处 理领域。
背景技术：
航天技术的日益发展有利地促进了世界各国的经济建设和社会进步，但与此同 时，由于人类对太空开发的狂热，致使大量的人造物体被发射到太空之中。随着航天器发射 数量的增加，失效的航天器，以及脱落、碰撞、爆炸等形成的空间碎片(space debris)的数 量也在逐年增加。大、中型空间碎片主要包括废弃的卫星、火箭及其脱落物等；小型空间碎 片主要包括卫星表面脱落物、发动机喷出的碎屑、碰撞或爆炸产生的碎片等等。探测是掌握空间碎片分布情况的前提和基础，也是空间碎片观测研究中至关重要 的工作。质量很大的大型碎片主要在地面进行观测，也较为容易识别，测量方法主要是可见 光观测、雷达观测等。数量众多的中型空间碎片(尺寸在Imm-IOcm之间)，由于尺度较小，地 基望远镜和雷达一般无法观测或测定其轨道，而航天器表面采样分析方法对其撞击后再进 行分析又可能为时已晚，因此他们常被称为危险空间碎片。目前探测危险空间碎片的研究 方法集中在天基可见光遥感探测，无论是探测方法还是识别方法都处于发展的初级阶段， 且可见光相机较高光谱成像仪在视域信息提供量上存在较大差距。在传统的危险空间碎片 (中型空间碎片)探测手段中，主要依靠天基可见光相机的成像序列并利用空间碎片的运动 特性进行探测，必须依赖有关空间碎片图像的样本信息，对疑似目标只能获取其位置及灰 度信息，天基可见光相机无法利用空间碎片围绕它的主轴进行自旋运动，算法适应性差。
发明内容
本发明目的是为了解决采用传统的危险空间碎片探测方法时，天基可见光相机无 法利用空间碎片围绕它的主轴进行自旋运动，必须依赖有关空间碎片图像的样本信息，算 法适应性差的问题，提供了基于希尔伯特_黄变换的空间碎片高光谱序列探测方法。本发明包括如下步骤
步骤一、连续T次采样疑似目标中心位置的高光谱曲线，并将获取的T段高光谱曲线处 理合成为待处理曲线；
步骤二、对步骤一获取的待处理曲线数据进行一维经验模态分解，获取两个本征模态 函数分量和残差，两个本征模态函数分量为一阶IMF分量和二阶IMF分量；
步骤三、对所述二阶IMF分量进行希尔伯特变换，得到二阶IMF分量的幅值和瞬时频
率；
步骤四、保留同时满足幅值高于均值一半和瞬时频率高于均值两个条件的二阶IMF分 量，并分割形成特征波段集合；
步骤五、循环搜索特征波段集合，判断目标是否旋转，进而确定疑似目标是否为空间碎 片，如旋转，则确定疑似目标为空间碎片，如不旋转，则确定疑似目标不是空间碎片。
6
本发明的优点
1)本发明利用高光谱成像仪进行天基采样，能够对目标因旋转所产生的光谱微弱变化 进行检测，克服了天基可见光相机无法利用空间碎片围绕它的主轴进行自旋运动这一事实 来进行探测，使得空间碎片的判定能够基于对象的特殊运动信息，而无需任何有关空间碎 片图像的样本信息，因而算法适应性强。2)本发明所提出的探测方法，在确定空间碎片的同时，不但可以通过高光谱图像 信息定位对象的位置，而且可以估计对象的自旋周期，可为空间碎片的天基探测搜集更多 的目标信息。
图1为本发明方法流程图2为一维经验模态分解流程图3为实施例中观测像素的20段原始采样数据合成后的高光谱曲线； 图4为合成高光谱曲线的一阶本征模态函数(IMFl)曲线； 图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进行自旋特征提取进而检测出空间碎片是非常 适合的。
7
本实施方式方法通过以下技术方案实现的选定观察对象中心位置处目标像素的 T段连续高光谱数据序列，合成1段数据后进行一维经验模态分解(Empirical Mode Decomposition, EMD)，对得到的二阶本征模态函数(Intrinsic Mode Function, IMF)进行 希尔伯特变换(Hilbert Transform，HT)获得其瞬时频率及振幅，进而综合判断特征位置并 分割为Γ段，各段反映的是每个采样时刻得到的高光谱数据的高频特征信息所在波段，通 过循环搜索判断目标是否旋转从而检测出空间碎片。具体方法包括以下步骤
步骤一、连续T次采样疑似目标中心位置的高光谱曲线，并将获取的T段高光谱曲线处 理合成为待处理曲线；
步骤二、对步骤一获取的待处理曲线数据进行一维经验模态分解，获取两个本征模态 函数分量和残差，两个本征模态函数分量为一阶IMF分量和二阶IMF分量；
步骤三、对所述二阶IMF分量进行希尔伯特变换，得到二阶IMF分量的幅值和瞬时频率； 步骤四、保留同时满足幅值高于均值一半和瞬时频率高于均值两个条件的二阶IMF分 量，并分割形成特征波段集合；
步骤五、循环搜索特征波段集合，判断目标是否旋转，进而确定疑似目标是否为空间碎 片，如旋转，则确定疑似目标为空间碎片，如不旋转，则确定疑似目标不是空间碎片。高光谱成像仪每次采样会得到观测区域的若干幅连续光谱图像，选择其中疑似目 标的中心位置像素，可以得到一条光谱曲线。连续Γ次观测采样，则可得到该位置像素的Γ 段光谱曲线。如果疑似目标为空间碎片，大多数情况下都会发生自旋，表现在同一位置像素 的光谱曲线上就是会有周期性的变化。下面在步骤一中详细说明。步骤一中的连续T次采样疑似目标中心位置的高光谱曲线，并将获取的T段高光 谱曲线处理合成为待处理曲线的过程为
步骤11、采用高光谱成像仪采集观测区域的多幅连续高光谱图像，选择其中疑似目标 的中心位置像素，获取一条疑似目标中心位置的高光谱曲线，所述疑似目标中心位置的高 光谱曲线的长度为m，
步骤12、按照步骤11的方法连续T次采样疑似目标中心位置的高光谱曲线，获取具有 连续数据的T段高光谱曲线，
步骤13、对每段高光谱曲线进行左右翻转，并连接在原始高光谱曲线之后，形成T段长 度为2m的新的高光谱曲线，
步骤14、将所述T段长度为2m的新的高光谱曲线依次拼接，合成为待处理曲线。待处 理曲线的长度为况=2mxr。步骤二中的获取两个本征模态函数分量和残差的过程为 设定输入的待处理曲线为χ Ο, = ^^--..,^ ,
步骤21、IMF分解过程初始化,且满足关系式^1(O = IG)成立，其中 ^1(Z)为第(11—1》次分解后剩余的残差函数
步骤22、筛选过程初始化，k=r且满足关系式知命々)=^1(Z)成立，其中为第n次本征模态函数分解中经过第P^i次筛选后的剩余函数 步骤23、根据筛选程序获取输入的待处理曲线％(/)经过第Η次本征模态函数分解的 剩余的残差函数中经过第次筛选后的剩余函数KO ；
步骤24、采用标准偏差准则判断输入的待处理曲线经过第Η次本征模态函数分 解的剩余的残差函数中经过第Jt次筛选后的剩余函数KO是否满足本征模态函数的条件，
BP(^ft-D(O - Kt(0Y /是否小 Fsd , 0.2</J5d <0.3 ；
判断结果为是，执行步骤25，判断结果为否，则t = t + l，然后执行步骤23， 步骤25、提取一个本征模态函数分= Htjo(I)；
步骤26、获取输入的待处理曲线经过第 次本征模态函数分解的剩余的残差函
数W) ω- ω ；
步骤27、判断是否满足下述关系式成立2 ；
判断结果为否，则η 二 η+1’然后执行步骤22，判断结果为是，完成提取过程，获得两 个本征模态函数分量一阶IMF分量、二阶IMF分.V.jq(/)，C2(T)^ Jl1个经过第2次本
征模态函数分解的剩余的残差RES :r2(i)。即祐)=^/)^ +*^)
步骤三所述的获取二阶IMF分量的幅值和瞬时频率的过程为
步骤31、对二阶IMF分量进行离散褶积，获得其希尔伯特变换；
步骤32、获取二阶IMF分量的幅值C2(I)的解析信号为e2(t)+jy(f)，将所述解析
信号的包络振幅作为二阶IMF分量的幅值，解析信号的包络振幅按如下公式 计算
步骤33、获取二阶IMF分量的瞬时频率/( )先求取所述解析信号的相角然后，根据下式获取二阶IMF分量的瞬时频率/(O
1 d
2π Μ合成高光谱曲线记录了同一目标位置像素的连续采样光谱序列，如果目标发生自 旋，则该序列会存在循环的光谱。但是同样环境下采集的光谱曲线是具有高度相关性的，只 在一部分高频信息所在波段存在差异，因此自适应地提取高频信息所在波段是非常有意义 的，可以用这些波段的集合作为各采样时刻的特征位置，下面在步骤四中详细说明。步骤四所述的特征波段集合的获取过程为
步骤41、计算二阶IMF分量的幅值平均值5和瞬时频率的平均值J" α = € ΛΓ},其中# = ^xr，
/ =IBBOTJ {/(/) I 1 </<
步骤42、舍弃二阶IMF分量中幅值小于幅值平均值《 —半的波段，同时舍弃瞬时频率 小于瞬时频率平均值Ji ,获取特征时刻集合ψ
f(t)>f
Ψ
a , 1</<JV
丨> —
步骤43、在范围内依次截取长度为的集合，获得Γ段采样数据各自对应的特征波段集
A W WW
口 I 1，I 2 …3 ι JT O如果对象为空间碎片，在自旋一周之后，相同位置的观察像素会具有一定相关性， 即便观察像素在碎片上的位置稍有偏移也是如此，只是相关性有所减弱。下面在步骤五进 行详细说明。步骤五所述的循环搜索特征波段集合，判断目标是否旋转，进而确定疑似目标是 否为空间碎片的过程为
步骤51、计算特征波段集合Y13Y2^Y1■两两之间特征集合的相异波段数目 D9 1< < j<T ,
步骤52、计算所述相异波段数目的均值
步骤53、计算等间隔采样数据之间的相异波段数目采用下式计算均值Ifl jAn ·—·· ΤΚ ΕΙΙ [d§ I ij MPS 为i，1 < i < j < r}3 0.1Γ <i < 0.9Γ ,
10由于采样数据量为Γ ,循环周期的有效数值经验证应在0.1Γ到0.9Γ之间。步骤54、获取Mjfc中最小值所对应的间隔jt* k* = aigmiii{Mt | ΟΛΓ <k< 0.9Γ}，
在0.1Γ到CXOtt^、日H取多个k值，i十算 每个k值对应白勺均值Mk ,并将其中最小
均值找出来，这个最小均值对应的间隔Jt*是用于衡量相关性的一个周期参数。步骤55、判断M/miii(Mt)>2足否成立，
判断结果为是，断定目标旋转导致观测像素在间隔Γ的周期具有较高的光谱相关性， 使得特征集合相异波段数目减少，进而判断疑似目标为空间碎片；判断结果为否，判断疑似 目标不是空间碎片。
具体实施方式
二、下面结图1至图9，结合某空间碎片模型的高光谱旋转成像序列 给出一个具体实施例
执行步骤一选定观察对象中心位置处目标像素的多个连续高光谱数据序列进行数据 合成，形成待处理曲线。对连续采集的20个高光谱数据，调用目标中心位置像素的高光谱曲线(共计20 段，每段长度m为220)，接着对每段高光谱曲线进行左右翻转，并连接在原始光谱曲线之 后，形成新的高光谱曲线(共计20段，每段长度为2m=440)，最后将20段新的光谱曲线依次 连接得到一条合成曲线为待处理曲线(如图3所示，长度==8800)作为后续步骤 的一维输入光谱数据——待处理曲线。执行步骤二 对合成的一维输入光谱数据进行一维经验模态分解，分解至一阶、二 阶本征模态函数以及残差函数即停止分解，具体描述如下
设定输入的待处理曲线为,/=1,2,.. ,8800，
步骤21、IMF分解过程初始化 = ,且满足关系式？^= 成立，其中 ‘ifr)为第(Λ—Ι)次分解后剩余的残差函数
步骤22、筛选过程初始化， =1，且满足关系式Pn(H)CO = ^1(Z)成立，其中
WO为第η次本征模态函数分解中经过第P^5次筛选后的剩余函数
步骤23、根据筛选程序获取输入的待处理曲线经过第π次本征模态函数分解的 剩余的残差函数中经过第it次筛选后的剩余函数KO ；
步骤24、采用标准偏差准则判断输入的待处理曲线经过第H次本征模态函数分 解的剩余的残差函数中经过第次筛选后的剩余函数Κ )是否满足本征模态函数的条件，即— Kt(O)2 /^fc-D(0 是否小于阈值Uso , Hsd 二0.25 ；
判断结果为是，执行步骤25，判断结果为否，则fc = t + l，然后执行步骤23， 步骤25、提取一个本征模态函数分= AntCO ；
步骤26、获取输入的待处理曲线经过第 次本征模态函数分解的剩余的残差函
数 C W = I1W- W ；
步骤27、判断是否满足下述关系式成立If >2 ；
判断结果为否，则《二《+1，然后执行步骤22，判断结果为是，完成提取过程，获得两 个本征模态函数分量一阶IMF分量(见图4)、二阶IMF分量(见图5) ;和j个经过第2次本征模态函数分解的剩余的残差RES (见图6) Γ2( )。即执行步骤三对IMF2进行希尔伯特变换，得到二阶IMF分量的幅值和瞬时频率。首先，按如下算式对二阶IMF分量c2(t% t = 1,2,...,8800进行离散褶积得到其
希尔伯特变换·Κ ) m
接着，按如下算式计算q 的解析信号巧的包络振幅《^ 然后，按如下算式计算解析信号的相角_和瞬时频率/( )
^(O=HTCfeOli ^ 1,
UjCOJ
rn-^^m
λπ at
该步骤完毕，得到二阶IMF分量的幅值《Θ和瞬时频率/(*)，分别如图7和图8所示。执行步骤四综合判断高光谱数据的高频信息所在波段，得到特征波段集合。对于二阶IMF分量的瞬时频率中小于其平均值=KMHSOOj )的波
段予以排除以保留具有特征意义的高频信息所在波段；同时由于二阶IMF分量的幅值和瞬 时频率是依波段对应的，而过小的幅值其相应波段以噪声成分居多，因此对于幅值小于其
平均值(S=WeaB二分之一的波段也要予以排除，最后得到特征时刻集
12
然后对IiiSSSOO依次截取长度为440的集合，即得到20段采样数据各自对应的特 征波段集〈ΥΨρΨ^^Ψ ,执行步骤五循环搜索判断目标是否旋转，进而确定是否为空间碎片。首先，计算Ψ^Ψ^,Ψ 两两之间特征集合的相异波段数目Iif. 1< </<20 ；
接着，计算相异波段数目的均值应二 mean |l<i<j<20} 47—79 ；
由于采样数据量为r=20,循环周期的有效数值经验证应在ojt=2到ο」9Γ=ι·之 间，对等间隔采样数据之间的相异波段数目采用下式计算均值 M1 =meaii{l3L| I5/间隔为λ, 1<I<J<20}, 2<i:<lS ；
找出Mfc (见表1)中最小值所对应的间隔JfJtSlg=IO ；
由于MZnda(Mi) = 47.79/1630=253大于阈值2，因此可判断目标旋转导致观测
像素在间隔10的周期具有较高的光谱相关性，使得特征集合相异波段数目减少，进而判断 目标为空间碎片。表1等间隔采样数据之间的相异波段数目均值
权利要求
基于希尔伯特 黄变换的空间碎片高光谱序列探测方法，其特征在于，它包括如下步骤步骤一、连续T次采样疑似目标中心位置的高光谱曲线，并将获取的T段高光谱曲线处理合成为待处理曲线；步骤二、对步骤一获取的待处理曲线数据进行一维经验模态分解，获取两个本征模态函数分量和残差，两个本征模态函数分量为一阶IMF分量和二阶IMF分量；步骤三、对所述二阶IMF分量进行希尔伯特变换，得到二阶IMF分量的幅值和瞬时频率；步骤四、保留同时满足幅值高于均值一半和瞬时频率高于均值两个条件的二阶IMF分量，并分割形成特征波段集合；步骤五、循环搜索特征波段集合，判断目标是否旋转，进而确定疑似目标是否为空间碎片，如旋转，则确定疑似目标为空间碎片，如不旋转，则确定疑似目标不是空间碎片。
2.根据权利要求1所述的基于希尔伯特-黄变换的空间碎片高光谱序列探测方法，其 特征在于，步骤一中的连续T次采样疑似目标中心位置的高光谱曲线，并将获取的T段高光 谱曲线处理合成为待处理曲线的过程为步骤11、采用高光谱成像仪采集观测区域的多幅连续高光谱图像，选择其中疑似目标 的中心位置像素，获取一条疑似目标中心位置的高光谱曲线，所述疑似目标中心位置的高 光谱曲线的长度为m，步骤12、按照步骤11的方法连续T次采样疑似目标中心位置的高光谱曲线，获取具有 连续数据的T段高光谱曲线，步骤13、对每段高光谱曲线进行左右翻转，并连接在原始高光谱曲线之后，形成T段长 度为2m的新的高光谱曲线，步骤14、将所述T段长度为2m的新的高光谱曲线依次拼接，合成为待处理曲线。
3.根据权利要求1所述的基于希尔伯特_黄变换的空间碎片高光谱序列探测方法，其 特征在于，步骤二中的获取两个本征模态函数分量和残差的过程为设定 输入的 待处理 曲线为 步骤21、IMF分解过程初始化η = 1 ,且满足关系式 成立，其中 次分解后剩余的残差函数；步骤22、筛选过程初始化，k=r且满足关系 成立，其中 为第η次本征模态函数分解中经过第(t-5次筛选后的剩余函数步骤23、根据筛选程序获取输入的待处理曲线经过第Η次本征模态函数分解的 剩余的残差函数中经过第Jt次筛选后的剩余函数KO ；步骤24、采用标准偏差准则判断输入的待处理曲线经过第H次本征模态函数分 解的剩余的残差函数中经过第Jt次筛选后的剩余函数Kk(J)足W满足本征模态函数的条 是否小于阈值丑迎 判断结果为是，执行步骤25，判断结果为否，则 然后执行步骤23， 步骤25、提取一个本征模态函数分 ；步骤26、获取输入的待处理曲线J<f)经过第n次本征模态函数分解的剩余的残差函 数 步骤27、判断是否满足下述关系式成立H> 2 ；判断结果为否，则 然后执行步骤22，判断结果为是，完成提取过程，获得两 个本征模态函数分量一阶IMF分量、二阶IMF分量{c丨(f_)，C2(f_；^ ；和个经过第2次本征模态函数分解的剩余的残差RES Γ2( )。
4.根据权利要求3所述的基于希尔伯特_黄变换的空间碎片高光谱序列探测方法，其 特征在于，步骤24的中丑皿=0.25。
5.根据权利要求1所述的基于希尔伯特-黄变换的空间碎片高光谱序列探测方法，其 特征在于，步骤三所述的获取二阶IMF分量的幅值和瞬时频率的过程为步骤31、对二阶IMF分量进行离散褶积，获得其希尔伯特变换Jft0 ；步骤32、获取二阶IMF分量的幅值C2(t)的解析信号为C2^) + Mf)，将所述解析信号的包络振幅φ)作为二阶IMF分量的幅值，解析信号的包络振幅按如下公式 计算步骤33、获取二阶IMF分量的瞬时频率/(if)先求取所述解析信号的相角 然后，根据下式获取二阶IMF分量的瞬时频率
6.根据权利要求1所述的基于希尔伯特-黄变换的空间碎片高光谱序列探测方法，其 特征在于，步骤四所述的特征波段集合的获取过程为步骤41、计算二阶IMF分量的幅值平均值5和瞬时频率的平均值J 步骤42、舍弃二阶IMF分量中幅值小于幅值平均值α —半的波段，同时舍弃瞬时频率 小于瞬时频率平均值7，获取特征时刻集八ψ步骤43、在1 Sf 范围内依次截取长度为的集合，获得段采样数据各自对应的特征波段集合
7.根据权利要求1所述的基于希尔伯特-黄变换的空间碎片高光谱序列探测方法，其 特征在于，步骤五所述的循环搜索特征波段集合，判断目标是否旋转，进而确定疑似目标是 否为空间碎片的过程为步骤51、计算特征波段集合 ·两两之间特征集合的相异波段数目 D9 \<i<j<T ,步骤52、计算所述相异波段数目的均值 步骤53、计算等间隔采样数据之间的相异波段数目采用下式计算均值 步骤54、获取Mjfc中最小值所对应的间隔 步骤55、判断 是否成立，判断结果为是，断定目标旋转导致观测像素在间「Rt*的周期具有较高的光谱相关性， 使得特征集合相异波段数目减少，进而判断疑似目标为空间碎片；判断结果为否，判断疑似 目标不是空间碎片。
全文摘要
基于希尔伯特-黄变换的空间碎片高光谱序列探测方法，属于图像处理领域，本发明为解决采用传统的危险空间碎片探测方法时，必须依赖有关空间碎片图像的样本信息，算法适应性差的问题，本发明方法包括一、连续T次采样疑似目标中心位置的高光谱曲线，并将获取的T段高光谱曲线处理合成为待处理曲线；二、一维经验模态分解，获取两个本征模态函数分量；三、对二阶IMF分量进行希尔伯特变换，得到其幅值和瞬时频率；四、保留幅值高于均值一半的部分，保留瞬时频率高于均值的部分，并分割形成特征波段集合；五、循环搜索特征波段集合，判断目标是否旋转，如旋转，则确定疑似目标为空间碎片，如不旋转，则确定疑似目标不是空间碎片。
文档编号G06T7/00GK101916439SQ20101023697
公开日2010年12月15日 申请日期2010年7月27日 优先权日2010年7月27日
发明者张淼, 沈毅, 王强, 王艳 申请人:哈尔滨工业大学