一种基于多子带融合处理的电离层空变TEC测量方法

一种基于多子带融合处理的电离层空变tec测量方法
技术领域
1.本发明属于信号处理技术领域，具体涉及一种基于多子带融合处理的电离层空变tec测量方法。


背景技术：

2.电离层是地球大气的重要组成部分之一，其形成是日地空间环境各种因素综合作用的结果，涉及太阳活动、微粒辐射、地磁场分布及高层大气运动等因素。利用电离层的分布及变化信息可以用于太阳活动监测、地磁场渐变、气候天气预测等自然科学研究；同时，电离层作为一种介电常数随机起伏的色散媒质，会严重影响电磁波的传播，导致无线电信号发生畸变，进而影响包括卫星通信、卫星导航和星载雷达等在内的众多空间信息系统的信息获取质量。因此，实现电离层精细结构探测对空间科学研究、国民经济发展等具有重要的科学意义和应用价值。
3.国内外现有电离层探测方法主要包括：借助电离层测高仪探测底层或顶层电离层电子密度分布结构；利用地基非相干散射雷达获取电离层电子、离子密度和温度等特征参数；通过卫星搭载朗缪尔探针等各类等离子体粒子浓度探测仪器来实现电离层电子、离子等粒子密度测量。然而，上述方法存在探测范围受地基站布设限制、需搭载专业电离层探测载荷以及探测分辨率较低等劣势。针对这些问题，近年来相关学者提出采用长波长星载sar实现电离层高精度探测，利用雷达回波在电离层中传输所受到的电离层效应误差反演实现全球电离层高精度、高分辨率探测。
4.现有针对长波长星载sar实现电离层探测的研究(guo w.,chen j.,liu w.,et al.time-variant tec estimation with fully polarimetric geo-sar data[j].electronics letters,2017,53(24).)中，利用全极化sar系统信号受到的法拉第旋转效应影响，可以通过估计四极化通道包含的法拉第旋转角，并进一步利用该角度与电离层tec(单位为tecu，1tecu＝10
16
el/m-2
)线性关系实现电离层tec反演。然而，全极化sar信号会受到热噪声以及极化定标误差影响，现有研究仅利用载频点处的时域信号进行tec反演，通过采用较大的平滑窗进行平均处理降低各类误差影响，降低了空变tec的分辨率；特别对于快速空变tec分布结构，采用传统方法将严重影响tec的探测精度与分辨率，无法实现空变tec高精度、高分辨率探测。


技术实现要素：

[0005]
为了克服上述现有技术的缺点，本发明的目的在于提供了一种基于多子带融合处理的电离层空变tec测量方法，利用长波长全极化星载sar对电离层效应误差敏感性高的特点，通过频域多子带法拉第旋转角估计及融合处理，能够实现全球范围内电离层空变tec精确测量，其测量精度优于0.003tecu，分辨率优于70m。
[0006]
为了达到上述目的，本发明采取的技术方案为：
[0007]
一种基于多子带融合处理的电离层空变tec测量方法，包括以下步骤：
[0008]
步骤一：距离向快速傅里叶变换；
[0009]
针对全极化sar获取的四个极化通道图像数据s
hh
、s
hv
、s
vh
、s
vv
，沿每个方位采样时刻(按行)进行快速傅里叶变换(fft)，得到距离向频谱数据s
hh-f
、s
hv-f
、s
vh-f
、s
vv-f
；
[0010]
步骤二：计算距离向带宽采样点数；
[0011]
根据雷达系统参数，计算距离向带宽采样点数num_bw；
[0012]
步骤三：计算距离向频谱子带数据带宽采样点数；
[0013]
根据频谱子带分割数n，计算频谱子带数据所包含带宽采样点数num_sub；
[0014]
步骤四：沿距离向频谱分割处理，计算各极化通道频谱子带数据矩阵；
[0015]
结合步骤二及步骤三计算得到的距离向带宽采样点数num_bw以及频谱子带带宽采样点数num_sub，针对四个极化通道的距离向频谱数据s
hh-f
、s
hv-f
、s
vh-f
、s
vv-f
，分别利用频谱分割处理，计算四个极化通道各自对应的n个频谱子带数据矩阵；
[0016]
步骤五：距离向快速逆傅里叶变换；
[0017]
根据步骤四计算得到的四个极化通道各自对应的n个频谱子带数据矩阵，分别沿每个方位采样时刻(按行)快速逆傅里叶变换(ifft)处理，获取四个极化通道各自对应的n个距离向子带图像数据；
[0018]
步骤六：估计各子带图像所含法拉第旋转角；
[0019]
分别选取每个子带对应的四个极化通道子带图像数据，通过线性组合得到一组正交圆极化波信号z
12
和z
21
，通过该组正交圆极化波信号共轭相乘取幅角fr
temp
，对fr
temp
进行平滑处理后得到法拉第旋转角估计值fr，遍历第1至第n组子带图像得到所有子带图像对应的法拉第旋转角估计值fr1、fr2、
…
frn；
[0020]
步骤七：获取空变tec测量结果；
[0021]
根据雷达系统参数，计算各子带图像所对应中心频率f0，利用法拉第旋转角与tec线性对应关系分别计算各子带图像所含空变tec值，采用平均处理获取最终空变tec测量结果。
[0022]
本发明的优点在于：
[0023]
(1)本发明具有空变tec测量精度高的特点，由于采用基于法拉第旋转角的tec反演模型，相比较传统测量方法，具有更高的测量精度(优于0.003tecu)。
[0024]
(2)本发明具有空变tec测量分辨率高的特点，由于采用频谱分割处理方法，利用子带信息融合，相比较传统测量方法，可以实现更高分辨率的空变tec测量(优于70m)。
[0025]
(3)本发明具有全球范围覆盖探测的特点，由于采用星载sar体制，相比较传统测量方法，突破地面站布设区域限制，可以借助卫星平台宽幅测量优势实现全球范围内空变tec测量。
[0026]
(4)本发明具有测量设备简单便捷的特点，利用全极化星载sar系统，相比较传统测量方法，无需搭载额外的电离层探测专业设备(例如：朗缪尔探针、质谱仪和光谱仪等)，仅利用星载sar获取的全极化数据进行反演计算获取空变tec，极大地降低了电离层探测载荷复杂度以便于卫星搭载。
[0027]
(5)本发明具有电离层实时观测的特点，由于采用全极化星载sar数据反演，相比较传统测量方法，可以通过星上实时处理设备对所获取的全极化数据进行处理得到空变tec实时处理结果。
附图说明
[0028]
图1是本发明的流程图。
[0029]
图2是本发明实施例中在仿真全极化数据中加入的空变tec示意图。
[0030]
图3是本发明实施例中估计得到的空变tec示意图。
具体实施方式
[0031]
下面将结合附图对本发明作进一步详细说明。
[0032]
本发明提出一种基于多子带融合处理的电离层空变tec测量方法，利用全极化星载sar的四个极化通道数据分别表示为s
hh
、s
hv
、s
vh
、s
vv
，其均为大小na×
nr的二维复数矩阵，二维复数矩阵一维为方位向，有na个采样点，每个采样点对应的时刻代表雷达开启回波接收窗接收回波信号；二维复数矩阵另一维为距离向，有nr个采样点，表示在每个方位采样时刻，以采样率fs对带宽大小为bw的回波信号进行采样，一次连续距离向采样包含nr个采样点。
[0033]
如图1所示，一种基于多子带融合处理的电离层空变tec测量方法，包括以下步骤：
[0034]
步骤一：距离向快速傅里叶变换；
[0035]
针对全极化sar获取的四个极化通道图像数据s
hh
、s
hv
、s
vh
、s
vv
，沿每个方位采样时刻(按行)进行快速傅里叶变换(fft)，得到距离向频谱数据s
hh-f
、s
hv-f
、s
vh-f
、s
vv-f
，具体表示为：
[0036]
设全极化sar获取的四个极化通道数据s
hh
、s
hv
、s
vh
、s
vv
分别为：
[0037][0038][0039][0040][0041]
式中，表示hh通道图像中第(1,1)个像素的数值，表示hh通道图像中第(1,nr)个像素的数值，表示hh通道图像中第(na,1)个像素的数值，表示hh通道图像中第(na,nr)个像素的数值；表示hv通道图像中第(1,1)个像素的数值，表示hv通道图像中第(1,nr)个像素的数值，表示hv通道图像中第(na,1)个像素的数值，表示hv
通道图像中第(na,nr)个像素的数值；表示vh通道图像中第(1,1)个像素的数值，表示vh通道图像中第(1,nr)个像素的数值，表示vh通道图像中第(na,1)个像素的数值，表示vh通道图像中第(na,nr)个像素的数值；表示vv通道图像中第(1,1)个像素的数值，表示vv通道图像中第(1,nr)个像素的数值，表示vv通道图像中第(na,1)个像素的数值，表示vv通道图像中第(na,nr)个像素的数值；
[0042]
针对四个极化通道图像数据，沿每个方位采样时刻(按行)进行快速傅里叶变换(fft)，得到距离向频谱数据s
hh-f
、s
hv-f
、s
vh-f
、s
vv-f
分别为：
[0043][0044][0045][0046][0047]
式中，表示hh通道距离向频谱数据中第(1,1)个像素的数值，表示hh通道距离向频谱数据中第(1,nr)个像素的数值，表示hh通道距离向频谱数据中第(na,1)个像素的数值，表示hh通道距离向频谱数据中第(na,nr)个像素的数值；表示hv通道距离向频谱数据中第(1,1)个像素的数值，表示hv通道距离向频谱数据中第(1,nr)个像素的数值，表示hv通道距离向频谱数据中第(na,1)个像素的数值，表示hv通道距离向频谱数据中第(na,nr)个像素的数值；表示vh通道距离向频谱数据中第(1,1)个像素的数值，表示vh通道距离向频谱数据中第(1,nr)个像素的数值，表示vh通道距离向频谱数据中第(na,1)个像素的数值，表示vh通道距离向频谱数据中第(na,nr)个像素的数值；表示vv通道距离向频谱数据中第(1,1)个像素的数值，表示vv通道距离向频谱数据中第(1,nr)个像素的数值，表示vv通道距离向频谱数据中第(na,1)个像素的数值，表示vv通道距离向频谱数据中第(na,nr)个像素的数值；fft
(
·
)表示对二维矩阵的每行数据进行快速傅里叶变换；
[0048]
快速傅里叶变换(fft)和快速逆傅里叶逆变换(ifft)的具体内容，请参考2012年6月电子工业出版社、ian g cumming等著、洪文等译的《合成孔径雷达成像-算法与实现》一书，第18页至第19页；
[0049]
步骤二：计算距离向带宽采样点数；
[0050]
根据雷达系统参数，计算距离向带宽采样点数num_bw：
[0051]
num_bw＝floor(fs·
bw)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)
[0052]
式中，fs表示距离向信号采样率，bw表示距离向信号带宽，floor(x)表示取不大于x的最大偶数运算；
[0053]
步骤三：计算距离向频谱子带数据带宽采样点数；
[0054]
根据频谱子带分割数n，计算频谱子带数据所包含带宽采样点数num_sub：
[0055][0056]
式中，n表示距离向频谱子带分割数量，其取值为大于1的正整数；
[0057]
步骤四：沿距离向频谱分割处理，计算各极化通道频谱子带数据矩阵；
[0058]
结合步骤二及步骤三计算得到的距离向带宽采样点数num_bw以及频谱子带带宽采样点数num_sub，针对四个极化通道的距离向频谱数据s
hh-f
、s
hv-f
、s
vh-f
、s
vv-f
，分别利用频谱分割处理，计算四个极化通道各自对应的n个频谱子带数据矩阵，具体为：
[0059]
4.1)针对hh极化通道数据，将原始数据带宽频谱等分为n份，构造n个频谱子带数据矩阵，其大小均为na×
nr，其中第n个频谱子带数据矩阵为(其中1≤n≤n)：
[0060][0061]
式中，表示hh极化通道第n个频谱子带数据矩阵中第(1，nr/2-num_bw/2+1+(n-1)
·
num_sub)个像素的数值，其值等于hh极化通道距离向频谱数据矩阵s
hh-f
中第(1，nr/2-num_bw/2+1+(n-1)
·
num_sub)个像素的数值；表示hh极化通道第n个频谱子带数据矩阵中第(1，nr/2-num_bw/2+n
·
num_sub)个像素的数值，其值等于hh极化通道距离向频谱数据矩阵s
hh-f
中第(1，nr/2-num_bw/2+n
·
num_sub)个像素的数值；表示hh极化通道第n个频谱子带数据矩阵中第(na，nr/2-num_bw/2+1+(n-1)
·
num_sub)个像素的数值，其值等于hh极化通道距离向频谱数据矩阵s
hh-f
中第(na，nr/2-num_bw/2+1+(n-1)
·
num_sub)个像素的数值；表示hh极化通道第n个频谱子带数据矩阵中第(na，nr/2-num_bw/2+n
·
num_sub)个像素的数值，其值等于hh极化通道距离向频谱数据矩阵s
hh-f
中第(na，nr/2-num_bw/2+n
·
num_sub)个像素的数值；
[0062]
4.2)针对hv极化通道数据，将原始数据带宽频谱等分为n份，构造n个频谱子带数据矩阵，其大小均为na×
nr，其中第n个频谱子带数据矩阵为(其中1≤n≤n)：
[0063][0064]
式中，表示hv极化通道第n个频谱子带数据矩阵中第(1，nr/2-num_bw/2+1+(n-1)
·
num_sub)个像素的数值，其值等于hv极化通道距离向频谱数据矩阵s
hv-f
中第(1，nr/2-num_bw/2+1+(n-1)
·
num_sub)个像素的数值；表示hv极化通道第n个频谱子带数据矩阵中第(1，nr/2-num_bw/2+n
·
num_sub)个像素的数值，其值等于hv极化通道距离向频谱数据矩阵s
hv-f
中第(1，nr/2-num_bw/2+n
·
num_sub)个像素的数值；表示hv极化通道第n个频谱子带数据矩阵中第(na，nr/2-num_bw/2+1+(n-1)
·
num_sub)个像素的数值，其值等于hv极化通道距离向频谱数据矩阵s
hv-f
中第(na，nr/2-num_bw/2+1+(n-1)
·
num_sub)个像素的数值；表示hv极化通道第n个频谱子带数据矩阵中第(na，nr/2-num_bw/2+n
·
num_sub)个像素的数值，其值等于hv极化通道距离向频谱数据矩阵s
hv-f
中第(na，nr/2-num_bw/2+n
·
num_sub)个像素的数值；
[0065]
4.3)针对vh极化通道数据，将原始数据带宽频谱等分为n份，构造n个频谱子带数据矩阵，其大小均为na×
nr，其中第n个频谱子带数据矩阵为(其中1≤n≤n)：
[0066][0067]
式中，表示vh极化通道第n个频谱子带数据矩阵中第(1，nr/2-num_bw/2+1+(n-1)
·
num_sub)个像素的数值，其值等于vh极化通道距离向频谱数据矩阵s
vh-f
中第(1，nr/2-num_bw/2+1+n-1)
·
num_sub)个像素的数值；表示vh极化通道第n个频谱子带数据矩阵中第(1，nr/2-num_bw/2+n
·
num_sub)个像素的数值，其值等于vh极化通道距离向频谱数据矩阵s
vh-f
中第(1，nr/2-num_bw/2+n
·
num_sub)个像素的数值；表示vh极化通道第n个频谱子带数据矩阵中第(na，nr/2-num_bw/2+1+(n-1)
·
num_sub)个像素的数值，其值等于vh极化通道距离向频谱数据矩阵s
vh-f
中第(na，nr/2-num_bw/2+1+(n-1)
·
num_sub)个像素的数值；表示vh极化通道第n个频谱子带数据矩阵中第(na，nr/2-num_bw/2+n
·
num_sub)个像素的数值，其值等于vh极化通道距离向频谱数据矩阵s
vh-f
中第(na，nr/2-num_bw/2+n
·
num_sub)个像素的数值；
[0068]
4.4)针对vv极化通道数据，将原始数据带宽频谱等分为n份，构造n个频谱子带数据矩阵，其大小均为na×
nr，其中第n个频谱子带数据矩阵为(其中1≤n≤n)：
[0069][0070]
式中，表示vv极化通道第n个频谱子带数据矩阵中第(1，nr/2-num_bw/2+1+(n-1)
·
num_sub)个像素的数值，其值等于vv极化通道距离向频谱数据矩阵s
vv-f
中第(1，nr/2-num_bw/2+1+n-1)
·
num_sub)个像素的数值；表示vv极化通道第n个频谱子带数据矩阵中第(1，nr/2-num_bw/2+n
·
num_sub)个像素的数值，其值等于vv极化通道距离向频谱数据矩阵s
vv-f
中第(1，nr/2-num_bw/2+n
·
num_sub)个像素的数值；表示vv极化通道第n个频谱子带数据矩阵中第(na，nr/2-num_bw/2+1+(n-1)
·
num_sub)个像素的数值，其值等于vv极化通道距离向频谱数据矩阵s
vv-f
中第(na，nr/2-num_bw/2+1+(n-1)
·
num_sub)个像素的数值；表示vv极化通道第n个频谱子带数据矩阵中第(na，nr/2-num_bw/2+n
·
num_sub)个像素的数值，其值等于vv极化通道距离向频谱数据矩阵s
vv-f
中第(na，nr/2-num_bw/2+n
·
num_sub)个像素的数值；
[0071]
步骤五：距离向快速逆傅里叶变换；
[0072]
根据步骤四计算得到的四个极化通道各自对应的n个频谱子带数据矩阵，分别沿每个方位采样时刻(按行)快速傅里叶逆变换(ifft)处理，获取四个极化通道各自对应的n个距离向子带图像数据，具体为：
[0073]
5.1)将步骤四得到的hh极化通道的n个频谱子带数据矩阵，沿每个方位采样时刻(按行)进行快速逆傅里叶变换(ifft)，得到距离向子带图像数据矩阵，其中第n个子带图像数据矩阵为：
[0074][0075]
式中，表示hh通道第n个子带图像数据矩阵中第(1,1)个像素的数值，表示hh通道第n个子带图像数据矩阵中第(1,nr)个像素的数值，表示hh通道第n个子带图像数据矩阵中第(na,1)个像素的数值，表示hh通道第n个子带图像数据矩阵中第(na,nr)个像素的数值；
[0076]
5.2)将步骤四得到的hv极化通道的n个频谱子带数据矩阵，沿每个方位采样时刻(按行)进行快速逆傅里叶变换(ifft)，得到距离向子带图像数据矩阵，其中第n个子带图像数据矩阵为：
[0077]
[0078]
式中，表示hv通道第n个子带图像数据矩阵中第(1,1)个像素的数值，表示hv通道第n个子带图像数据矩阵中第(1,n)个像素的数值，表示hv通道第n个子带图像数据矩阵中第(na,1)个像素的数值，表示hv通道第n个子带图像数据矩阵中第(na,nr)个像素的数值；
[0079]
5.3)将步骤四得到的vh极化通道的n个频谱子带数据矩阵，沿每个方位采样时刻(按行)进行快速逆傅里叶变换(ifft)，得到距离向子带图像数据矩阵，其中第n个子带图像数据矩阵为：
[0080][0081]
式中，表示vh通道第n个子带图像数据矩阵中第(1,1)个像素的数值，表示vh通道第n个子带图像数据矩阵中第(1,nr)个像素的数值，表示vh通道第n个子带图像数据矩阵中第(na,1)个像素的数值，表示vh通道第n个子带图像数据矩阵中第(na,nr)个像素的数值；
[0082]
5.4)将步骤四得到的vv极化通道的n个频谱子带数据矩阵，沿每个方位采样时刻(按行)进行快速逆傅里叶变换(ifft)，得到距离向子带图像数据矩阵，其中第n个子带图像数据矩阵为：
[0083][0084]
式中，表示vv通道第n个子带图像数据矩阵中第(1,1)个像素的数值，表示vv通道第n个子带图像数据矩阵中第(1,nr)个像素的数值，表示vv通道第n个子带图像数据矩阵中第(na,1)个像素的数值，表示vv通道第n个子带图像数据矩阵中第(na,nr)个像素的数值；
[0085]
步骤六：估计各子带图像所含法拉第旋转角；
[0086]
分别选取每个子带对应的四个极化通道子带图像数据，通过线性组合得到一组正交圆极化波信号z
12
和z
21
，通过该组正交圆极化波信号共轭相乘取幅角fr
temp
，对fr
temp
进行平滑处理后得到法拉第旋转角估计值fr，遍历第1至第n组子带图像得到所有子带图像对应的法拉第旋转角估计值fr1、fr2、
…
、frn，具体为：
[0087]
6.1)选取第1个子带对应的四个极化通道子带图像数据进行线性组合处理，计算一组正交圆极化波信号z
12_1
和z
21_1
：
[0088]z12_1
＝s
vh_sub_1-s
hv_sub_1
+j
·
(s
hh_sub_1
+s
vv_sub_1
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(19)
[0089]z21_1
＝s
hv_sub_1-s
vh_sub_1
+j
·
(s
hh_sub_1
+s
vv_sub_1
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(20)
[0090]
式中，z
12_1
和z
21_1
分别代表第1个子带对应计算得到的一组正交圆极化波信号矩
阵，其均为大小为na×
nr的二维复数矩阵，j为虚数单位，即
[0091]
6.2)对步骤6.1)中一组正交圆极化波信号矩阵共轭相乘取幅角fr
temp_1
：
[0092][0093]
式中，fr
temp_1
表示利用正交圆极化波信号共轭相乘计算得到的幅角值，其为大小为na×
nr的二维实数矩阵，
*
表示共轭，angle(x)为取复数x的幅角运算；
[0094]
6.3)对上述计算得到的幅角值fr
temp_1
(二维实数矩阵)，利用矩形平滑窗对fr
temp_1
进行平滑处理得到法拉第旋转角估计值fr1，其结果为大小na×
nr的二维实数矩阵：
[0095]
fr1＝《fr
temp_1
》
ꢀꢀ
(22)
[0096]
式中，《x》表示对x采用窗口大小为(2m+1，2m+1)的窗函数进行平滑处理，即以x所在坐标(p，q)为中心(其中1≤p≤na，1≤q≤nr)，求周围第p-m行到p+m行，第q-m列到q+m列内的所有分量均值作为x值，需要注意的是，对于靠近矩阵边缘的元素，会存在平滑处理窗口超出矩阵边界的情况，即p-m≤0，p+m≥na，q-m≤0，q+m≥nr，此时，该元素不做平滑处理直接以计算出x值作为该位置处结果；
[0097]
6.4)遍历第1至第n组子带图像，利用上述步骤6.1)至6.3)子步骤分别计算每组子带图像对应的法拉第旋转角fr1、fr2、
…
、frn；
[0098]
步骤七：获取空变tec测量结果；
[0099]
根据雷达系统参数，计算各子带图像所对应中心频率f0，利用法拉第旋转角与tec线性对应关系分别计算各子带图像所含空变tec值，采用平均处理获取最终空变tec测量结果，具体为：
[0100]
7.1)根据sar系统工作载频fc、带宽bw以及频谱子带分割数n，分别计算第1至第n组子带图像各自所对应的中心频率f0(1)、f0(2)、
…
、f0(n)，其中第n个子带图像对应的中心频率f0(n)为：
[0101][0102]
式中，fc表示雷达系统工作载频，1≤n≤n；
[0103]
7.2)根据步骤六得到的各子带图像对应法拉第旋转角值fr，利用tec与法拉第旋转角间的线性转换关系，获得各子带图像所含空变tec值，其中第n个子带图像对应的空变tec为：
[0104][0105]
式中，frn表示计算得到的第n个子带图像对应的法拉第旋转角，k
ω
为常数且k
ω
＝2.365
×
104a
×
m2/kg，b表示平行于雷达波束方向的地球磁场强度，单位为wb/m2；
[0106]
7.3)根据n个子带图像估计得到的空变tec值，利用平均处理获取最终空变tec为：
[0107][0108]
经过以上七个步骤，利用全极化星载sar信号完成了对空变tec的精确测量。
[0109]
下面将结合实施例对本发明作进一步的详细说明。
[0110]
由于现有全极化星载sar系统所获取的数据均已受到电离层效应影响，本身包含法拉第旋转角，因此本实施例中采用仿真得到全极化星载sar数据，并在四个极化通道信号中加入电离层空变tec(如图2所示)所引起的法拉第旋转角，处理过程中所涉及的参数如表1所示：
[0111]
表1实施例参数
[0112][0113][0114]
本实施例具体包含以下步骤：
[0115]
步骤一：距离向快速傅里叶变换；
[0116]
针对全极化sar获取的四个极化通道图像数据s
hh
、s
hv
、s
vh
、s
vv
，均为大小为2048
×
2048的二维复数矩阵，按式(5)(6)(7)(8)沿每个方位采样时刻(按行)进行快速傅里叶变换(fft)，得到距离向频谱数据s
hh-f
、s
hv-f
、s
vh-f
、s
vv-f
，均为大小为2048
×
2048的二维复数矩阵；
[0117]
步骤二：计算距离向带宽采样点数；
[0118]
根据雷达系统参数，按式(9)计算距离向带宽采样点数1486；
[0119]
步骤三：计算距离向频谱子带数据带宽采样点数；
[0120]
根据频谱子带分割数8，按式(10)计算频谱子带数据所包含带宽采样点数186；
[0121]
步骤四：沿距离向频谱分割处理，计算各极化通道频谱子带数据矩阵；
[0122]
结合步骤二及步骤三计算得到的距离向带宽采样点数1486以及频谱子带带宽采样点数186，针对四个极化通道的距离向频谱数据s
hh-f
、s
hv-f
、s
vh-f
、s
vv-f
，分别利用频谱分割处理，计算四个极化通道各自对应的8个频谱子带数据矩阵，具体为：
[0123]
4.1)针对hh极化通道数据，将原始数据带宽频谱等分为8份，按式(11)构造8个频谱子带数据矩阵，其均为大小为2048
×
2048的二维复数矩阵；
[0124]
4.2)针对hv极化通道数据，将原始数据带宽频谱等分为8份，按式(12)构造8个频谱子带数据矩阵，其均为大小为2048
×
2048的二维复数矩阵；
[0125]
4.3)针对vh极化通道数据，将原始数据带宽频谱等分为8份，按式(13)构造8个频谱子带数据矩阵，其均为大小为2048
×
2048的二维复数矩阵；
[0126]
4.4)针对vv极化通道数据，将原始数据带宽频谱等分为8份，按式(14)构造8个频谱子带数据矩阵，其均为大小为2048
×
2048的二维复数矩阵；
[0127]
步骤五：距离向快速逆傅里叶变换；
[0128]
根据步骤四计算得到的四个极化通道各自对应的8个频谱子带数据矩阵，分别沿每个方位采样时刻(按行)快速傅里叶逆变换(ifft)处理，获取四个极化通道各自对应的8个距离向子带图像数据，具体为：
[0129]
5.1)将步骤四得到的hh极化通道的8个频谱子带数据矩阵，按式(15)沿每个方位采样时刻(按行)进行快速逆傅里叶变换(ifft)，得到8个距离向子带图像数据矩阵，其均为大小为2048
×
2048的二维复数矩阵；
[0130]
5.2)将步骤四得到的hv极化通道的8个频谱子带数据矩阵，按式(16)沿每个方位采样时刻(按行)进行快速逆傅里叶变换(ifft)，得到8个距离向子带图像数据矩阵，其均为大小为2048
×
2048的二维复数矩阵；
[0131]
5.3)将步骤四得到的vh极化通道的8个频谱子带数据矩阵，按式(17)沿每个方位采样时刻(按行)进行快速逆傅里叶变换(ifft)，得到8个距离向子带图像数据矩阵，其均为大小为2048
×
2048的二维复数矩阵；
[0132]
5.4)将步骤四得到的vv极化通道的8个频谱子带数据矩阵，按式(18)沿每个方位采样时刻(按行)进行快速逆傅里叶变换(ifft)，得到8个距离向子带图像数据矩阵，其均为大小为2048
×
2048的二维复数矩阵；
[0133]
步骤六：估计各子带图像所含法拉第旋转角；
[0134]
分别选取每个子带对应的四个极化通道子带图像数据，通过线性组合得到一组正交圆极化波信号z
12
和z
21
，通过该组正交圆极化波信号共轭相乘取幅角fr
temp
，对fr
temp
进行平滑处理后得到法拉第旋转角估计值fr，遍历第1至第8组子带图像得到所有子带图像对应的法拉第旋转角估计值fr1、fr2、
…
、fr8，具体为：
[0135]
6.1)选取第1个子带对应的四个极化通道子带图像数据进行线性组合处理，按式(19)、式(20)计算一组正交圆极化波信号z
12_1
和z
21_1
，其均为大小为2048
×
2048的二维复数矩阵；
[0136]
6.2)按式(21)对该组正交圆极化波信号矩阵共轭相乘取幅角fr
temp_1
，其为大小为2048
×
2048的二维实数矩阵；
[0137]
6.3)对上述计算得到的幅角值fr
temp_1
，按式(22)采用大小为(2m+1＝17，2m+1＝17)矩形平滑窗对fr
temp_1
进行平滑处理得到法拉第旋转角估计值fr1，即以fr
temp_1
矩阵中第(p，q)项为中心(其中1≤p≤2048，1≤q≤2048)，求周围第p-8行到p+8行，第q-8列到q+8列内的所有分量均值作为第(p，q)项值，需要注意的是，对于靠近矩阵边缘的元素，即第1至8行、第2041至2048行、第1至8列以及第2041至2048列中所有元素，会存在平滑处理窗口超出矩阵边界的情况，因此这些项不做平滑处理，直接以平滑前数值作为该项结果；由此，计算得到法拉第旋转角估计值fr1，其为大小为2048
×
2048的二维实数矩阵；
[0138]
6.4)遍历第1至第8组子带图像，利用上述步骤6.1)至6.3)子步骤分别计算每组子带图像对应的法拉第旋转角fr1、fr2、
…
、frn，其为大小为2048
×
2048的二维实数矩阵；
[0139]
步骤七：获取空变tec测量结果；
[0140]
根据雷达系统参数，计算各子带图像所对应中心频率f0，利用法拉第旋转角与tec
线性对应关系分别计算各子带图像所含空变tec值，采用平均处理获取最终空变tec测量结果，具体为：
[0141]
7.1)根据sar系统工作载频1236.5mhz、带宽38mhz以及频谱子带分割数8，按式(23)分别计算第1至第8组子带图像各自所对应的中心频率f0(1)、f0(2)、
…
、f0(n)；
[0142]
7.2)根据步骤六得到的各子带图像对应法拉第旋转角值fr，利用tec与法拉第旋转角间的线性转换关系，按式(24)计算获得各子带图像所含空变tec值；
[0143]
7.3)根据8个子带图像估计得到的空变tec值，按式(25)利用平均处理计算获取最终空变tec，如图3所示。
[0144]
对比图2所示仿真中加入的原始空变tec与根据本实施例计算得到的图3所示估计得到的空变tec，可以看出两者极为接近，计算两者相关系数为0.9924；进一步统计采用本发明方法计算得到的空变tec估计误差标准差为0.0022tecu，即本发明方法可实现测量精度优于0.003tecu的空变tec高精度测量；另外，结合本实施例中所采用的雷达系统带宽以及平滑处理窗口参数，计算得到空变tec测量分辨率为67.11m，即本发明方法可实现分辨率优于70m的空变tec高分辨率测量。
[0145]
综上，上述实施例说明本发明一种基于多子带融合处理的电离层空变tec测量方法可以实现空变tec高精度、高分辨率测量。