星载微波辐射计观测亮温的动态滤波重采样方法及装置

摘要

本发明提供了一种星载微波辐射计观测亮温的动态滤波重采样方法及装置，包括根据预设载荷参数建立载荷天线方向图；根据载荷天线方向图、标准载荷参数以及地球基本参数计算预设观测点的观测几何，得到预设观测点的地表投影天线方向图的数据矩阵；根据所述数据矩阵计算预设观测点的滤波系数；根据预设观测点的滤波系数进行关于观测带计数的线性插值计算，得到所有观测点的滤波系数，以便根据所有观测点的滤波系数得到的重采样天线方向图观测亮温。本发明的上述处理方式首先考虑了观测四边形形变的动态滤波，可以有效的减小形变带来的误差，更准确的模拟了真实滤波的情况，提高了所有观测点的滤波系数的准确度且具有业务化应用的潜力。

说明书

技术领域

本发明涉及卫星观测及数据处理领域，具体而言，涉及星载微波辐射计观测亮温的动态滤波重采样方法及装置。

背景技术

星载微波辐射计具有功耗低、体积小、质量轻和工作稳定可靠等特点，应用领域非常宽广。从大的方面来说，星载微波辐射计主要应用于大气探测、海洋观测、对地观测微波遥感3个方面；从具体探测模版来说，星载微波辐射计主要应用于气相、农林、地质、海洋环境检测和军事侦察等，还可用于天文、医疗和导弹的末制导等方面。

星载微波辐射计的观测亮温(BrightnessTemperature)重采样(Resample)技术是指，将不同频率下不同分辨率的观测亮温处理为不同频率下相同分辨率的观测亮温的方法。观测亮温是卫星数据接收单位分发的辐射计产品中最核心的部分，其数值大小依赖于被观测者的物理信息以及观测的天线方向图，由于被观测者的物理信息是客观确定的，而观测的天线方向图是可改变的仪器属性。因此重采样技术的实质就是对天线方向图(通常使用归一化天线方向图)进行一系列计算，使不同分辨率的天线方向图统一。

目前常用的重采样方法，均是使用Backus-Gilbert方法对某一观测点的地表投影天线方向图和该观测点地理位置上邻近的观测点的地表投影天线方向图构成的数据矩阵(在计算机图形学概念中，图和矩阵是等价的)进行计算。这一计算过程的结果包括一组滤波系数和由该系数对数据矩阵进行滤波得到的重采样天线方向图。由于天线方向图和观测亮温的线性一致性，由该组系数能够直接获得重采样的观测亮温。

实际应用中，目前常用的业务化重采样方法是使用特定区域(通常是低纬度地区，如赤道地区)的某一数据矩阵结合Backus-Gilbert方法获得一组确定的滤波系数，用该系数对N个所有的数据矩阵进行滤波处理，然后根据滤波的结果对整体数据进行分析。

发明人在研究中发现，现有技术采用特定区域获得的某一组确定的滤波系数对N个所有的数据矩阵进行滤波处理。这种方法的滤波系数是不变化的，没有考虑不同地理区域的观测几何(观测者和被观测者的相对地理关系)的差异。

发明内容

本发明的目的在于提供一种星载微波辐射计观测亮温的动态滤波重采样方法及装置，通过计算预设观测点的滤波系数，得到所有观测点的滤波系数，用以提高所有观测点的滤波系数的准确度并且本发明的技术具有业务化应用的潜力。

第一方面，本发明实施例提供了一种星载微波辐射计观测亮温的动态滤波重采样方法，包括：

根据预设载荷参数建立载荷天线方向图；所述预设载荷参数包括：载荷的观测频率、载荷的主波束效率、载荷的半功率束宽；

根据所述载荷天线方向图、标准载荷参数以及地球基本参数计算预设观测点的观测几何，得到所述预设观测点的地表投影天线方向图的数据矩阵；

根据所述预设观测点的地表投影天线方向图的数据矩阵，计算所述预设观测点的滤波系数；

根据所述预设观测点的滤波系数进行关于观测带计数的插值计算，得到所有观测点的滤波系数，以便根据所述所有观测点的滤波系数得到的重采样天线方向图观测亮温。

结合第一方面，本发明实施例提供了第一方面的第一种可能的实施方式，其中，所述根据预设载荷参数建立载荷天线方向图包括：

根据天线方向图公式$>F(\theta ,n,C)=\frac{Cf(\theta ,n=0)+\frac{1-C}{n+1}f(\theta ,n)}{C+\frac{1-C}{n+1}}$>建立渐削圆形口径模型；其中，F(θ,n,C)表示天线方向图；f(θ,n)表示天线方向图变换因子且J_n(x)是n阶第一类贝塞尔函数；C是边缘照射且C越小渐削程度越大；n是分布线型且n越大渐削程度越大；a是口径半径；β是平面波的相位常数且空气中θ为立体角分布的天线方向图里的角度自变量；

根据公式$>F(\theta ,n,C)=\frac{Cf(\theta ,n=0)+\frac{1-C}{n+1}f(\theta ,n)}{C+\frac{1-C}{n+1}},$>>和$>{ϵ}_r=\frac{{[C+\frac{1-C}{n+1}]}^2}{C^2+\frac{2C(1-C)}{n+1}+\frac{{(1-C)}^2}{2n+1}}$>计算渐削圆形口径模型与渐削效率、口径半径和边缘照射的第一函数关系；其中，ε_r表示渐削效率；

根据边缘照射C和口径半径a与渐削效率的函数关系，将所述第一函数关系将转换为渐削圆形口径模型与渐削效率的第二函数关系；

根据所述第二函数关系，获取与一个渐削圆形口径模型对应的一个渐削效率和一组模型主波束效率；

将与所述模型主波束效率相同的主波束效率对应的渐削圆形口径模型确定为所述载荷天线方向图。

结合第一方面的第一种可能的实施方式，本发明实施例提供了第一方面的第二种可能的实施方式，其中，所述根据所述载荷天线方向图、标准载荷参数以及地球基本参数计算预设观测点的观测几何，包括：

根据载荷天线方向图、卫星运动的6个开普勒轨道根数和地球半径，计算所述载荷天线方向图在地表的投影；

根据卫星运行速度及卫星的圆锥扫描速度，计算卫星和预设观测点在轨道平面坐标系的坐标；所述预设观测点为一整轨的卫星观测中的一次圆锥扫描的中心点过南极点、升轨赤道、北极点和降轨赤道四个带状区域的观测点；

通过欧拉角旋转的方式，将所述卫星和预设观测点在轨道平面坐标系的坐标，转换为第二赤道坐标系内的坐标；

根据地球的自转角速度和观测时间间隔，对第二赤道坐标系内的赤经进行修正，得到卫星和预设观测点在地理坐标系的坐标；

根据所述卫星和预设观测点在地理坐标系的坐标以及空间坐标与平面坐标的投影规则，计算卫星观测的观测方位角和类正轴等距投影坐标系的坐标；

根据所述载荷天线方向图在地表的投影、卫星的观测方位角和类正轴等距投影坐标系的坐标，计算预设观测点的地表投影天线方向图，以得到所述预设观测点的地表投影天线方向图的数据矩阵。

结合第一方面的第二种可能的实施方式，本发明实施例提供了第一方面的第三种可能的实施方式，其中，所述根据所述预设观测点的地表投影天线方向图的数据矩阵，计算所述预设观测点的滤波系数，包括：

调制Backus-Gilbert方法对应的公式中β从0向无穷大变化，并在首次出现构造天线方向图非负时，计算所述载荷天线方向图中所述预设观测点的滤波系数；其中：V＝G+Eβ；E是单位矩阵，代表一个观测点；G是一个n行n列的矩阵；u、v是n个元素的列向量；且G_ij＝∫G_iG_jdA、u_i＝∫G_idA、v_i＝∫G_iG₀dA；其中，A表示重采样区域；i、j分别表示矩阵元素编号；β表示平滑参数。

结合第一方面的第三种可能的实施方式，本发明实施例提供了第一方面的第四种可能的实施方式，其中，所述根据所述预设观测点的滤波系数进行关于观测带计数的插值计算，得到所有观测点的滤波系数，包括：

将所述预设观测点作为插值结点，沿经向对所述滤波系数进行关于观测带计数的线性插值计算，得到所有观测点的滤波系数。

第二方面，本发明实施例还提供了一种星载微波辐射计观测亮温的动态滤波重采样装置，包括：

建立单元，用于根据预设载荷参数建立载荷天线方向图；所述预设载荷参数包括：载荷的观测频率、载荷的主波束效率、载荷的半功率束宽；

第一计算单元，用于根据所述建立单元建立的所述载荷天线方向图、标准载荷参数以及地球基本参数计算预设观测点的观测几何，得到所述预设观测点的地表投影天线方向图的数据矩阵；

第二计算单元，用于根据所述第一计算单元计算的所述预设观测点的地表投影天线方向图的数据矩阵，计算所述预设观测点的滤波系数；

插值计算处理单元，用于根据所述第二计算单元计算的所述预设观测点的滤波系数进行关于观测带计数的插值计算，得到所有观测点的滤波系数，以便根据所述所有观测点的滤波系数得到的重采样天线方向图观测亮温。

结合第二方面，本发明实施例提供了第二方面的第一种可能的实施方式，其中，所述建立单元包括：

建立子单元，用于根据天线方向图公式建立渐削圆形口径模型；其中，F(θ,n,C)表示天线方向图；f(θ,n)表示天线方向图变换因子且J_n(x)是n阶第一类贝塞尔函数；C是边缘照射且C越小渐削程度越大；n是分布线型且n越大渐削程度越大；a是口径半径；β是平面波的相位常数且空气中θ为立体角分布的天线方向图里的角度自变量；

第一计算子单元，用于根据公式$>F(\theta ,n,C)=\frac{Cf(\theta ,n=0)+\frac{1-C}{n+1}f(\theta ,n)}{C+\frac{1-C}{n+1}},$>>和$>{ϵ}_r=\frac{{[C+\frac{1-C}{n+1}]}^2}{C^2+\frac{2C(1-C)}{n+1}+\frac{{(1-C)}^2}{2n+1}}$>计算所述建立子单元建立的所述渐削圆形口径模型与渐削效率、口径半径和边缘照射的第一函数关系；其中，ε_r表示渐削效率；

第一转换子单元，用于根据边缘照射C和口径半径a与渐削效率的函数关系，将所述第一计算子单元计算得到的所述第一函数关系将转换为渐削圆形口径模型与渐削效率的第二函数关系；

获取子单元，用于根据所述第一转换子单元转换得到的所述第二函数关系，获取与一个渐削圆形口径模型对应的一个渐削效率和一组模型主波束效率；

确定子单元，用于将与所述获取子单元获取的所述模型主波束效率相同的主波束效率对应的渐削圆形口径模型确定为所述载荷天线方向图。

结合第二方面的第一种可能的实施方式，本发明实施例提供了第二方面的第二种可能的实施方式，其中，所述第一计算单元包括：

第二计算子单元，用于根据载荷天线方向图、卫星运动的6个开普勒轨道根数和地球半径，计算所述载荷天线方向图在地表的投影；

第三计算子单元，用于根据卫星运行速度及卫星的圆锥扫描速度，计算卫星和预设观测点在轨道平面坐标系的坐标；所述预设观测点为一整轨的卫星观测中的一次圆锥扫描的中心点过南极点、升轨赤道、北极点和降轨赤道四个带状区域的观测点；

第二转换子单元，用于通过欧拉角旋转的方式，将所述卫星和预设观测点在轨道平面坐标系的坐标，转换为第二赤道坐标系内的坐标；

修正子单元，用于根据地球的自转角速度和观测时间间隔，对第二赤道坐标系内的赤经进行修正，得到卫星和预设观测点在地理坐标系的坐标；

第四计算子单元，用于根据所述卫星和预设观测点在地理坐标系的坐标以及空间坐标与平面坐标的投影规则，计算卫星观测的观测方位角和类正轴等距投影坐标系的坐标；

第五计算子单元，用于根据所述载荷天线方向图在地表的投影、卫星的观测方位角和类正轴等距投影坐标系的坐标，计算预设观测点的地表投影天线方向图，以得到所述预设观测点的地表投影天线方向图的数据矩阵。

结合第二方面的第二种可能的实施方式，本发明实施例提供了第二方面的第三种可能的实施方式，其中，所述第二计算单元包括：

调制子单元，用于调制Backus-Gilbert方法对应的公式中β从0向无穷大变化；

第六计算子单元，用于在所述调制子单元在调制过程中首次出现构造天线方向图非负时，计算所述载荷天线方向图中所述预设观测点的滤波系数；其中：V＝G+Eβ；E是单位矩阵，代表一个观测点；G是一个n行n列的矩阵；u、v是n个元素的列向量；且G_ij＝∫G_iG_jdA、u_i＝∫G_idA、v_i＝∫G_iG₀dA；其中，A表示重采样区域；i、j分别表示矩阵元素编号；β表示平滑参数。

结合第二方面的第三种可能的实施方式，本发明实施例提供了第二方面的第四种可能的实施方式，其中，所述插值计算处理单元包括：

设置子单元，用于将所述预设观测点设置为插值结点；

线性插值计算子单元，用于根据所述设置子单元设置的所述插值结点沿经向对所述滤波系数进行关于观测带计数的线性插值计算，得到所有观测点的滤波系数。

本发明实施例提供的星载微波辐射计观测亮温的动态滤波重采样方法及装置，其根据预设载荷参数建立载荷天线方向图，通过建立的载荷天线方向图、标准载荷参数以及地球基本参数计算预设观测点的观测几何，得到预设观测点的地表投影天线方向图的数据矩阵；然后根据得到的数据矩阵计算预设观测点的滤波系数，最后对根据预设观测点的滤波系数进行关于观测带计数的线性插值计算，得到所有观测点的滤波系数，以便根据所有观测点的滤波系数得到的重采样天线方向图观测亮温。本发明的上述处理方式考虑了观测四边形形变，有效的减小形变带来的误差，并且通过计算预设观测点的滤波系数，计算所有观测点的滤波系数，提高了所有观测点的滤波系数的准确度并且本发明的技术具有业务化应用的潜力。

为使本发明的上述目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附附图，作详细说明如下。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

图1示出了本发明实施例所提供的一种星载微波辐射计观测亮温的动态滤波重采样方法的流程图；

图2示出了本发明实施例所提供的一组29×29观测点以及其观测时星下点在类等距正轴投影坐标系中的坐标的结构示意图；

图3示出了本发明实施例所提供的星载微波辐射计圆锥扫描模式示意图；

图4示出了本发明实施例所提供的卫星位置与观测位置的相互关系示意图；

图5示出了本发明实施例所提供的升轨和降轨的示意图；其中，a、b、c、d为选则的四个插值节点区域；e、f、g分别是由每次圆锥扫描观测起始点、中间点和结束点分别构成的数列，且上述升轨为从南极到北极的轨道，上述降轨为从北极到南极的轨道；

图6示出了本发明实施例所提供的一组观测点之间的位置关系示意图；

图7示出了本发明实施例所提供的滤波系数的变化规律示意图；其中，滤波系数是基于图5中的e、f、g三个数列分别进行计算得到的，且a、b、c、d和e分别对应图6中的point113、point1、point15、point211和point255的位置；

图8示出了本发明实施例所提供的预设载荷在一次圆锥扫描中(如图5所示的a、b、c、d四个插值节点区域)以及在升轨和降轨过程中(如图5所示的e、f、g三个数列)，本发明的动态滤波的结果和传统动态滤波结果的对比示意图；

图9示出了本发明实施例所提供的一种星载微波辐射计观测亮温的动态滤波重采样装置的结构示意图；

图10示出了本发明实施例所提供的一种星载微波辐射计观测亮温的动态滤波重采样装置中监测单元和第一计算单元的结构示意图；

图11示出了本发明实施例所提供的一种星载微波辐射计观测亮温的动态滤波重采样装置中第二计算单元和第一计算单元的结构示意图；

主要元件符号说明：

11、建立单元；12、第一计算单元；13、第二计算单元；14、线性插值处理单元；111、建立子单元；112、第一计算子单元；113、第一转换子单元；114、获取子单元；115、确定子单元；121、第二计算子单元；122、第三计算子单元；123、第二转换子单元；124、修正子单元；125、第四计算子单元；126、第五计算子单元；131、调制子单元；132、第六计算子单元；141、设置子单元；142、线性插值计算子单元142。

具体实施方式

下面将结合本发明实施例中附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

辐射计主要用于测量地球表面的物质通过热辐射辐射出的电磁波，其是备用遥感传感器，只接收海面或大气的辐射，从中提取物理信息，而不发射探测电磁波。辐射计主要用于获取海面温度，海面风速，大气中的水汽含量，云液态水含量等特定信息。目前，常用的辐射计为微波辐射计，将其搭载在载荷上(优选为卫星)，用于观测上述特定信息的亮温数据，然后对所述亮温数据进行反演处理，以获得准确的海洋大气参数。

考虑到现有技术中，星载微波辐射计的观测亮温(BrightnessTemperature)重采样(Resample)技术中，常用的业务化重采样方法是使用特定区域(通常是低纬度地区，如赤道地区)的某一数据矩阵结合Backus-Gilbert方法获得一组确定的滤波系数，用该系数对N个所有的数据矩阵进行滤波处理。而该方法中滤波系数是不变的，其并未考虑不同地理区域的观测几何(观测者和被观测者的相对地理关系)的差异，如地球自转、卫星运动与观测方式等因素带来的观测几何变化，使得滤波系数的误差较大，以致于后续不能准确的建立天线方向图以及进行亮温观测。基于此，本发明实施例提供了一种星载微波辐射计观测亮温的动态滤波重采样方法及装置，下面通过实施例进行描述。

本发明提供了一种星载微波辐射计观测亮温的动态滤波重采样方法，该方法可以通过计算机等终端设备执行，该终端设备上预先装有相应的处理软件，如Matlab或Fortran等编程软件。参见图1所示的星载微波辐射计观测亮温的动态滤波重采样流程图，该方法包括以下步骤：

步骤S101、根据预设载荷参数建立载荷天线方向图；所述预设载荷参数包括：载荷的观测频率、载荷的主波束效率和载荷的半功率束宽。

首先，考虑到现有技术获得观测几何的方法，通常是直接读取卫星观测数据中特定区域的经纬度值和观测方向，而由于现有技术中直接获得的数据携带有大量的隐函数，无法对所述隐函数进行分析判断，导致通过该方法获得的数据在某一个位置(即获取的数据对应的位置)的精度准确，而在所有其他位置的精度均较差(即偏离选取位置的其他位置的精度较差)，进而不方便后续建模处理及理论分析研究的问题，即现有技术直接获取数据的方法存在泛用性较差的问题)；故为了解决该问题，本实施例中，根据载荷的观测频率、载荷的主波束效率和载荷的半功率束宽等预设载荷参数建立天线方向图，即建立卫星观测模型，能够直观的模拟星载微波辐射计。

其中，本实施例中的预设载荷为卫星及其搭载的微波辐射计这个整体，用于观测亮温数据；而所述预设载荷参数则包括：载荷的观测频率为搭载在卫星上的微波辐射计上的观测频率(即载荷的观测频率)，搭载在卫星上的微波辐射计的主波束效率(即载荷的主波束效率)，以及搭载在卫星上的微波辐射计的半功率束宽(即载荷的半功率束宽)。

本实施例中的星载微波辐射计(AMSR-E)的的典型参数如下表表1和表2所示：

轨道模式太阳同步回归卫星高度(km)705半长轴a(km)7072轨道偏心率e0轨道倾角i98.15°近心点幅角ω0轨道周期T(min)99扫描周期t(s)1.5扫描一周内观测间隔243(观测)/581(总共)

表1

表2

步骤S102、根据上述载荷天线方向图、标准载荷参数以及地球基本参数计算预设观测点的观测几何，得到预设观测点的地表投影天线方向图的数据矩阵。

具体的，所述数据矩阵表示不同预设观测点的地表投影天线方向图之间的位置关系和数量关系；其中，卫星数据通常采用矩阵的形式进行保存，比如半轨2000×234的数据，2000可以视为经向数据(近似沿卫星飞行方向分布)，243可以视为纬向数据(沿一次圆锥扫描的扫描点的分布，近似纬向)。

所述标准载荷参数包括卫星运动的6个开普勒轨道根数、单次观测时间和圆锥扫描速度(其中，星载微波辐射计使用圆锥扫描雷达的扫描方式，完成在角度上对目标自动跟踪的雷达)；所述地球基本参数包括：地球半径和地球自转角速度。

优选的，本实施例中的预设观测点选取一整轨的卫星观测中的一次圆锥扫描的中心点过南极点、升轨赤道、北极点、降轨赤道四个带状区域的观测点。其中，上述一整轨总共4000个左右观测带(本实施例中选择上述四个区域的观测带即可)，而一条观测带有243个观测点。

需要说明的是，上述四个观测区域为具有代表性的观测区域，同样选择上述四个区域的预设观测带和对应的观测点仅仅是本发明实施例中选择的一种优选的选择，而本发明中的观测区域并不仅限制于上述四个位置，还可以为8个(即包括这四个位置以及这个四个位置中间的四个位置)或10、12个等。

步骤S103、根据得到的预设观测点的地表投影天线方向图的数据矩阵，计算该预设观测点的滤波系数。

本实施例中，优选通过Backus-Gilbert(连续介质的反演理论)方法简称BG方法对预设观测点的数据矩阵进行计算，得到预设观测点的滤波系数。其中，BG方法是由Backus和Gilbert建立的，目前已形成一套完整、系统的理论。

本实施例中，将BG方法应用到设置的上述预设观测点计算得到数据矩阵中，得到上述四个位置(南极点、升轨赤道、北极点和降轨赤道)的观测点的滤波系数，具体计算过程如下：

1、选取某一圆锥扫描观测带(即上述四个区域中任意一个区域中的观测带，其中，每个观测带包含243个观测点)，及其前后各14个圆锥扫描观测带的地表投影天线方向图构成数据矩阵，共243×29个地表投影天线方向图；

2、依据选取的圆锥扫描观测带，逐次计算上述243个数据点(即243个观测点)的滤波系数。其中，计算上述每一个观测点时的方法如下：

A.建立平面坐标系。具体的，将观测点按正轴等距投影至平面坐标原点，正北方向为y轴正向，圆锥扫描方向为x轴负向。

B.计算以该观测点为中心，29×29个数据点在该平面坐标系中的坐标和观测方向。

C.根据BG方法计算上述四个位置(南极点、升轨赤道、北极点和降轨赤道)的观测点的数滤波系数。

D.注意，圆锥扫描带起始和结束的各14个位置的数据矩阵不是29×29。在某一维度上会相应减少。例如第一个和第二个点为15×29、16×29，倒数第二个和倒数第一个点为16×29、15×29。

需要说明的是，BG方法仅仅是本发明实施例中计算滤波系数的一种优选方法，而本发明中的滤波方法并不仅限制于该方法，如还可以为图像反卷积技术和散射计图像重建算法等。

步骤S104、根据上述预设观测点的滤波系数进行关于观测带计数的插值计算，得到所有观测点的滤波系数，以便根据所有观测点的滤波系数得到的重采样天线方向图观测亮温。

本实施例中的插值计算则是根据上述四个预设观测点进行线性插值计算，得到所有观测点的滤波系数；其中，线性插值方法是插值方法的一种，这种方法最简便。如果增加插值节点(分8段)或采用别的插值函数(如样条插值等)也可以得到近似的结果，故本实施例对插值方法和插值预分的段数不做具体限制。

具体的，在得到所有观测点的滤波系数以后，根据该所有观测点的滤波系数，对所有观测点的观测几何(即所有观测点的地表投影天线方向图的数据矩阵)进行滤波得到重采样天线方向图，并根据重采样天线方向图观测亮温数据。通过该方法，较准确的得到了除了预测观测点之外的其他观测点的滤波系数，且通过该滤波系数得到的天线方向图更准确的模拟了真实滤波的情况，从而使得根据重采样天线方向图观测的亮温数据更加真实可靠。

本实施例中，基于Backus-Gilbert算法的动态滤波的观测亮温重采样算法，能够有效地降低观测几何的变化带来的误差。

本发明实施例提供的星载微波辐射计观测亮温的动态滤波重采样方法，其根据预设载荷参数建立载荷天线方向图，通过建立的载荷天线方向图、标准载荷参数以及地球基本参数计算预设观测点的观测几何，得到预设观测点的地表投影天线方向图的数据矩阵；然后根据得到的数据矩阵计算预设观测点的滤波系数，最后对根据预设观测点的滤波系数进行关于观测带计数的线性插值计算，得到所有观测点的滤波系数，以便根据所有观测点的滤波系数得到的重采样天线方向图观测亮温。本发明的上述处理方式考虑了观测四边形形变，有效的减小形变带来的误差，并且通过计算预设观测点的滤波系数，计算所有观测点的滤波系数，提高了所有观测点的滤波系数的准确度并且本发明的技术具有业务化应用的潜力。

考虑到现有技术中直接获取观测几何的方法，具有不方便后续建模处理及理论分析研究的问题(即存在泛用性较差的问题)，本发明实施例中首先要建立天线方向图，本发明实施例中优选为立体角分布的天线方向图，以更好模拟星载微波辐射计的观测模型；实际中，星载微波辐射计的反射器系统决定了天线方向图的分布，本发明中通过建立渐削圆形口径模型，用以较好的模拟星载微波辐射计的反射器系统。基于此，本实施例中的天线方向图的建立过程如下：

首先根据天线方向图公式$>F(\theta ,n,C)=\frac{Cf(\theta ,n=0)+\frac{1-C}{n+1}f(\theta ,n)}{C+\frac{1-C}{n+1}},$>建立渐削圆形口径模型；其中，所述建立渐削圆形口径模型，用以较好的模拟星载微波辐射计的反射器系统，以决定了天线方向图的分布。上述天线方向图公式中，F(θ,n,C)表示天线方向图；f(θ,n)表示天线方向图变换因子且J_n(x)是n阶第一类贝塞尔函数；C是边缘照射，C越小渐削程度越大；n是分布线型，n越大渐削程度越大；a是口径半径；β是平面波的相位常数，空气中θ为立体角分布的天线方向图里的角度自变量；然后，根据第一函数关系将所述渐削圆形口径模型表示为渐削效率的第二函数(即根据边缘照射C和口径半径a与渐削效率的函数关系，将第一函数关系将转换为渐削圆形口径模型与渐削效率的第二函数关系)，上述第一函数关系为天线方向图与渐削效率、口径半径和边缘照射的函数关系，且上述第一函数关系是根据公式>和计算得到的；其中，ε_r表示渐削效率；边缘照射C和口径半径a可以根据所述观测频率和半功率束宽通过公式用渐削效率代替，故上述函数关系则表示为天线方向图与渐削效率的函数关系；于是，由一个渐削效率就能获得一个天线方向图，进而得到一组模型主波束效率。最后，根据所述主波束效率和上述步骤得到的模型主波束效率，找到二者相同时的天线方向图即为所述载荷天线方向图，即将与所述模型主波束效率相同的主波束效率对应的天线方向图确定为所述载荷天线方向图。

本实施例中，要计算出预设观测点的地表投影天线方向图，以获取预设观测点的地表投影天线方向图的数据矩阵，首先需要根据载荷天线方向图、卫星运动的6个开普勒轨道根数和地球半径，计算天线方向图在地表的投影。

而为了方便观测亮温重采样的计算，并使地图投影误差最小，优选采用等距正轴投影的一种变形方式：选取某一观测点为图像中心，经圈仍为子午圈，纬圈为过图像中心的大圆，其余和正轴等距投影相同。重采样选取300km×300km的范围足够，将这样的区域的投影中心规定为类似标准世界地图的赤道位置，其正轴等距投影的图像变形是正轴等距投影中最小的。图2示出了一组29×29观测点以及其观测时星下点在类等距正轴投影坐标系中的坐标的示意图。

而要想得到上述类正轴等距投影坐标系的坐标，需要首先计算卫星和预设观测点在轨道平面坐标系的坐标，在将上述坐标进行一系列转换，即可得到卫星和预设观测点类正轴等距投影坐标系的坐标。

转换的具体过程如下：首先，根据卫星运行速度及卫星的圆锥扫描速度，计算卫星和预设观测点在轨道平面坐标系的坐标；其中，该预设观测点为一整轨的卫星观测中的一次圆锥扫描的中心点过南极点、升轨赤道、北极点、降轨赤道4个带状区域的观测点。然后，通过欧拉角旋转的方式，将卫星和预设观测点在轨道平面坐标系的坐标，转换为第二赤道坐标系内的坐标；在根据地球的自转角速度和观测时间间隔，对第二赤道坐标系内的赤经进行修正，得到卫星和预设观测点在地理坐标系的坐标；最后，根据卫星和预设观测点在地理坐标系的坐标以及空间坐标与平面坐标的投影规则(其中，该投影规则是地图学中的概念，是一种将三维空间坐标变换为二维平面坐标的方法，为公知技术，本发明对比不在详细说明)，计算卫星观测的观测方位角和类正轴等距投影坐标系的坐标。通过上述转换方式，即可得到卫星的观测方位角和卫星和预设观测点的类正轴等距投影坐标系的坐标。

如图3所示，轨道平面坐标系：轨道平面坐标系的“平面”二字，意思是卫星的运动轨迹在一个平面内。考虑了星载辐射计的圆锥扫描模式，卫星运动的平面坐标可以方便的转写为三维空间坐标系的模式。这种模式可以认为是两种运动的叠加：卫星的前进运动和辐射计的圆锥扫描。这需要设定卫星运动的速度以及圆锥扫描的速度以及方式。

天球坐标系：关于地球卫星的轨道研究，通常采用第二赤道坐标系。卫星位置和观测点在轨道平面坐标系内坐标可以方便的通过欧拉角旋转，变为第二赤道坐标系内的坐标。这个过程需要设定卫星的6个开普勒轨道根数。欧拉角旋转可以方便的描述卫星轨道平面在天球坐标系中的变换：第一欧拉角即为升交点经度，第二欧拉角即为轨道倾角，第三欧拉角即为近心点幅角。

地理坐标系：包含地球经纬度的坐标系称为地理坐标系。其坐标和第二赤道坐标系中的赤经和赤纬仅仅相差一个常数。将第二赤道坐标系中的坐标转变为地理坐标系的方法，需要在赤经上进行关于地球自转角速度和观测时间间隔的修正。

图4示出了本发明实施例所提供的卫星位置与观测位置的相互关系示意图；如图4所示，建立空间坐标系xyz和平面坐标系xoy。卫星观测时，通常保持不变的是天底角α。由于轨道平面和地球曲率以及地形的变化，观测的入射角β并不是固定的。假设地球半径R为6367km，卫星高度H(开普勒根数中的轨道半径与地球半径之差)为705km。已知三角形SoC中，CS＝R+H，oC＝R，α＝47.5°，则三角形SoC确定。这时β＝55.0°。地表P的增益G_P(ρ_P＝SP)，可以通过SP与So的夹角和天线方向图的立体角分布求出。

其中，观测方位角的确定过程如下：准确的观测方向是如图2中所示的向量So在地球表面的投影。星下点在地理坐标系中的投影位置和观测点位置之间是有一定距离的。受限于地图投影的规则，直接使用这两点计算观测方向并不可靠。实际观测方向与经圈切向的夹角越大，投影规则带来的误差就越明显。在向量So上选取距o非常近的一点，用这点在投影平面的投影点计算的观测方位角更为准确。

最后，根据得到的天线方向图在地表的投影、类正轴等距投影坐标系的坐标以及卫星观测的观测方位角，计算预设观测点的地表投影天线方向图，以得到预设观测点的地表投影天线方向图的数据矩阵。

目前重采样的业务化算法大多使用一个固定数值的滤波窗口对高分辨率亮温数据进行滤波处理。这样的方法计算简便，却忽略了不同的观测足印构成的观测四边形的形变造成的误差。而本发明实施例使用考虑了观测四边形形变的动态滤波，可以有效的减小形变带来的误差。为了方便起见，我们约定了计数规则，展示在图5和图6中。其中，图5示出了本发明实施例所提供的升轨和降轨的示意图；其中，a、b、c、d为选则的四个插值节点区域；e、f、g分别是由每次圆锥扫描观测起始点、中间点和结束点分别构成的数列，且上述升轨为从南极到北极的轨道，上述降轨由于卫星运动方向从南-北变为北-南，其观测计数是升轨观测计数关于观测中心Point113做中心反射变换；在图5中，左侧为升轨，右侧为降轨；图6示出了本发明实施例所提供的一组观测点之间的位置关系示意图。

在使用Backus-Gilbert方法进行重采样计算时，对平滑参数的选取有这样的规则：平滑参数的增大导致了重采样结果的变坏。结合天线方向图的非负特性，一个合理的选定平滑参数的方法是：当平滑参数从0向无穷大变化时，首次出现构造天线方向图非负的平滑参数即为最佳平滑参数。在这样的最佳平滑参数的选择下，能过获得一组滤波系数。

根据这样的规则，我们计算了具有代表性的7个区域(a-e)的滤波系数。滤波系数的纬向分布由于受到地球自转和圆锥扫描投影的双重影响，并没有明显的规律。滤波系数的经向分布主要受地球自转的影响，规律性较明显。图7中横坐标计数方式是从南极开始的一整轨观测计数，1000为升轨赤道附近，2000为北极区域，3000为降轨赤道附近，4000为南极区域。在图7中，直线表示图6中obs122点在上述1000、2000和3000位置的中滤波系数的变化情况；点划线表示图6中的obs18点在上述1000、2000和3000位置的中滤波系数的变化情况；细虚线表示图6中的obs236点在上述1000、2000和3000位置的中滤波系数的变化情况。

因此，设计了这样的滤波窗口：以升轨赤道、北极、降轨赤道、南极四个观测带为插值节点，沿经向对系数进行关于观测带计数的线性插值。这主要是考虑了滤波系数纬向的变化规律不如经向的变化规律简单。图7中的上b图和上c图看起来不那么“平滑”，是因为图像轴范围较窄和计算最优平滑参数时算法的最小间隔。

平滑参数是BG方法中的一个可以调节的数，具体如下：

BG方法的滤波系数表示为如下形式：

其中：V＝G+Eβ，E是单位矩阵。滤波范围内一个矩阵元素代表一个观测点。通常选取矩阵的行列数一致(假设为n)。G是一个n行n列的矩阵，u、v是n个元素的列向量。G_ij＝∫G_iG_jdA、u_i＝∫G_idA、v_i＝∫G_iG₀dA。A表示重采样区域；角标i、j表示矩阵元素编号。

对于确定的观测几何，G、u、v都是确定的。而β是平滑参数，选取范围0到无穷大。

根据Backus-Gilbert方法对应的公式计算滤波系数；其中：V＝G+Eβ，E是单位矩阵。滤波范围内一个矩阵元素代表一个观测点；通常选取矩阵的行列数一致(假设为n)；G是一个n行n列的矩阵，u、v是n个元素的列向量；G_ij＝∫G_iG_jdA、u_i＝∫G_idA、v_i＝∫G_iG₀dA；A表示重采样区域；角标i、j表示矩阵元素编号；β表示平滑参数；

调制Backus-Gilbert方法中的平滑参数β从0向无穷大变化，并在首次出现构造天线方向图非负时，计算天线方向图的滤波系数。在计算了上述四个预设观测点的位置的滤波系数后，还需要沿经向对上述代表性的四个预设观测点的滤波系数进行关于观测带计数的线性插值计算，得到重采样天线方向图；其中，插值节点为预设区域的预设观测点。其中，用滤波系数进行滤波计算(每一个系数与对应位置的天线方向图相乘，再进行累加运算)，所得到的结果就是构造天线方向图。

参考图5，图5中沿着f，在a、b中间的某一点p处的滤波系数可以根据下式计算：其中，大写字母P、A、B表示该位置的滤波系数，其中A、B是已知的。ab、pb、pa是两点间的圆锥扫描带计数间隔，这可以从卫星数据的保存格式中方便的获得。这就是一种线性插值的表示方式。

其中，P、A、B分别表示对应位置的滤波系数，其中，A、B是已知常量；ab、pb和pa是两点间的圆锥扫描带计数间隔。

图8示出了本发明实施例所提供的预设载荷在一次圆锥扫描中(如图5所示的a、b、c、d四个插值节点区域)以及在升轨和降轨过程中(如图5所示的e、f、g三个数列)，本发明的动态滤波的结果和传统动态滤波结果的对比示意图，即本发明的技术和传统滤波重采样技术的结果比较，一个量化的方式是比较重采样之后的天线方向图和重采样目标的天线方向图之间的差异。图8中是使用传统方法——即固定系数的滤波(如图8中的虚线所示)以及动态滤波(即如图8中的实线所示)对7个具有典型特征的区域的重采样结果。传统方法在选取的滤波系数的区域的周围区域，结果较好。距离这一区域越远，效果越差。

对应上述方法，本发明还提供了一种星载微波辐射计观测亮温的动态滤波重采样装置，所述动态滤波重采样装置可运行于计算机等终端设备，所述终端设备上设置有Matlab或Fortran等编程软件，用于实现星载微波辐射计观测亮温的动态滤波重采样计算，参考图9，所述装置包括包括：

建立单元11，用于根据预设载荷参数建立载荷天线方向图；预设载荷参数包括：载荷的观测频率、载荷的主波束效率、载荷的半功率束宽；

第一计算单元12，用于根据建立单元11建立的载荷天线方向图、标准载荷参数以及地球基本参数计算预设观测点的观测几何，得到预设观测点的地表投影天线方向图的数据矩阵；

第二计算单元13，用于根据第一计算单元12计算的预设观测点的地表投影天线方向图的数据矩阵，计算预设观测点的滤波系数；

插值计算处理单元14，用于根据第二计算单元13计算的预设观测点的滤波系数进行关于观测带计数的线性插值计算，得到所有观测点的滤波系数，以便根据所有观测点的滤波系数得到的重采样天线方向图观测亮温。

本实施例为防止直接获取观测几何的方法，由于直接获得的数据携带有大量的隐函数，无法对隐函数进行分析判断，导致通过该方法获得的数据在某一个位置(即获取的数据对应的位置)的精度准确，而在所有其他位置的精度均较差(即偏离所有其他位置较大)，进而使得不方便后续建模处理及理论分析研究的问题(即泛用性较差的问题)，本发明采取了建立天线方向图模拟星载微波辐射计的方式，参考图10，本实施例中，建立单元11包括：

建立子单元111，用于根据天线方向图公式建立渐削圆形口径模型；其中，F(θ,n,C)表示天线方向图；f(θ,n)表示天线方向图变换因子且J_n(x)是n阶第一类贝塞尔函数；C是边缘照射且C越小渐削程度越大；n是分布线型且n越大渐削程度越大；a是口径半径；β是平面波的相位常数且空气中θ为立体角分布的天线方向图里的角度自变量；

第一计算子单元112，用于根据公式>和计算建立子单元111建立的渐削圆形口径模型与渐削效率、口径半径和边缘照射的第一函数关系；其中，ε_r表示渐削效率；

第一转换子单元113，用于根据边缘照射C和口径半径a与渐削效率的函数关系，将第一计算子单元112计算得到的第一函数关系将转换为渐削圆形口径模型与渐削效率的第二函数关系；

获取子单元114，用于根据第一转换子单元113转换得到的第二函数关系，获取与一个渐削圆形口径模型对应的一个渐削效率和一组模型主波束效率；

确定子单元115，用于将与获取子单元114获取的模型主波束效率相同的主波束效率对应的渐削圆形口径模型确定为载荷天线方向图。

进一步的，参考图10，第一计算单元12包括：

第二计算子单元121，用于根据载荷天线方向图、卫星运动的6个开普勒轨道根数和地球半径，计算载荷天线方向图在地表的投影；

第三计算子单元122，用于根据卫星运行速度及卫星的圆锥扫描速度，计算卫星和预设观测点在轨道平面坐标系的坐标；预设观测点为一整轨的卫星观测中的一次圆锥扫描的中心点过南极点、升轨赤道、北极点和降轨赤道四个带状区域的观测点；

第二转换子单元123，用于通过欧拉角旋转的方式，将卫星和预设观测点在轨道平面坐标系的坐标，转换为第二赤道坐标系内的坐标；

修正子单元124，用于根据地球的自转角速度和观测时间间隔，对第二赤道坐标系内的赤经进行修正，得到卫星和预设观测点在地理坐标系的坐标；

第四计算子单元125，用于根据卫星和预设观测点在地理坐标系的坐标以及空间坐标与平面坐标的投影规则，计算卫星观测的观测方位角和类正轴等距投影坐标系的坐标；

第五计算子单元126，用于根据载荷天线方向图在地表的投影、卫星的观测方位角和类正轴等距投影坐标系的坐标，计算预设观测点的地表投影天线方向图，以得到预设观测点的地表投影天线方向图的数据矩阵。

进一步的，参考图11，第二计算单元13包括：

调制子单元131，用于调制Backus-Gilbert方法对应的公式中β从0向无穷大变化；

第六计算子单元132，用于在调制子单元在调制过程中首次出现构造天线方向图非负时，计算载荷天线方向图中预设观测点的滤波系数；其中：V＝G+Eβ；E是单位矩阵，代表一个观测点；G是一个n行n列的矩阵；u、v是n个元素的列向量；且G_ij＝∫G_iG_jdA、u_i＝∫G_idA、v_i＝∫G_iG₀dA；其中，A表示重采样区域；i、j分别表示矩阵元素编号；β表示平滑参数。

进一步的，参考图11，插值计算处理单元14，包括：

设置子单元141，用于将预设观测点设置为插值结点；

线性插值计算子单元142，用于根据设置子单元设置的插值结点沿经向对滤波系数进行关于观测带计数的线性插值计算，得到所有观测点的滤波系数。

本发明实施例提供的星载微波辐射计观测亮温的动态滤波重采样装置，其根据预设载荷参数建立载荷天线方向图，通过建立的载荷天线方向图、标准载荷参数以及地球基本参数计算预设观测点的观测几何，得到预设观测点的地表投影天线方向图的数据矩阵；然后根据得到的数据矩阵计算预设观测点的滤波系数，最后对根据预设观测点的滤波系数进行关于观测带计数的线性插值计算，得到所有观测点的滤波系数，以便根据所有观测点的滤波系数得到的重采样天线方向图观测亮温。本发明的上述处理方式考虑了观测四边形形变，有效的减小形变带来的误差，并且通过计算预设观测点的滤波系数，计算所有观测点的滤波系数，提高了所有观测点的滤波系数的准确度并且本发明的技术具有业务化应用的潜力。

本发明实施例所提供的进行星载微波辐射计观测亮温的动态滤波重采样方法的计算机程序产品，包括存储了程序代码的计算机可读存储介质，所述程序代码包括的指令可用于执行前面方法实施例中所述的方法，具体实现可参见方法实施例，在此不再赘述。

所属领域的技术人员可以清楚地了解到，为描述的方便和简洁，上述描述的系统、装置和单元的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。

在本申请所提供的几个实施例中，应该理解到，所揭露的系统、装置和方法，可以通过其它的方式实现。以上所描述的装置实施例仅仅是示意性的，例如，所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，又例如，多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些通信接口，装置或单元的间接耦合或通信连接，可以是电性，机械或其它的形式。

所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。

另外，在本发明各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。

所述功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-OnlyMemory)、随机存取存储器(RAM，RandomAccessMemory)、磁碟或者光盘等各种可以存储程序代码的介质。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。