基于动态散斑时频特性的合作目标微动参数反演

王利国; 李亚清; 巩蕾; 王谦

doi:doi:10.3788/AOS202141.0228001

光学学报, 2021, 41 (2): 0228001, 网络出版: 2021-02-27

基于动态散斑时频特性的合作目标微动参数反演下载： 839次

Inversion Algorithm for Micro-Motion Parameters of a Cooperative Target Based on Time-Frequency Feature of Dynamic Speckle

论文大纲

王利国 ^*李亚清巩蕾王谦

作者单位

西安工业大学光电工程学院, 陕西西安 710021

AI 词云图 AI语音精读 AI语音超短摘要

注：本部分内容由 AI 自动生成，请您知悉。

摘要

基于激光动态散斑的时频信号,建立了一套适用于角反射器阵列的微动参数反演算法。首先利用物理光学方法推导了角反射器及角反射器阵列激光散斑的实时强度公式,然后基于短时傅里叶变换研究了散斑功率谱形成机制及其数字特征,最后提出了频谱相关法及时频-相幅变换算法,并采用该方法提取了动态散斑时频谱线周期及振幅分布,反演了三种典型运动状态下目标的自旋周期及旋转轴指向。结果表明:基于几个周期的散斑强度序列,所提反演算法可以得到高精度的旋转周期及视线角,但对旋转轴方位角的反演精度相对较差,需要更多的观测数据才能得到满意的结果。

Abstract

In this paper, a set of inversion algorithms for the micro-motion parameters in a retroreflector array was established based on the time-frequency signals of laser dynamic speckles. Firstly, the real-time intensitys formula of laser speckle from a retroreflector and a retroreflector array were derived by a physical optics method. Then, based on the short-time Fourier transform, the formation mechanism and digital characteristics of the speckle power spectrum were studied. Finally, the spectral correlation method and time-frequency to phase-amplitude transform algorithm were proposed, and the algorithm was used to extract the period and amplitude distribution of the time-frequency spectral lines of dynamic speckles. Besides, we made an inversion on the spin period and the direction of the rotation axis of a target in three typical motion states. The results show that the proposed inversion algorithm can obtain the rotation period and line-of-sight angle with high accuracy through speckle intensity signals in a few periods, but the inversion accuracy for the azimuth angle of the rotation axis is relatively poor. It follows that more observation data are required for satisfactory results.

1 引言

对空间目标的运动状态进行观测是航天器测控的重要组成部分,也是目标识别的基本手段之一,其中以目标微动状态的识别和测量最为关键,也最为困难。这是因为目标距离遥远,表面细节辨认困难,无法通过一般的成像技术观察目标自身的旋转、翻滚等运动状态。

目前,基于地基观测的目标微动状态测量技术主要有三种,分别是基于光变曲线、微多普勒频谱以及激光动态散斑的测量技术。

光变曲线反映的是目标姿态变化导致入射角、散射角发生变化,从而导致散射光强度的变化。基于接收到的光强变化的周期性可以反演卫星、火箭等空间目标的旋转周期和旋转轴指向,这一方法在理论和实践上都比较成熟^[1-7],但不适用于散射光强对目标转动不敏感的情况,而且测量精度也不够理想。

目标表面各点运动速度的差异会导致各点回波多普勒频移不同,从而形成多普勒频谱。基于这一原理的微动反演技术被广泛用于直升飞机旋翼转速的判断以及导弹目标识别等^[8-10]。其优点是反演精度高,信息量大,缺点是光路复杂不易调整,对硬件要求较高。目前,该技术多用于中程目标的探测,对于空间目标等远程目标的探测尚处于理论研究和实验探索阶段^[11-13]。

目标的运动会导致接收点上激光散斑强度发生变化,动态散斑的变化频率比普通光变曲线的频率更高,携带的信息也更精确。激光散斑常被用于测量平动、转动和振动等多种运动形式的速度和频率。具体测量方法又分为两类:一类是基于时间相关函数的测量方法,这类方法常用于测量圆柱旋转角度^[14]及速度^[15]、球半径及旋转速度^[16]、人的心率^[17]等;第二类是基于动态散斑强度频谱的测量方法,但这类方法多用来测量流体的速度,如利用散斑对比度成像方法测量血流速度^[18-20]。此外,还有直接使用动态散斑功率谱的方法,如:1985年,Fujii等^[21]开发了一种利用动态散斑功率谱斜率测量血流平均速度的技术;2005年,Barton等^[22]在光学相干层析成像(OTC)系统上,利用动态散斑频谱宽度测量了眼底的血流速度。目前,研究人员对动态散斑形成机制还缺少深入理解^[23],而上述方法一方面依赖于具体的目标形态及其运动状态,灵活性不足,另一方面不能全面挖掘散斑携带的信息,可测物理量有限,测量过程中有些重要参数只能通过实验得到。

本文以自旋运动的合作目标作为研究对象,提出了一种基于动态散斑功率谱的目标微动参数反演方法。首先推导了动态散斑功率谱公式,该公式提示了散斑功率谱的形成机制;然后参考基于微多普勒的微动反演算法,利用短时傅里叶变换得到阵列回波时频信号;之后提出频谱相关函数及时频-相幅变换算法,并采用该算法提取时频信号周期及幅度等信息;最后以三种典型运动状态下的角反射器阵列为例,反演了目标的转动周期和旋转轴指向角度。

2 角反射器阵列激光动态散斑模型

本节以图1(a)所示的角反射器阵列为例建立动态合作目标的激光回波模型。图1(a)是一个典型的半球形角反射器阵列,它曾装载于日本的ADEOS-II卫星上^[24]。它共有9个反射单元,单元孔径中心分布于半径为8.57 cm的半球上。中心为9号单元(其法线为阵列中心轴),1~8号单元均匀分布在阵列周围,相邻单元的法线在水平面内的夹角为45°,每个单元的法线与9号单元法线之间的夹角均为50°。

图 1. 角反射器阵列结构及探测系统光路几何示意图。(a)角反射器阵列结构;(b)探测系统光路几何

Fig. 1. Structure of retroreflector array and optical geometry of detecting system. (a) Structure of retroreflector array; (b) optical geometry of detecting system

下载图片查看所有图片

目标位于空中,并以角速度ω绕某固定轴(此轴称为旋转轴)旋转,利用单站激光雷达对其进行观测。目标与激光雷达的几何关系如图1(b)所示。L为观测距离,θ为视线方向与旋转轴的夹角(视线角),α为阵列中心轴与旋转轴的夹角。此外,定义φ为旋转轴在xy平面内的投影方向与x轴的夹角(称为指向方位角),旋转轴在阵列坐标系x'y'z'中的方位角为φ'(此角度不随角反射器的旋转而改变)。以上4个角度可唯一确定角反射器阵列的姿态。

2.1 角反射器单元激光回波模型

本文研究的角反射器阵列所用的角反射器单元为熔石英材质,采用内切圆切割而成,折射率为n=1.455,底面半径为r=1.586 cm。激光以任意入射角i₀照射角反射器,其后向散射截面可以表示为^[25]

σ (i_{0}) = βS (i_{0}), (1)

其中,

S (i_{0}) = 2 r^{2} [\arcsin \sqrt[]{1 - 2 ta n^{2} (i)} - \sqrt[]{2} \tan (i) \sqrt[]{1 - 2 ta n^{2} (i)}], (2)

式中:β是一个系数,受激光波长、角反射器反射效率及接收点位置的影响;S(i₀)是角反射器的有效反射面积^[26];i=arcsin[sin(i₀)/n]。在单站模式下,可认为β是一个时不变的值。

根据散射截面的定义,角反射器阵列上第j个反射单元的回波远场幅度可以表示为

|U_{j} (p)| = σ_{j} |U_{i} (ρ_{j})| / L, (3)

式中:p和ρ_j分别为接收点和角反射器中心的三维坐标;σ_j为第j个反射单元的后向散射截面;U_i(ρ_j)为入射光场。本文采用高斯波束,故角反射器孔径上的入射光场可以表示为

U_{i} (ρ) = \exp (- \frac{ρ^{2}}{w^{2}} + i k \frac{ρ^{2}}{2 F} + i kL), (4)

式中:ρ为入射点的横向坐标;ρ为入射点的横向距离;k=2π/λ为波数,λ为波长;w和F分别为入射波束的波束半径和波前曲率半径。

当目标距激光雷达的距离L很大时,发射场和反射场的波前都近似为球面,其波前曲率半径约等于L,此时角反射器的散射场可以表示为

U_{j} (p, t) = \frac{σ_{j}}{L} \exp (- \frac{ρ_{j}^{2}}{w^{2}}) \cdot \exp [i kL (ρ_{j}, 0, t) + i kL (ρ_{j}, p, t) + i ϕ_{j}], (5)

式中:ρ_j为第j个反射单元孔径中心的横向距离;L(ρ_j,0,t)和L(ρ_j,p,t)分别为t时刻第j个反射单元孔径中心到发射点和接收点的距离;ϕ_j是由角反射器单元安装误差引起的附加随机相位。

当目标处于运动状态时,反射器单元的回波幅度和相位都处于变化之中,但幅度相对于相位是缓变函数。因此,在短时间内,本文只考虑运动对相位的影响。设角反射器单元的运动速度为v_j,在时刻t=t₀+Δt时,激光双程传输的距离为

\{\begin{array}{l} L (ρ_{j}, 0, t) = L (ρ_{j}, 0, t_{0}) + v_{j 1} Δ t \\ L (ρ_{j}, p, t) = L (ρ_{j}, p, t_{0}) + v_{j 2} Δ t \end{array}, (6)

式中:L(ρ_j,0,t₀)和L(ρ_j,p,t₀)分别为t₀时刻第j个反射单元孔径中心到发射点和接收点的距离;v_j₁和v_j₂为角反射器速度的径向分量。对于单站激光雷达,有v_j₁=v_j₂=v_j。此时,回波场可以表示为

U_{j} (p, t) = \frac{σ_{j}}{L} \exp (- \frac{ρ_{j}^{2}}{w^{2}}) \exp [i kL (ρ_{j}, 0, t_{0}) + i kL (ρ_{j}, p, t_{0}) + i ϕ_{j}] \exp (i 2 k v_{j} Δ t) 。 (7)

令

\{\begin{array}{l} E_{j} = \frac{σ_{j}}{L} \exp (- \frac{ρ_{j}^{2}}{w^{2}}) \\ Θ_{j} = kL (ρ_{j}, 0, t_{0}) + kL (ρ_{j}, p, t_{0}) + ϕ_{j} \\ f_{j} = 2 v_{j} / λ \end{array}, (8)

式中:E_j、Θ_j和f_j分别为角反射器的回波幅度、初始相位和多普勒频移;λ为激光波长。

2.2 角反射器阵列动态散斑强度模型

总接收场由各反射单元的回波相干叠加而成,即

U (p, t) = \overset{m}{\sum_{j = 1}} U_{j} (p, t) = \overset{m}{\sum_{j = 1}} E_{j} \exp (Θ_{j}) \exp (i 2 π f_{j} Δ t), (9)

则接收场的散斑强度为

I (p, t) = {|U (p, t)|}^{2} = \overset{m}{\sum_{j = 1}} \overset{m}{\sum_{j' = 1}} I_{jj'} \exp (i Φ_{jj'}) \exp [i 2 π (f_{j} - f_{j'}) Δ t], (10)

式中:m为角反射器单元的总数;I_jj'=E_jE_j'和Φ_jj'=Θ_j-Θ_j'分别为两反射单元回波的互强度和初始相位差;Φ_jj'在j=j'时是常量0,否则是一个均方根远大于波长的随机变量。

由于Φ_jj'的存在,接收到的光强是一个随时间变化的随机变量。图2为数值仿真得到的旋转角反射器阵列回波强度的序列样本。可见,接收的光强I(0,t)随着时间的变化有着高度的随机性,无法直接从强度序列中获得有效信息。

图 2. 自旋运动角反射器阵列散斑强度的时间序列

Fig. 2. Intensity time sequence of speckle field of a rotating retroreflector

下载图片查看所有图片

2.3 角反射器阵列动态散斑时频信号模型

为了在接收光强时得到目标运动状态的信息,本文对接收光强时间序列作傅里叶变换,得到了其频谱函数。频谱函数表达式为

G (p, f) = F [I (p, t)] = \overset{m}{\sum_{j = 1}} \overset{m}{\sum_{j' = 1}} I_{jj'} \exp [i Φ_{jj'}] δ [f - (f_{j} - f_{j'})], (11)

式中:δ为狄拉克函数;f为光强信号的频率。G(p,f)在f≠0时仍是一个随机变量。进一步计算频谱强度,得

G_{2} (p, f) = G (p, f) G^{*} (p, f) = \overset{m}{\sum_{j = 1}} \overset{m}{\sum_{j' = 1}} \overset{m}{\sum_{l = 1}} \overset{m}{\sum_{l' = 1}} I_{jj'} I_{ll'} \exp [i Φ_{jj'} - i Φ_{ll'}] δ [f - (f_{j} - f_{j'})] δ [f - (f_{l} - f_{l'})] = \overset{m}{\sum_{j = 1}} \overset{m}{\sum_{j' = 1}} \overset{m}{\sum_{l = 1}} \overset{m}{\sum_{l' = 1}} I_{jj'} I_{ll'} \exp [i Φ_{jj'} - i Φ_{ll'}] δ [(f_{j} - f_{j'}) - (f_{l} - f_{l'})] δ [f - (f_{j} - f_{j'})], (12)

式中:G₂(p,f)为散斑功率谱密度函数。由(12)式可知,G₂(p,f)只在f=f_j-f_j'且f_j-f_j'=f_l-f_l'时非零,或者说,非零值的位置是由任意两反射单元回波多普勒频移之差决定的。由于散射单元的离散性,功率谱并不是连续的,而是存在有限数量的频点。假设功率谱内共存在G个非零频点,其中第g个频点对应的频差为Δf_g=f_j-f_j'。由于阵列的对称性等因素,可能存在N_g(N_g≥1)对(j,j')组合对应同一个Δf_g(实际中,当两个频差之差的绝对值小于频谱分辨率时,就可以认为它们相同)。

当N_g=1时,功率谱在f=Δf_g= $f_{j_{g}}$ - $f_{j'_{g}}$ 处的值为

K_{g} = G_{2} (p, f = Δ f_{g}) = I_{j_{g} j'_{g}} I_{j_{g} j'_{g}} \exp [i Φ_{j_{g} j'_{g}} - i Φ_{j_{g} j'_{g}}] = I_{j_{g} j'_{g}}^{2} \geq 0 。 (13)

当N_g>1时,功率谱的值为

K_{g} = G_{2} (p, f = Δ f_{g}) = \overset{N_{g}}{\sum_{n = 1}} \overset{N_{g}}{\sum_{n' = 1}} I_{j_{gn} j'_{gn}} I_{j_{gn'} j'_{gn'}} \exp [i Φ_{j_{gn} j'_{gn}} - i Φ_{j_{gn'} j'_{gn'}}] = \overset{N_{g}}{\sum_{n = 1}} \overset{N_{g}}{\sum_{n' = 1}} I_{j_{gn} j'_{gn}} I_{j_{gn'} j'_{gn'}} \exp [i k Φ_{j_{gn} j'_{gn}} - i k Φ_{j_{gn'} j'_{gn'}}] \geq 0, (14)

式中:j_gnj'_gn∈{j_g₁j'_g₁,j_g₂j'_g₂,…,j_gNgj'_gNg}。

角反射器阵列回波散斑功率谱密度函数可以表示为

G_{2} (p, f) = \overset{G}{\sum_{g = 1}} K_{g} δ (f - Δ f_{g}) 。 (15)

综上,散斑实时功率谱上存在G个非零离散频点,每个频点的强度都大于或等于零。由于功率谱来自强度的傅里叶变换,所以它是一个偶函数。

由于一次傅里叶变换的时间序列不能太长,因此需要在不同的时刻对目标进行探测,以全面反映目标的运动信息。从信号处理角度,利用时间窗口滑动截取散斑强度信号,作短时傅里叶变换,就可得到动态散斑的时频信号G₂(p,f,t)。

3 典型运动状态下角反射器阵列的微动参数反演

时频信号虽然仍含有一定的随机性,但已经可以提供足够的信息。本节将展示如何由时频信号反演得到目标的运动参数。反演所需的散斑强度序列由计算机程序依据(10)式数值仿真得到。

3.1 同轴旋转阵列

首先考虑最简单的情况,即阵列旋转轴与阵列中心轴共轴。设置仿真参数如下:波长λ=1.06 μm,角反射器阵列到观察点的距离L=1000 km,旋转周期T=2 s,视线角θ=5°,采样周期dt=10 μs,总仿真时长t_e=30 s,短时傅里叶变换所用矩形时间窗长度N=1024。

对强度序列作短时傅里叶变换,可以得到接收强度的时频信号,如图3所示。离散的频谱随着时间变化呈周期性变化。根据谱线的幅度,可以将它们划分为5组。其中强度最大的8条谱线分别对应着中心反射单元与周围8个反射单元的组合,本文称之为0阶谱线。其他4组分别由间距不同的周围反射单元反射光干涉形成,其谱线数目不便计数,依幅度从小到大称为1~4阶谱线。

图 3. 同轴旋转阵列动态散斑的时频图

Fig. 3. Time-frequency feature of dynamic speckle of coaxial rotating array

下载图片查看所有图片

为了更准确地得到谱线周期,本文对(16)式所示的频谱相关函数进行了计算,计算结果如图4所示。

图 4. 同轴旋转阵列动态散斑的频谱相关函数

Fig. 4. Spectral correlation function of dynamic speckle of coaxial rotating array

下载图片查看所有图片

b (Δ t) = \int_{0}^{t_{e}} d t \int_{-}^{\infty} d f G_{2} (p, f, t) G_{2} (p, f, t + Δ t) 。 (16)

通过提取曲线的峰值位置可以得到曲线具有两种周期性:一种周期为2 s,其对应最大相关峰值的位置,通过最大峰值的位置可以反演目标旋转周期;另一种周期为0.25 s,即旋转周期的1/8,这是由角反射器阵列的旋转对称性带来的。但由于每个反射单元的安装误差是随机的,角反射器阵列的旋转对称性并不完美,所以旋转1/8周后的功率谱并不完全一致,其相关系数小于1。通过仿真数据反演的周期没有显示出存在误差,但考虑到有限采样时间对时间分辨率的限制,可以认为旋转周期的反演精度即为采样时间,在本算例中反演精度为10^-5 s(2009年,Kirchner等^[1]基于角反射器阵列光变曲线反演旋转周期的精度只能达到约0.6 s,相对精度为1%)。

在转速已确定的情况下,一条谱线的幅度由对应的两反射单元的相对位置决定,第g条谱线的实时频率可以表示为

Δ f_{g} (t) = 2 |v_{gj} (t) \cdot e_{j} - v_{gj'} (t) \cdot e_{j'}| / λ = 2 ω r_{g} \sin θ \cos (ωt + ϕ_{g}) / λ, (17)

其中:e_j和e_j'分别是形成第g条谱线的一对反射单元j和j'所在的径向方向,对于遥远的空间目标,可近似认为e_j=e_j';ω=2π/T为阵列旋转的角速度;r_g= $|r_{g}|$ 为单元j和j'之间的横向距离,r_g为两反射中心矢量差垂直于旋转轴的分量;ϕ_g为第g条谱线的相位,由相应反射单元对的相对位置决定。当r_g垂直于入射面(旋转轴和视线轴所在的平面)时,频差达到最大值,即

F_{g} = 2 ω r_{g} \sin θ / λ 。 (18)

为了从仿真数据中提取角反射器阵列的指向角θ,需要从图3中得到各阶谱线的幅度,但直接判读会造成较大误差。考虑到谱线在时频图中的轨迹是余弦函数,因此把时频图沿时间向累加,累加和在谱线零斜率点(即极值点)形成峰值,提取所有正半轴的峰值位置,可得到各阶谱线的振幅。

利用公式

Γ (f) = \int_{0}^{t_{e}} G_{2} (p, f, t) d t, (19)

可以得到图5所示的频谱累积强度曲线。(19)式中,Γ(f)为频谱累积强度。提取曲线的峰值位置,结合(18)式可以反推得到角反射阵列的视线角,如表1所示。结果显示,依据单谱线反演的视线角精度在1%以内,误差主要来源于图5中曲线峰值位置的定位误差。综合多个曲线的反演结果,得到视线角反演值为θ=4.98°,相对误差约为-0.4%。这一反演精度也远优于Kirchner等^[1]基于光变曲线的目标倾角的反演精度(约1.6°)。

图 5. 同轴旋转阵列动态散斑的频谱累积强度分布

Fig. 5. Spectral integrated intensity of dynamic speckle of the coaxial rotating array

下载图片查看所有图片

表 1. 各阶谱线的振幅值及其对应的视线角反演值

Table 1. Amplitude of each spectral lines and corresponding line-of-sight angles

g	r_g /cm	F_g /kHz	θ /(°)	Δθ /(°)
1	4.22	10.83	4.99	-0.01
2	5.51	14.06	4.96	-0.04
3	7.79	19.92	4.97	-0.03
4	10.18	26.07	4.98	-0.02
5	11.02	28.32	4.99	-0.01

查看所有表

3.2 异轴旋转阵列

目标并不一定总是绕阵列中心轴旋转,设旋转轴与阵列轴的夹角为α=20°,视线角仍为θ=5°,旋转轴在阵列坐标系中的方位角φ'=5°,其他参数不变,利用计算机仿真得到散斑强度序列,进一步计算可以得到异轴旋转阵列动态散斑的时频图,如图6所示。与图3相比,图6主要发生了三点变化:一是同阶谱线的幅度不再一致,二是同阶谱线的平均强度也产生了较大差异,三是同一条谱线在不同时刻的强度发生了明显变化。这些变化都是由异轴旋转过程造成的不对称性引起的。

图 6. 异轴旋转阵列动态散斑的时频图

Fig. 6. Time-frequency feature of dynamic speckle of non-coaxial rotating array

下载图片查看所有图片

图7给出了异轴旋转阵列功率谱的时间相关函数。可见:最大相关峰值仍出现在2 s处(由此可以得到旋转周期为2 s),其他各峰仍然存在且位置不变,但峰值明显变小(这意味着时频信号的1/8周期性受到了削弱)。可以预料,旋转轴与阵列轴的夹角越大,频谱1/8周期性越弱,非整数周期位置的相关峰值就会越小。

图 7. 异轴旋转阵列动态散斑的频谱相关函数

Fig. 7. Spectral correlation function of dynamic speckle of non-coaxial rotating array

下载图片查看所有图片

为了进一步反演角反射器阵列的旋转轴指向,需要提取各谱线的幅度。但原来振幅相同的同阶谱线发生了分裂,各谱线振幅值多且相互接近,若仍以频谱时间累加的方法,则各谱线对应的峰值区会交错重叠,导致峰值消失或峰值密集不易识别。因此,本文提出了一种新的变换算法——时频-相幅变换,其表达式为

γ (f', τ) = \int_{-}^{\infty} \int_{0}^{t_{e}} G_{2} (p, f, t) \cdot δ [f - f'h (t - τ)] d f d t, (20)

式中:h(t)为谱线的基函数,所有谱线随时间的变化都遵循此函数,只是每条曲线的振幅f'和时延τ不同。对于散射单元绕固定轴旋转的情况,谱线的基函数为

h (t) = \cos (\frac{2 π}{T} t), (21)

式中:T为阵列旋转周期。当f'、τ与时频图中第k条谱线的参数f'_k、τ_k相同时,振幅-相位函数γ(f',τ)在(f'_k,τ_k)处会出现峰值。提取出γ(f',τ)函数所有峰值的位置,就可以得到所有谱线的幅度。

在实际计算过程中需要把(20)式离散化。结合图6对应的时频数据,可以得到图8所示的振幅-相位函数图像,其纵坐标为Φ=2πτ/T,为相位延迟。图中存在一系列分散的峰值,提取出所有峰值的横坐标位置,就可以得到各谱线的幅度,如表2所示。

图 8. 异轴旋转阵列动态散斑时频谱线的振幅-相位分布图

Fig. 8. Amplitude-phase distribution of time-frequency spectral lines of dynamic speckle of non-coaxial rotating array

下载图片查看所有图片

表 2. 异轴旋转阵列时频谱线振幅的提取结果

Table 2. Extracted amplitude of time-frequency spectral lines of non-coaxial rotating array

g	f' /kHz
1	10.27
2	10.34
3	10.74
4	10.80
5	11.59
6	14.12
7	14.48
8	14.94
9	15.10
10	15.13
11	18.83
12	19.47
13	20.07
14	24.88
15	24.90
16	26.00
17	26.01
18	26.70
19	27.49
20	27.54
21	28.41

查看所有表

利用(18)式可以从理论上计算角反射器阵列在任意运动状态下的谱线振幅。记同一状态下所有振幅的理论值为集合Y,它是阵列姿态(θ,α,φ')的函数。将振幅-相位函数提取得到的所有振幅值记为集合Y'。二者之间的偏差定义为

δ = \sum_{g = 1} \min [Y' {(g) - Y]}^{2} 。 (22)

由(22)式可知δ也是阵列姿态的函数。利用粒子群优化算法搜寻到使δ最小化的(θ,α,φ'),即为反演得到的阵列姿态参数。本文基于表2中的数据进行了反演,反演结果为θ=5.01°,α=19.92°,φ'=3.44°。可见,在异轴旋转情况下,本文方法反演的θ、α的精度优于Kirchner等^[1]基于光变曲线方法反演的精度,但对方位角φ'的反演误差较大,这是因为阵列散斑功率谱对方位角较不敏感。

3.3 复合运动阵列

空间目标一般不会停在固定位置,除了绕轴旋转外,目标的质心位置也在作高速运动。在一个较短的时间段内,可以认为目标质心作平动。设目标质心以5000 m/s的速度(v)平行于x轴负方向移动,同时目标作绕阵列中心轴的旋转运动。激光雷达位于目标正下方,视线方向与旋转轴夹角仍为θ=5°,旋转轴方位角φ=0°,其他参数不变。

图9显示了复合运动角反射器阵列动态散斑时频图。与图3相比,图9最大的不同在于谱线的幅度随着时间推移而增大。平移运动从两个角度影响了散斑时频图。一是由于阵列上各反射单元相对激光雷达的方向略有不同,因此各单元平移运动在接收点处造成的多普勒频移也不同。这一效应造成的最大频差为f_max= $\frac{v r_{\max}}{Lλ}$ ,其中r_max是各反射单元之间的最大距离。对于本次仿真参数设置,f_max=0.5 kHz。即平动造成的频差远小于转动造成的频差,因此可以忽略这一因素的影响。平动对散斑时频信号的第二个影响是,阵列位置随时间的变化导致视线角发生了变化,从而导致谱线振幅随时间而变。图9中谱线幅度的增长是由第二个因素引起的。

图 9. 复合运动阵列动态散斑的时频图

Fig. 9. Time-frequency feature of dynamic speckle of the array with compound movement

下载图片查看所有图片

图10显示了复合运动阵列散斑时频信号的频谱相关函数。虽然时频信号的时间周期性被破坏,但代表1/8周期性和整数周期性的峰值仍然存在,只是峰值的整体趋势是随着时间差的增大而减小。结合角反射器阵列回波频谱相关函数的特性,仍然可以反演得到阵列旋转周期为2 s。

图 10. 复合运动阵列动态散斑的频谱相关函数

Fig. 10. Spectral correlation function of dynamic speckle of the array with compound movement

下载图片查看所有图片

考虑到本文的目的是要从散斑中提取阵列旋转运动的信息,而平动的信息可以由其他测量方式较为简单地测得,所以本文只利用动态散斑来反演旋转轴初始时刻的指向角度(θ₀,φ)。为此,仍需要提取各谱线的初始振幅。参考(18)式,当视线角θ发生变化时,谱线振幅可以表示为

f' (t) = F' \sin [θ (t)], (23)

式中:F'=2ωr_g/λ为谱线相应反射单元对造成的最大频差;视线角θ是时间的函数,且受观测系统多个因素的影响。

观察图9可以发现谱线幅度与观测时间大致呈线性关系,这意味着在观测时间不太长,视线角及谱线振幅变化不大的情况下,(23)式可以取一阶线性近似,即

f' (t) = F' [\sin (θ_{0}) + \cos (θ_{0}) Δ θ (t)] (24)

式中:θ₀为初始时刻的视线角;Δθ(t)=θ(t)-θ₀。在一阶近似下,与入射平面平行的目标速度分量会改变视线角,而垂直分量对视线角没有影响,即

Δ θ (t) = - vt \cos (Δ φ) / L, (25)

式中:Δφ=φ″-φ,φ″为目标平动方向的方位角。当目标沿x轴负向运动时,φ″=180°,因此本次仿真中Δφ的理论值为180°。

综合(24)式和(25)式可以得谱线幅度的近似表达式为

f' (t) = F' \sin (θ_{0}) - F' \cos (θ_{0}) \cos (Δ φ) vt / L 。 (26)

从图9所示的时频图中可以提取谱线上包络线的斜率,从而可以得到谱线振幅随时间的变化率s=0.114。谱线的数学表达式为

f' (t) = f'_{0} + f'_{0} \frac{s}{T} t 。 (27)

比较(26)和(27)式可得

Δ φ = \arccos [- \frac{sL}{vT} \tan (θ_{0})] 。 (28)

对于复合运动的目标,谱线的基函数为

h (t) = f' (t) \cos (\frac{2 π}{T} t) = (f'_{0} + f'_{0} \frac{s}{T} t) \cos (\frac{2 π}{T} t) 。 (29)

利用时频-相幅变换法得到了复合运动目标时频谱线的振幅-相位分布,如图11所示。可以发现,各谱线峰值点明显分为5列,分别对应同轴旋转阵列的5阶谱线。提取各峰值点位置,对同阶谱线峰值频率值求平均,可以得到各阶谱线的振幅。利用(18)式可反演初始时刻的视线角,结果如表3所示。对5个视线角的反演值取平均可得θ₀=5.01°,将其代入(28)式得Δφ=176.9°。比较仿真参数可知,视线角的反演精度较高,而方位角的反演精度较低。

图 11. 复合运动阵列散斑时频谱线的振幅-相位分布图

Fig. 11. Amplitude-phase distribution of time-frequency spectral lines of the array with compound movement

下载图片查看所有图片

表 3. 复合运动阵列时频谱线的振幅以及对应的反演视线角

Table 3. Amplitude of time-frequency spectral lines of the array with compound movement and corresponding inverted line-of-sight angle

g	F_g /kHz	θ₀ /(°)
1	10.87	5.00
2	14.24	5.02
3	20.08	5.01
4	26.27	5.01
5	28.43	5.01

查看所有表

4 结论

本文利用短时傅立叶变换方法对角反射器阵列激光回波时频信号特征进行了研究,揭示了动态散斑的形成机制及散斑功率谱的组成。针对角反射器阵列的三种典型运动状态,本文提出了基于动态散斑时频信号反演其旋转周期及旋转轴指向的方法。得到如下的具体结论:

1) 角反射器阵列动态散斑是由所有反射单元激光回波两两干涉叠加形成的,其功率谱是由各单元回波多普勒频移两两频差值组成的,功率谱随时间的变化形成了有限数目的谱线。

2) 当目标平动速度相对可以忽略时,旋转目标散斑时频谱线也是周期性振荡的,根据谱线的周期即可反演目标的旋转周期,其反演误差约为散斑信号的采样周期。即使考虑目标平动的影响,这一算法依然有效。

3) 根据各谱线的幅度及角反射器阵列的结构,可以反演阵列的旋转轴指向。对于质心不动的目标,根据单个探测点的散斑只能反演视线角。依靠多点探测或对平动目标进行连续探测,可以同时反演视线角与旋转轴方位角。所提反演算法对视线角的反演精度较高,达到0.1°以内,对方位角的反演精度较差。

4) 本文反演过程只利用了各谱线的幅度,但时频图中还有其他可利用的信息,包括各谱线的强度(亮度)和相位,综合利用这些信息,有望进一步提高反演精度。

本文提出的基于动态散斑时频信号的反演方法,具有设备简单、反演精度高的特点。这一方法不仅可以推广到其他类型的角反射器上,还可以推广到包含多个离散散射点的目标上。

参考文献

[1] Kirchner G, Kucharski D, Cristea E. Gravity probe-B: new methods to determine spin parameters from kHz SLR data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2009, 47(1): 370-375.

[2] Kucharski D, Kirchner G, Koidl F. Spin parameters of nanosatellite BLITS determined from Graz 2 kHz SLR data[J]. Advances in Space Research, 2011, 48(2): 343-348.

[3] Kucharski D, Lim H C, Kirchner G, et al. Spin axis precession of LARES measured by satellite laser ranging[J]. IEEE Geoscience and Remote Sensing Letters, 2014, 11(3): 646-650.

[4] 田琪琛, 李智, 徐灿, 等. 基于光学散射特性的失稳空间目标旋转分析[J]. 光散射学报, 2017, 29(3): 266-270.

Tian Q C, Li Z, Xu C, et al. Rotation analysis of unstable spaces based on optical scattering[J]. The Journal of Light Scattering, 2017, 29(3): 266-270.

[5] 刘通, 陈浩, 沈鸣, 等. 旋转卫星激光测距数据分析与处理[J]. 中国激光, 2017, 44(5): 0504001.

Liu T, Chen H, Shen M, et al. Spinning satellite laser ranging data analysis and processing[J]. Chinese Journal of Lasers, 2017, 44(5): 0504001.

[6] 刘通, 沈鸣, 高鹏骐, 等. 基于漫反射激光测距的火箭残骸翻滚姿态估算[J]. 中国激光, 2019, 46(1): 0104007.

Liu T, Shen M, Gao P Q, et al. Tumbling motion estimation of rocket body based on diffuse reflection laser ranging[J]. Chinese Journal of Lasers, 2019, 46(1): 0104007.

[7] 单斌, 梁勇奇, 李恒年. 基于光度观测的空间目标姿态与角速度估计[J]. 光学学报, 2017, 37(5): 0512002.

Shan B, Liang Y Q, Li H N. Attitude and angular speed estimation of spacial objects based on photometric observation[J]. Acta Optica Sinica, 2017, 37(5): 0512002.

[8] Setlur P, Ahmad F, Amin M. Helicopter radar return analysis: estimation and blade number selection[J]. Signal Processing, 2011, 91(6): 1409-1424.

[9] Clemente C, Soraghan J J. GNSS-based passive bistatic radar for micro-Doppler analysis of helicopter rotor blades[J]. IEEE Transactions on Aerospace and Electronic Systems, 2014, 50(1): 491-500.

[10] 王云鹏, 胡以华, 雷武虎, 等. 典型旋翼形状参数微多普勒激光探测计算方法[J]. 红外与激光工程, 2018, 47(9): 118-126.

Wang Y P, Hu Y H, Lei W H, et al. Algorithm of typical rotor shape parameters by micro-Doppler laser detection[J]. Infrared and Laser Engineering, 2018, 47(9): 118-126.

[11] 韩勋. 基于窄带微动特征的空间锥体目标识别方法研究[D]. 西安: 西安电子科技大学, 2015.

HanX. Research on recognition of space cone-shaped targets based on narrowband radar feature[D]. Xi'an:Xidian University, 2015.

[12] 王明军, 宫彦军, 柯熙政, 等. 地基观测空间在轨运行圆锥目标的激光多普勒频谱[J]. 中国科学:技术科学, 2018, 48(4): 424-432.

Wang M J, Gong Y J, Ke X Z, et al. The laser radar Doppler spectrum of moving conical target in its spaced orbit observed by ground-based station[J]. Scientia Sinica: Technologica, 2018, 48(4): 424-432.

[13] 郭力仁, 胡以华, 王云鹏, 等. 基于粒子滤波的高阶运动目标激光探测微动参数估计[J]. 光学学报, 2018, 38(9): 0912006.

Guo L R, Hu Y H, Wang Y P, et al. Laser micro-motion parameter estimation of high-order moving target based on improved particle filter[J]. Acta Optica Sinica, 2018, 38(9): 0912006.

[14] Hayashi A, Kitagawa Y. High-resolution rotation-angle measurement of a cylinder using speckle displacement detection[J]. Applied Optics, 1983, 22(22): 3520-3525.

[15] 张耿. 粗糙目标激光散斑统计特性及微运动特征分析[D]. 西安: 西安电子科技大学, 2013.

ZhangG. Statistical properties of laser speckle from rough objects and analysis on micro-motion characteristic[D]. Xi'an:Xidian University, 2013.

[16] Nakamura T, Asakura T. Statistical properties of integrated dynamic speckles produced by a rotating spheroidal object[J]. Journal of Optics, 1993, 24(3): 135-140.

[17] 李翰卿. 基于生物散斑的非接触式人体心率和脉搏波测量方法研究[D]. 厦门: 厦门大学, 2018.

Li HQ. Study on non-contact measuring method of human heart rate and pulse wave based on bio-speckle[D]. Xiamen: Xiamen University, 2018.

[18] Briers D. Time-varying laser speckle for measuring motion and flow[J]. Proceedings of SPIE, 2001, 4242: 25-39.

[19] Kazmi S M, Richards L M, Schrandt C J, et al. Expanding applications, accuracy, and interpretation of laser speckle contrast imaging of cerebral blood flow[J]. Journal Cerebral Blood Flow Metabolism, 2015, 35(7): 1076-1084.

[20] 陈俊波, 曾亚光, 袁治灵, 等. 基于动态散斑的光学相干层析成像技术[J]. 光学学报, 2018, 38(1): 0111001.

Chen J B, Zeng Y G, Yuan Z L, et al. Optical coherence tomography based on dynamic speckle[J]. Acta Optica Sinica, 2018, 38(1): 0111001.

[21] Fujii H, Asakura T, Nohira K, et al. Blood flow observed by time-varying laser speckle[J]. Optics Letters, 1985, 10(3): 104-106.

[22] Barton J, Stromski S. Flow measurement without phase information in optical coherence tomography images[J]. Optics Express, 2005, 13(14): 5234-5239.

[23] Humeau-Heurtier A, Guerreschi E, Abraham P, et al. Relevance of laser Doppler and laser speckle techniques for assessing vascular function: state of the art and future trends[J]. IEEE Transactions on Biomedical Engineering, 2013, 60(3): 659-666.

[24] Sawabe M, Kashimoto M. ADEOS-II precise orbit determinations with GPS and SLR[J]. Advances in Space Research, 1999, 23(4): 763-766.

[25] 钟声远, 李长桢, 陈念江, 等. 导航卫星激光后向反射器研究[J]. 激光与红外, 2011, 41(8): 834-839.

Zhong S Y, Li C Z, Chen N J, et al. Study of the laser retro-reflector on navigation satellites[J]. Laser & Infrared, 2011, 41(8): 834-839.

[26] 刘万里, 欧阳健飞, 曲兴华. 激光光束入射角度变化对角锥棱镜测量精度的影响[J]. 光学精密工程, 2009, 17(2): 286-291.

Liu W L, Ouyang J F, Qu X H. Effect of incident laser beam angle varying on cube corner retro-reflector measurement accuracty[J]. Optics and Precision Engineering, 2009, 17(2): 286-291.

3 典型运动状态下角反射器阵列的微动参数反演

王利国, 李亚清, 巩蕾, 王谦. 基于动态散斑时频特性的合作目标微动参数反演[J]. 光学学报, 2021, 41(2): 0228001. Liguo Wang, Yaqing Li, Lei Gong, Qian Wang. Inversion Algorithm for Micro-Motion Parameters of a Cooperative Target Based on Time-Frequency Feature of Dynamic Speckle[J]. Acta Optica Sinica, 2021, 41(2): 0228001.

基于动态散斑时频特性的合作目标微动参数反演 下载： 839次

1 引言

2 角反射器阵列激光动态散斑模型

图 1. 角反射器阵列结构及探测系统光路几何示意图。(a)角反射器阵列结构;(b)探测系统光路几何

Fig. 1. Structure of retroreflector array and optical geometry of detecting system. (a) Structure of retroreflector array; (b) optical geometry of detecting system

2.1 角反射器单元激光回波模型

2.2 角反射器阵列动态散斑强度模型

图 2. 自旋运动角反射器阵列散斑强度的时间序列

Fig. 2. Intensity time sequence of speckle field of a rotating retroreflector

2.3 角反射器阵列动态散斑时频信号模型

3 典型运动状态下角反射器阵列的微动参数反演

3.1 同轴旋转阵列

图 3. 同轴旋转阵列动态散斑的时频图

Fig. 3. Time-frequency feature of dynamic speckle of coaxial rotating array

图 4. 同轴旋转阵列动态散斑的频谱相关函数

Fig. 4. Spectral correlation function of dynamic speckle of coaxial rotating array

图 5. 同轴旋转阵列动态散斑的频谱累积强度分布

Fig. 5. Spectral integrated intensity of dynamic speckle of the coaxial rotating array

表 1. 各阶谱线的振幅值及其对应的视线角反演值

Table 1. Amplitude of each spectral lines and corresponding line-of-sight angles

3.2 异轴旋转阵列

图 6. 异轴旋转阵列动态散斑的时频图

Fig. 6. Time-frequency feature of dynamic speckle of non-coaxial rotating array

图 7. 异轴旋转阵列动态散斑的频谱相关函数

Fig. 7. Spectral correlation function of dynamic speckle of non-coaxial rotating array

图 8. 异轴旋转阵列动态散斑时频谱线的振幅-相位分布图

Fig. 8. Amplitude-phase distribution of time-frequency spectral lines of dynamic speckle of non-coaxial rotating array

表 2. 异轴旋转阵列时频谱线振幅的提取结果

Table 2. Extracted amplitude of time-frequency spectral lines of non-coaxial rotating array

3.3 复合运动阵列

图 9. 复合运动阵列动态散斑的时频图

Fig. 9. Time-frequency feature of dynamic speckle of the array with compound movement

图 10. 复合运动阵列动态散斑的频谱相关函数

Fig. 10. Spectral correlation function of dynamic speckle of the array with compound movement

图 11. 复合运动阵列散斑时频谱线的振幅-相位分布图

Fig. 11. Amplitude-phase distribution of time-frequency spectral lines of the array with compound movement

表 3. 复合运动阵列时频谱线的振幅以及对应的反演视线角

Table 3. Amplitude of time-frequency spectral lines of the array with compound movement and corresponding inverted line-of-sight angle

4 结论

Article Outline

相关论文

相关资讯

关于本站 Cookie 的使用提示

全站搜索

基于动态散斑时频特性的合作目标微动参数反演下载： 839次