一种基于全极化高轨SAR的电离层时变TEC测量方法与流程
本发明涉及一种基于全极化高轨合成孔径雷达(sar)的电离层时变总电子量(tec)测量方法，属于信号处理技术领域。
背景技术：
近年来，随着科技发展与社会进步，卫星通信、卫星导航和天基雷达系统已广泛应用于军事和民用的各个方面，并成为人类生活必不可少的工具，这使得包括电离层在内的空间环境监测和技术保障愈显重要和迫切。同时，biomass、nisar等卫星计划将对全球包括电离层在内的大气环境进行探测，该探测数据将为地球空间气象科学以及空间信息系统气象保障等方面提供重要帮助，具有极为重大科学意义和应用价值。
当前国内外针对电离层探测仍主要依靠传统电离层探测手段，即地面站垂测、基于gps信号掩星反演等。但是，由于存在探测区域受地基站布设网络限制，需要搭载专门的电离层探测载荷等问题，导致无法实现全球电离层高精度测量。针对该问题，近年来国外学者提出了利用低频段sar信号在电离层传播过程中会引入电离层效应误差，通过误差估计进而对电离层特征参数进行反演，实现电离层高精度探测。现有研究中，利用全极化sar系统，可以通过估计电离层所引入的法拉第旋转角，根据法拉第旋转角与电离层tec的映射关系，实现电离层空变tec(单位为tecu)反演。然而，由于现有研究中所涉及的均为中低轨星载全极化sar系统，卫星飞行速度较快且合成孔径时间较短，仅能针对不同照射地区的电离层tec分布(即空变tec)进行探测，而无法实现固定区域的时变tec精确探测。
同时，目前针对高轨sar的研究仅仅停留在电离层所引起的成像误差补偿上，多数均采用相位梯度自聚焦的处理方法估计电离层时变tec所引起的误差相位。但由于该方法处理过程中不采用固定的误差模型，仅基于像素间相位梯度变化估计误差相位，从而完全丢失掉电离层时变tec所引起的绝对相位误差，因此无法实现时变tec精确探测。
技术实现要素：
本发明的目的是为了解决利用星载sar系统实现电离层时变tec测量的瓶颈技术问题，提出一种基于全极化高轨sar的电离层时变tec测量方法，充分利用高轨sar飞行速度较慢且具有较长合成孔径时间，对于同一地区可实现较长时间持续观测，以及采用全极化信号可精确估计法拉第旋转角的特点，实现基于全极化高轨sar的电离层时变tec测量新体制。该方法能够实现全球实时电离层时变tec精确测量，获取毫秒级时变tec测量结果，扩展了sar在遥感科学与地球空间气象科学中的应用。
一种基于全极化高轨sar的电离层时变tec测量方法，包括以下几个步骤：
步骤一：计算方位向合成孔径点数；
根据雷达系统参数，计算方位向合成孔径点数num_a；
步骤二：方位向数据补零；
根据步骤一所计算得到的方位向合成孔径点数num_a，分别在每个极化通道图像数据每一列的头部和尾部补零，补零数目均为num_a/2，得到新的四个极化通道图像数据分别为mhh-ins、mhv-ins、mvh-ins、mvv-ins；
步骤三：方位向傅里叶变换；
将步骤二得到的补零后的四个极化通道图像数据mhh-ins、mhv-ins、mvh-ins、mvv-ins沿每个距离门(按列)进行快速傅里叶变换(fft)，得到方位向频谱数据mhh-fft、mhv-fft、mvh-fft、mvv-fft；
步骤四：方位向信号解压缩，获取方位向回波信号；
根据雷达系统参数，计算每个距离门(每列)对应的参考斜距rref，以此为参考计算出方位向信号解压缩因子φdecom，利用该解压缩因子乘以步骤三所得的方位向频谱数据mhh-fft、mhv-fft、mvh-fft、mvv-fft得到等效方位向回波频谱信号mhh-decomfft、mhv-decom-fft、mvh-decom-fft、mvv-decom-fft；
步骤五：方位向逆傅里叶变换；
将步骤四得到的等效方位向回波频谱信号mhh-decom-fft、mhv-decom-fft、mvh-decom-fft、mvv-decom-fft沿每个距离门(按列)进行快速逆傅里叶变换(ifft)，得到方位向回波信号mhh-echo、mhv-echo、mvh-echo、mvv-echo；
步骤六：估计法拉第旋转角；
根据步骤五得到的四个极化通道方位向回波信号mhh-echo、mhv-echo、mvh-echo、mvv-echo，通过线性组合得到一组正交圆极化波信号z12和z21，通过该组正交圆极化波信号共轭相乘取幅角frtemp，对frtemp进行平滑处理后得到法拉第旋转角估计值fr；
步骤七：获取合成孔径时间内时变tec；
根据步骤六得到的法拉第旋转角估计值fr，利用时变tec与法拉第旋转角间的线性转换关系，获得合成孔径时间内每个方位时刻的tec，即合成孔径时间内的时变tec。
本发明的优点在于：
(1)本发明提出一种基于全极化高轨sar的电离层时变tec测量方法，该方法具有时变tec测量精度高的特点。由于采用法拉第旋转角模型与tec反演模型，相对于传统测量方法，具有更高的测量精度。
(2)本发明提出一种基于全极化高轨sar的电离层时变tec测量方法，该方法具有可持续测量固定区域tec时变特性的特点。由于采用高轨sar体制，相比较传统测量方法，利用雷达合成孔径时间长的特性，可在分钟量级的合成孔径时间段内对固定区域进行时变tec持续测量。
(3)本发明提出一种基于全极化高轨sar的电离层时变tec测量方法，该方法具有可测量时变tec短时间间隔时变特性的特点。由于采用高轨sar体制，相比较传统测量方法，利用雷达脉冲重复间隔短的特性，可实现变化时间间隔为毫秒量级的时变tec测量。
(4)本发明提出一种基于全极化高轨sar的电离层时变tec测量方法，该方法具有实时观测的特点。由于采用全极化高轨sar体制，相比较传统测量体制，可通过星上实时处理sar系统获取的全极化数据得到时变tec实时测量结果。
(5)本发明提出一种基于全极化高轨sar的电离层时变tec测量方法，该方法具有全球覆盖观测的特点。由于采用全极化高轨sar体制，相比较传统测量体制，不受地面站布设网络的限制，可利用高轨sar全球覆盖观测的特点实现全球时变tec测量。
(6)本发明提出一种基于全极化高轨sar的电离层时变tec测量方法，该方法具有观测设备简单的特点。由于采用全极化高轨sar体制，相比较传统测量体制，无需搭载专门用于电离层探测的探测设备(如：无线电等离子成像仪、电离层顶层探测器等)，可直接利用全极化高轨sar获取的全极化数据进行时变tec反演，极大地节省探测成本且便于卫星搭载。
附图说明
图1是本发明提出的一种基于全极化高轨sar的电离层时变tec测量方法的方法流程图；
图2是本发明实施例中在仿真全极化数据中加入的合成孔径时间内的时变tec；
图3是本发明实施例中采用本发明提出方法估计得到的合成孔径时间内的时变tec。
具体实施方式
下面将结合附图对本发明作进一步详细说明。
本发明的一种基于全极化高轨sar的电离层时变tec测量方法，该方法基于全极化高轨sar的四个极化通道数据(二维复数矩阵)分别表示为mhh、mhv、mvh、mvv，设其大小均为na×nr，一维是方位向，有na个采样点，prf表示雷达的脉冲重复频率，1/prf称为雷达的脉冲重复周期，不同的方位向采样点对应不同的方位时刻，相邻两个方位向采样点相差时间间隔为1/prf；另一维是距离向，有nr个采样点，表示每个方位时刻开启一次回波接收窗，对回波信号进行采样，一次连续距离向采样有nr个采样点，采样率为fs。
本发明是一种基于全极化高轨sar的电离层时变tec测量方法，方法流程图如图1所示，具体包括以下几个步骤：
步骤一：计算方位向合成孔径点数；
根据雷达系统参数，计算方位向合成孔径点数num_a：
式中，λ表示雷达系统工作波长，ro表示场景中心点的参考斜距，la表示方位向天线长度，vref表示以场景中心点为参考时的等效速度，prf表示雷达脉冲重复频率。floor(x)表示取不大于x的最大偶数运算。
步骤二：方位向数据补零；
根据步骤一计算得到的方位向合成孔径点数num_a，分别在每个极化通道图像数据每一列的头部和尾部补零，补零数目均为num_a/2，得到新的四个极化通道图像数据分别为mhh-ins、mhv-ins、mvh-ins、mvv-ins，具体为：
设原始四个极化通道数据mhh、mhv、mvh、mvv分别为：
式中，
结合步骤一计算得到的方位向合成孔径点数num_a，分别在四个极化通道数据每一列的头部和尾部分别补零num_a/2，得到补零后新的四个极化通道图像数据mhh-ins、mhv-ins、mvh-ins、mvv-ins，其均扩大为大小na-ins×nr-ins的二维复数矩阵，且na-ins＝na+num_a，nr-ins＝nr，四个通道的信号分别为：
式中，
步骤三：方位向傅里叶变换；
将步骤二得到的补零后的四个极化通道图像数据mhh-ins、mhv-ins、mvh-ins、mvv-ins沿每个距离门(按列)进行快速傅里叶变换(fft)，得到方位向频谱数据mhh-fft、mhv-fft、mvh-fft、mvv-fft，具体表示为：
式中，
在本发明中，快速傅里叶变换(fft)和快速逆傅里叶逆变换(ifft)的具体内容，请参考2012年6月电子工业出版社、iangcumming等著、洪文等译的《合成孔径雷达成像-算法与实现》一书，第18页至第19页。
步骤四：方位向信号解压缩，获取等效方位向回波频谱信号；
根据雷达系统参数，计算每个距离门(每列)对应的参考斜距rref，以此为参考计算出方位向信号解压缩因子φdecom，利用该解压缩因子乘以步骤三所得的方位向频谱数据mhh-fft、mhv-fft、mvh-fft、mvv-fft得到等效方位向回波频谱信号mhh-decomfft、mhv-decom-fft、mvh-decom-fft、mvv-decom-fft，具体为：
(1)根据雷达系统参数，计算每个距离门(每列)对应的参考斜距rref：
式中，c表示光速。需要注意的是，式(14)为每个距离门(每列)对应的参考斜距rref二维矩阵表示方式，式中处于同一距离门(每列)不同方位向(每行)的参考斜距数值是相同的，即每个距离门(每列)采用相同的参考斜距。
(2)雷达的脉冲重复频率为prf，计算方位向频谱数据每个方位向(每行)的方位频率fa为：
(3)根据上述计算所得每个距离门(每列)对应的参考斜距rref以及每个方位向(每行)的方位频率fa，可计算得到方位向信号解压缩因子φdecom:
式中，
(4)将计算得到的方位向信号解压缩因子φdecom乘以方位向频谱数据mhh-fft、mhv-fft、mvh-fft、mvv-fft得到等效方位向回波频谱信号mhh-decom-fft、mhv-decom-fft、mvh-decom-fft、mvv-decom-fft：
式中，j表示虚数单位，即
步骤五：方位向逆傅里叶变换；
将步骤四得到的等效方位向回波频谱信号mhh-decom-fft、mhv-decom-fft、mvh-decom-fft、mvv-decom-fft沿每个距离门(按列)进行快速逆傅里叶变换(ifft)，得到方位向回波信号mhh-echo、mhv-echo、mvh-echo、mvv-echo，具体表示为：
式中，
步骤六：估计法拉第旋转角；
根据步骤五得到的四个极化通道方位向回波信号mhh-echo、mhv-echo、mvh-echo、mvv-echo，通过线性组合得到一组正交圆极化波信号z12和z21，通过该组正交圆极化波信号共轭相乘取幅角frtemp，对frtemp进行平滑处理后得到法拉第旋转角估计值fr。具体为：
(1)对四个极化通道方位向回波信号进行线性组合处理，计算一组正交圆极化波信号z12和z21：
z12＝mvh-echo-mhv-echo+j×(mhh-echo+mvv-echo)(25)
z21＝mhv-echo-mvh-echo+j×(mhh-echo+mvv-echo)(26)
式中，z12和z21分别代表计算得到的一组正交圆极化波信号，其均为大小为na-ins×nr-ins的二维复数矩阵。
(2)对该组正交圆极化波信号共轭相乘取幅角frtemp：
式中，frtemp表示利用正交圆极化波信号共轭相乘计算得到的幅角值，其为大小为na-ins×nr-ins的二维实数矩阵，*表示共轭，angle(x)为取复数x的幅角运算。
(3)对上述计算得到的幅角值frtemp(二维实数矩阵)，在仅考虑时变tec的条件下首先将frtemp按列叠加求平均值，再对得到的结果(一维实数数组)进行平滑处理得到法拉第旋转角估计值fr，其为长度为na-ins的一维实数数组：
式中，frtemp(:,i)表示幅角值矩阵frtemp第i列所有元素，<x>表示对x采用窗口大小为(2m+1)的窗函数进行平滑处理，即以x所在坐标p为中心，求前后第p-m到p+m个内的所有元素均值作为x值，需要注意的是，对于靠近数组边缘的元素，会存在平滑处理窗口超出数组边界的情况，即p-m≤0，p+m≥na-ins，此时，该元素不做平滑处理直接以计算出x值作为该处结果。
步骤七：获取合成孔径时间内的时变tec；
根据步骤六得到的估计法拉第旋转角值fr，利用时变tec与法拉第旋转角间的线性转换关系，获得合成孔径时间内每个方位时刻的tec，即合成孔径时间内的时变tec：
式中，fr(n)表示估计法拉第旋转角值fr中第n个元素值，kω为常数且kω＝2.365×104a×m2/kg，b表示平行于雷达波束方向的地球磁场强度，单位为wb/m2。
经过以上七个步骤，利用全极化高轨sar信号完成了对合成孔径时间内的时变tec精确测量。
实施例
下面将结合附图和实施例对本发明作进一步的详细说明。
本实施例提供一种基于全极化高轨sar的电离层时变tec测量方法。由于现有星载sar系统没有直接获取的全极化高轨sar数据，因此实施例中采用仿真得到全极化高轨sar数据，并在信号方位向合成孔径时间内加入电离层时变tec(如图2所示)所引起的法拉第旋转角，处理过程中所涉及的参数如表1所示：
表1实施例参数
本实施例具体包含以下步骤：
步骤一：计算方位向合成孔径点数；
根据雷达系统参数，按式(1)计算方位向合成孔径点数为7768。
步骤二：方位向数据补零；
根据步骤一所计算得到的方位向合成孔径点数7768，按式(6)(7)(8)(9)分别在每个极化通道图像数据mhh、mhv、mvh、mvv每一列的头部和尾部补零，补零数目均为3884，得到新的四个极化通道图像数据分别为mhh-ins、mhv-ins、mvh-ins、mvv-ins。原始每个极化通道图像数据mhh、mhv、mvh、mvv均为大小为1024×1024的二维复数矩阵，经过方位向数据补零处理后，得到新的四个极化通道图像数据mhh-ins、mhv-ins、mvh-ins、mvv-ins均为大小为8792×1024的二维复数矩阵。
步骤三：方位向傅里叶变换；
将步骤二得到的补零后的四个极化通道图像数据mhh-ins、mhv-ins、mvh-ins、mvv-ins沿每个距离门(按列)进行快速傅里叶变换(fft)，得到方位向频谱数据mhh-fft、mhv-fft、mvh-fft、mvv-fft。
步骤四：方位向信号解压缩，获取等效方位向回波频谱信号；
根据表1给出的雷达系统参数，计算每个距离门(每列)对应的参考斜距rref，以此为参考计算出方位向信号解压缩因子φdecom，利用该解压缩因子乘以步骤三所得的方位向频谱数据mhh-fft、mhv-fft、mvh-fft、mvv-fft得到等效方位向回波频谱信号mhh-decom-fft、mhv-decom-fft、mvh-decom-fft、mvv-decom-fft，具体为：
(1)根据雷达系统参数，计算每个距离门(每列)对应的参考斜距rref：
根据式(14)计算得到的每个距离门(每列)对应的参考斜距rref，式(14)中处于同一距离门(每列)不同方位向(每行)的参考斜距数值是相同的，即每个距离门(每列)采用相同的参考斜距。
(2)雷达的脉冲重复频率prf为120hz，按式(15)计算方位向频谱数据每个方位向(每行)的方位频率fa。
(3)根据上述计算所得每个距离门(每列)对应的参考斜距rref以及每个方位向(每行)的方位频率fa，按式(16)计算得到方位向信号解压缩因子φdecom。
(4)按式(17)(18)(19)(20)将计算得到的方位向信号解压缩因子φdecom乘以方位向频谱数据mhh-fft、mhv-fft、mvh-fft、mvv-fft得到等效方位向回波频谱信号mhh-decom-fft、mhv-decom-fft、mvh-decom-fft、mvv-decom-fft。
步骤五：方位向逆傅里叶变换；
将步骤四得到的等效方位向回波频谱信号mhh-decom-fft、mhv-decom-fft、mvh-decom-fft、mvv-decom-fft沿每个距离门(按列)进行快速逆傅里叶变换(ifft)，得到方位向回波信号mhh-echo、mhv-echo、mvh-echo、mvv-echo。
步骤六：估计法拉第旋转角；
根据步骤五得到的四个极化通道方位向回波信号mhh-echo、mhv-echo、mvh-echo、mvv-echo，通过线性组合得到一组正交圆极化波信号z12和z21，对该组正交圆极化波信号共轭相乘取幅角frtemp，进行平滑处理后得到所估计法拉第旋转角fr。具体为：
(1)对四个极化通道方位向回波信号进行线性组合处理，按式(25)和(26)计算得到一组正交圆极化波信号z12和z21，其均为大小为8792×1024的二维复数矩阵。
(2)按式(27)对该组正交圆极化波信号z12和z21进行共轭相乘并取幅角frtemp。
(3)对上述计算得到的幅角值frtemp(大小为8792×1024的二维实数矩阵)，在仅考虑时变tec的条件下，按式(28)首先将frtemp按列叠加求平均值，再对得到的结果(长度为8792的一维实数数组)采用长度为2m+1＝9的平滑窗进行平滑处理，即以frtemp按列叠加求平均值中以第p项为中心，求前后第p-4到p+4个内的所有元素均值作为第p项值，需要注意的是，对于该数组第1项到第4项，第8789项到8792项，由于平滑处理窗口超出数组边界，因此这些项不做平滑处理，直接以平滑前数值作为该项结果。由此，计算得到法拉第旋转角估计值fr，其为长度为8792的一维实数数组。
步骤七：获取合成孔径时间内的时变tec；
根据步骤六得到的法拉第旋转角估计值fr，利用时变tec与法拉第旋转角间的线性转换关系，按式(29)计算得到合成孔径时间内每个方位时刻的tec，即合成孔径时间内的时变tec，如图3所示。
对比图1所示仿真中原始加入的合成孔径时间内时变tec与根据本发明计算得到的图2所示估计得到的合成孔径时间内时变tec，可以看出两者极为接近。统计采用本发明方法计算得到的时变tec估计误差标准差为0.0287tecu，即可实现高精度时变测量。同时，由于所实施例中时变tec在合成孔径时间内以雷达的脉冲重复周期0.0083s为间隔变化，因此可见本发明提出的测量方法可实现毫秒级的时变tec测量。
综上，上述实施例说明本发明提出的一种基于全极化高轨sar的电离层时变tec测量方法，可以实现时变tec高精度测量。
- 还没有人留言评论。精彩留言会获得点赞！