CN105844000A - 一种mcc 表面海流反演方法 - Google Patents

Abstract

本发明涉及一种MCC表面海流反演方法，所述方法包括对MODIS卫星进行预处理；利用综合优化云检测算法，进行云覆盖检测；对于晴空区，分别选取海表温度和550nm反射率数据为示踪物，采用MMC反演得到表面海流场；对于缺测的海表流数据，利用DINEOF方法进行插值补缺；进行变分同化优化调整获得短时平均表面流。本发明的有益效果为：本发明提出的云检测方法，提高了海面低云、碎云以及部分被云覆盖像元的识别准确率，为下一步计算表面海流更精确的去除了干扰，保证了后续结果的精度。本发明在插值基础上，提出的一种表面海流场全局优化调整方法，大大增加了反演结果的覆盖范围，同时也使得反演精度得到明显提高，增强了该方法的使用范围。

Description

一种MCC表面海流反演方法

技术领域

本发明属于卫星遥感应用领域，具体涉及一种MCC表面海流反演方法。

背景技术

海流对海洋中多种物理、化学、生物和地质过程，以及海洋上空的气候和天气的形成及变化，都有影响和制约作用，故了解和掌握海流规律，对渔业、航运和排污等都有重要意义。例如：海流使低纬度热量向高纬度传输，暖流上空有热量和水汽向上输送，使得层结不稳定、空气湿度增大而易产生降水，寒流产生逆温，层结稳定，水汽不易向上输送，蒸发又弱，下层相对湿度有时虽然很大，但只能成雾，不能成雨；寒暖流交汇的海区，海水受到扰动，可以将下层营养盐类带到表层，有利于鱼类大量繁殖，为鱼类提供诱饵；两种海流还可以形成“水障”，阻碍鱼类活动，使得鱼群集中，往往形成较大的渔场；海流还可以使污染物因迅速扩散而加快其稀释和净化的速度，也相应地使污染范围扩大。

海流在军事上也有重要的意义。对于海军作战，海流是考虑作战的重要因素之一，恰当地利用海流，对作战会起到降低战争消耗，增加胜利因素的一个筹码，反之，会对战争起到不利影响，甚至会导致失败。在第二次世界大战中有这方面的成功的战例。在战争中，各沿海国都把海洋水文资料作为军事情报进行收集、整理，了解海洋水文状况对于海军进行防御、进攻起着重要作用。舰艇顺水航行，有利于提高航行速度，节约时间和燃料；下降流利于潜艇的下潜，上升流利于潜艇的出击；贸易航线因船只多，容易暴露军事目标；风浪大的航线(特别是南北半球的西风漂流)利于偷袭等等。

目前利用的常规手段对大面积的海域的表面海流实施高频率实时观测几乎是不可能的，例如，江苏沿海类似908大范围的海洋调查相隔30年一次，局部地区5-10年一次，重点港口2-3年一次。其他非常规手段由于种种原因难以广泛使用，例如：SAR覆盖范围小，重复观测周期长；岸基雷达成本高，范围有限；高度计资料仅适用于大范围的地转流，且同样存在重复观测周期长的问题。而被动的卫星遥感影像资料具有重复观测周期短、覆盖范围大、资料容易获取等优点，因此，利用卫星影像资料进行表面海流遥感反演，具有非常重要理论和实践意义。

目前，最大相关系数方法(MCC方法)是利用卫星影像资料进行表面海流遥感反演的主要手段。

最初MCC法是用在AVHRR图像的亮温上来获得表面海流。后来，一些学者用SST(海表温度)或SST的水平梯度来代替亮温作为示踪物得到流速。Garciaand Robinson和Tokmakian将MCC方法用到CZCS的可见光图像上，而不是红外、近红外图像来反演海表面流场。Isern-Fontanet利用微波反演的SST跟踪海流。Crocker分别利用AVHRR反演的海表温度和modis、SeaWiFS反演得到的海色图像得到海表面流场。Bowen做了一项很复杂的MCC方法反演研究，他们选取澳大利亚东部海域七年的AVHRR数据，用MCC法估算海表面流的精度达到了0.08到0.2m/s。Tokmakian选取California附近海域，用ADCP测量的结果与MCC得到的结果进行的比较，均方差在0.14-0.25m/s。

但这些方法都存在以下问题：

无论是以SST，还是水色要素浓度为示踪物，都只能得到晴空区结果；

由于方法本身的局限性，即使晴空区也会出现无法反演的现象；

因为这些严格的条件，限制了MCC技术的广泛应用。

发明内容

为了解决现有技术存在的上述问题，本发明提供了一种MCC表面海流反演方法。本发明针对如何解决MCC方法获得的表面海流在云覆盖区或由于方法本身造成缺测的问题提出了一种优化方法，从而提高反演覆盖度及精度。

本发明所采用的技术方案为：

一种MCC表面海流反演方法，其改进之处在于：所述方法包括

对MODIS卫星进行预处理；

利用综合优化云检测算法，进行云覆盖检测；

对于晴空区，分别选取海表温度和550nm反射率数据为示踪物，采用MMC反演得到表面海流场；

对于缺测的海表流数据，利用DINEOF方法进行插值补缺；

进行变分同化优化调整获得短时平均表面流。

进一步地，所述对MODIS卫星进行预处理包括

(1.1)进行辐射值计算：将卫星观测到所储存的计数值转换为可用的物理值；

其中，辐亮度L和反射率数据定标公式如下：

R_B,T,FS＝reflectance_scales_B(SI_B,T,FS-reflectance_offsets_B)；

式中，R_B,T,FS为反射通道的反射率值，SI_B,T,FS是数据集中16位计数值，reflectance_scales_B和reflectance_offsets_B是对应波段的缩放比和截距参数；

L_B,T,FS＝radiance_scales_B(SI_B,T,FS-radiance_offsets_B)；

式中，L_B,T,FS为辐射通道的辐亮度值，SI_B,T,FS是数据集中16位计数值，radiance_scales_B和radiance_offsets_B是对应波段的缩放比和截距参数；

由普朗克公式：

$L=\frac{2{\mathrm{hc}}^2}{{\mathrm{\λ}}^5(e^{\frac{hc}{\mathrm{\λ}kBT}}-1)};$

可以将辐亮度L_B.T.F.S转换为亮温BT(Brightness Temperature)：

$BT=\frac{hc}{\mathrm{\λ}kln(\frac{2{\mathrm{hc}}^2}{{\mathrm{L\λ}}^5}+1)};$

其中h＝6.626×10^-34J·s为普朗克常数，c＝2.9979×10⁸m·s^-1为真空光速，k＝1.38×10^-23J·K^-1为玻尔兹曼常数，λ代表波段中心波长；

(1.2)进行条带噪声去除：定位噪声带行数据，对其数据进行滤波平滑处理，对噪声带采取空间域的超限邻域平均-选择式掩膜平滑的方法来减弱噪声带的影响；

(1.3)进行边缘数据重叠的纠正：采用Kriging插值方法和分形插值算法对MODIS资料进行几何纠正；

(1.4)进行太阳天顶角纠正：将不同时刻太阳天顶角下的探测数据换算成相当于太阳处于天顶时的探测值；

其订正公式为：

R_c＝R_B.T.F.S/cos(θ)；

式中，R_c为订正后的反射率，R_B.T.F.S为定标后的反射率，θ为太阳天顶角；

(1.5)进行临边效应的订正：将各探测值换算成相当于卫星在天顶探测时的测值；

其订正方法为：

T＝T_B＋ΔT；

其中，T：订正后的温度；T_B：Planck亮度温度；

$T_B=C_2\mathrm{\υ}/In(\frac{c_1{\mathrm{\υ}}^3}{E}+1);$

式中，υ：探测波段的中心波数；E：定标后的辐射率；C₁,C₂为常数，C₁＝1.1910659·10^-5,C₂＝1.438833；ΔT：温度订正值；

$\mathrm{\Δ}T=(e^{0.00012{\mathrm{\θ}}^2}-1)\·(001072\·T_B-26.81);$

其中，θ：探测点的卫星天顶角；

$\mathrm{\θ}={\sin }^{-1}\left(\frac{\mathrm{Re}+H}{\mathrm{Re}}sin\mathrm{\σ}\right);$

其中，Re：地球半径；H：卫星高度；σ：探测点的卫星探测角；

(1.6)进行图像的投影：采用的投影方式为等经纬度投影。

进一步地，所述利用综合优化云检测算法，进行云覆盖检测包括

(2.1)进行时间识别：白天云检测所限定的太阳高度角要小于85°，大于85°视为夜间；

(2.2)进行海陆识别：区分出海洋、陆地与内陆湖泊，通过获取MOD03里的海陆标识数据集Land/SeaMask实现海陆识别；

(2.3)进行雪识别：采用归一化雪盖指数测量雪对可见光和短波红外波段的反射特性和反射差的相对大小；

(2.4)判别洋面太阳耀斑区；

当太阳反射角度θ_r位于0°～36°时，作太阳耀斑区处理；太阳反射角按下式确定：

其中，θ为太阳天顶角，θ₀为视角，为方位角；

(2.5)针对不同的下垫面，进行单一像元的分组检测，将检测项目分为高云、中云、低云三组检测。

进一步地，所述利用最大相关系数方法反演得到表面海流场包括

在特征量跟踪法中用于估计海流的第一幅图设为模板窗口，第二幅图设为搜索窗口，利用最大相关系数法来跟踪一幅图到下一幅图的海表特征的位移，在搜索窗口中找出与模板窗口大小相同且相关性的子窗口，对模板窗口进行重复的操作，得到完整的表面海流区域；

其中，示踪物选取：以海表面温度和叶绿素浓度来作为估计海流的示踪物，所述示踪物包括荧光素、叶绿素、悬浮泥沙、表面温度及盐度。

进一步地，所述对于缺测的海表流数据，利用DINEOF方法进行插值补缺包括通过经验正交函数来推测出缺值点的数据，其基本原理包括：

(5.1)构造观测数据矩阵；

(X₀)_m×n原始观测矩阵为m×n的数据矩阵，m为数据的空间维数，n为数据时间上的维数，缺值点用0表示；(X₀)_m×n矩阵的每一行表示当前时刻流场在某一方向上的流矢量，流矢量个数为m；每一列代表有n个连续相邻时次的流场结果用于流场重构；

(5.2)对X₀进行SVD分解；

选择EOF的最优数n，X₀的SVD分解可表示如下X₀＝UDV^T，U的维数为m×r，表示X₀的空间分解量，V的维数为r×n，表示X₀的时间分解量，D为对角特征值矩阵，维数为r×r，在得到U,D,V的值后，可按照X₀的重构算法(X_a)_ij＝(U_ND_NV_N ^T)_ij，来得到；

(5.3)缺值点数值由新的值代替；

X_a＝UDV^T，(X_a)_ij＝(U_ND_NV_N ^T)_ij直到收敛为止；

(5.4)EOF的最优选择数n的确定；

通过交叉验证方法选择n，采用MC客观分析的交叉验证方法以决定最优EOF阶数。

进一步地，所述进行优化调整包括

采用变分同化方法并结合正则化思想对表面海流场全局优化：

设特征量跟踪法得到的表面流场为(u^b,v^b)，经过变分同化的质量控制得到流场(u,v)；

寻找表面流场的分析值(u,v)使得(u,v)满足散度场：

$\frac{\∂u}{\∂x}+\frac{\∂v}{\∂y}=D(x,y)---\left(1\right)$

(u^b,v^b)为推导得到的表面流场，使得分析值(u,v)与表面流场充分接近，满足式(1)时，同时满足下式：

$J\[u,v\]=\underset{\mathrm{\Ω}}{\∫}\[{(u-u^b)}^2+{(v-v^b)}^2\]dxdy=min!---\left(2\right)$

利用Lagrange乘子法，把式(1)作为约束条件，即条件变分法，引进新的泛函如下：

$J_1\[u,v\]=\underset{\mathrm{\Ω}}{\∫}\[{(u-u^b)}^2+{(v-v^b)}^2-2\mathrm{\λ}(x,y)(\frac{\∂u}{\∂x}+\frac{\∂v}{\∂y}-D(x,y\left)\right)\]dxdy=min---\left(3\right)$

利用欧拉方程得：

$0=\frac{{\mathrm{\δJ}}_1}{2}=\underset{\mathrm{\Ω}}{\∫}\[(u-u^b)\mathrm{\δ}u+(v-v^b)\mathrm{\δ}v-\mathrm{\λ}(x,y)(\frac{\∂\left(\mathrm{\δ}u\right)}{\∂x}+\frac{\∂\left(\mathrm{\δ}v\right)}{\∂y})\]dxdy---\left(4\right)$

$\underset{\mathrm{\Ω}}{\∫}\[(u-u^b+\frac{\∂\mathrm{\λ}}{\∂x})\mathrm{\δ}u+(\[(v-v^b+\frac{\∂\mathrm{\λ}}{\∂y})\mathrm{\δ}v\]dxdy-\underset{\∂\mathrm{\Ω}}{\∫}\mathrm{\λ}(\mathrm{\δ}u,\mathrm{\δ}v)nds=0---\left(5\right)$

得出：

$\left\{\begin{array}{c}u-u^b+\frac{\∂\mathrm{\λ}}{\∂x}=0\\ v-v^b+\frac{\∂\mathrm{\λ}}{\∂y}=0\end{array}\right.---\left(6\right)$

和边界条件λ|_Γ＝0，即：

$\left\{\begin{array}{c}u=u^b-\frac{\∂\mathrm{\λ}}{\∂x}\\ v=v^b-\frac{\∂\mathrm{\λ}}{\∂y}\end{array}\right.---\left(7\right)$

另一方面：(u,v)满足：

$\frac{\∂u}{\∂x}+\frac{\∂v}{\∂y}=\frac{\∂u^b}{\∂x}+\frac{\∂v^b}{\∂y}=D(x,y)---\left(8\right)$

将(7)代入(8)，

$\frac{\∂u^b}{\∂x}+\frac{\∂v^b}{\∂y}-\mathrm{\Δ}\mathrm{\λ}=D(x,y)---\left(9\right)$

得出：

$\left\{\begin{array}{c}\mathrm{\Δ}\mathrm{\λ}=D^b(x,y)-D(x,y)\\ \mathrm{\λ}{|}_{\mathrm{\Γ}}=0\end{array}\right.---\left(10\right)$

解出λ，即可得到(u,v)；

利用松驰法得到λ，从而得到分析值u、v；

令f(x,y)＝D^b(x,y)-D(x,y)，将(9)式离散得：

λ_i+1,j＋λ_i-1,j+λ_i,j＋1＋λ_i,j-1-4λ_i,j＝f_i,j (11)

其中，每一格点上给定λ的初始猜值它不满足方程(10)，经过m次叠代后，余量重新求出(i,j)点值，使得余量等于或接近于0，从而得到λ_i,j的最后叠代值；将得到的λ(x,y)代入(7)式，就可得到表面海流的分析值(u,v)。

本发明的有益效果为：

(1)提出的云检测方法，提高了海面低云、碎云以及部分被云覆盖像元的识别准确率，为下一步计算表面海流更精确的去除了干扰，保证了后续结果的精度。

(2)在插值基础上，提出的一种表面海流场全局优化调整方法，大大增加了反演结果的覆盖范围，同时也使得反演精度得到明显提高，增强了该方法的使用范围。

本发明主要从工程的角度，依据长期的实践经验及技巧，提出了一种适合业务运行的表面海流反演算法。

附图说明

图1是本发明提供的一种MCC表面海流反演方法流程示意图；

图2是本发明提供的一种MCC表面海流反演方法中云检测流程示意图；

图3是本发明提供的一种MCC表面海流反演方法中最大相关系数法计算海流实施例示意图；

图4是本发明提供的一种MCC表面海流反演方法中松驰迭代示意图。

具体实施方式

下面结合附图对本发明的具体实施方式作进一步的详细说明。

以下描述和附图充分地示出本发明的具体实施方案，以使本领域的技术人员能够实践它们。其他实施方案可以包括结构的、逻辑的、电气的、过程的以及其他的改变。实施例仅代表可能的变化。除非明确要求，否则单独的组件和功能是可选的，并且操作的顺序可以变化。一些实施方案的部分和特征可以被包括在或替换其他实施方案的部分和特征。本发明的实施方案的范围包括权利要求书的整个范围，以及权利要求书的所有可获得的等同物。在本文中，本发明的这些实施方案可以被单独地或总地用术语“发明”来表示，这仅仅是为了方便，并且如果事实上公开了超过一个的发明，不是要自动地限制该应用的范围为任何单个发明或发明构思。

如图1所示，本发明提供了一种MCC表面海流反演方法。

本发明首先对EOS/MODIS卫星进行投影变换、几何校正、辐射变换，利用综合优化云检测算法，进行云覆盖检测。对于晴空区，分别选取海表温度和550nm反射率数据为示踪物，利用最大相关系数方法反演得到表面海流场。对于缺测的海表流数据，利用DINEOF方法进行插值补缺，利用变分同化方法并结合正则化思想的表面海流场全局优化方法进行优化调整，调整后的数据在覆盖度及反演精度上都有较大提高。具体步骤包括：

第一步：对MODIS卫星进行预处理；

(1)辐射计算

MODIS数据在使用前，必须经过辐射值计算，将卫星观测到所储存的计数值转换为可用的物理值。

其中，辐亮度L和反射率数据定标公式如下：

R_B,T,FS＝reflectance_scales_B(SI_B,T,FS-reflectance_offsets_B)；

式中，R_B,T,FS为反射通道的反射率值，SI_B,T,FS是数据集中16位计数值，reflectance_scales_B和reflectance_offsets_B是对应波段的缩放比和截距参数；

L_B,T,FS＝radiance_scales_B(SI_B,T,FS-radiance_offsets_B)；

式中，L_B,T,FS为辐射通道的辐亮度值，SI_B,T,FS是数据集中16位计数值，radiance_scales_B和radiance_offsets_B是对应波段的缩放比和截距参数；

由普朗克公式：

$L=\frac{2{\mathrm{hc}}^2}{{\mathrm{\λ}}^5(e^{\frac{hc}{\mathrm{\λ}kBT}}-1)};$

可以将辐亮度L_B.T.F.S转换为亮温BT(Brightness Temperature)：

$BT=\frac{hc}{\mathrm{\λ}kln(\frac{2{\mathrm{hc}}^2}{{\mathrm{L\λ}}^5}+1)};$

其中h＝6.626×10^-34J·s为普朗克常数，c＝2.9979×10⁸m·s^-1为真空光速，k＝1.38×10^-23J·K^-1为玻尔兹曼常数，λ代表波段中心波长。

(2)条带噪声去除

由于MODIS的条带噪声非常有规律的，条带噪声线平行于扫描方向，而且相邻条带噪声线之间的距离为扫描宽度，但是产生条带噪声线的位置并不总是在扫描带的分界处。所以可以先定位噪声带行数据，在对其进行滤波平滑。对噪声带采取空间域的超限邻域平均-选择式掩膜平滑的方法来减弱噪声带的影响；

(3)边缘数据重叠的纠正

由于像素的实际地面大小随着观测角的增大而增大，除了星下点数据外，相邻两行的像素存在彼此的重叠本发明采用Kriging插值方法和分形插值技术对MODIS资料进行了几何纠正。

(4)太阳天顶角纠正

把不同时刻的太阳天顶角下的探测数据换算成相当于太阳处于天顶时的探测值。

其订正公式为：

R_c＝R_B.T.F.S/cos(θ)

式中，R_c为订正后的反射率，R_B.T.F.S为定标后的反射率，θ为太阳天顶角。

(5)临边效应的订正

将各探测值换算成相当于卫星在天顶探测时的测值。

其订正方法为：

T＝T_B＋ΔT；

其中，T：订正后的温度；T_B：Planck亮度温度；

$T_B=C_2\mathrm{\υ}/In(\frac{c_1{\mathrm{\υ}}^3}{E}+1);$

式中，υ：探测波段的中心波数；E：定标后的辐射率；C₁,C₂为常数，C₁＝1.1910659·10^-5,C₂＝1.438833；ΔT：温度订正值；

$\mathrm{\Δ}T=(e^{0.00012{\mathrm{\θ}}^2}-1)\·(001072\·T_B-26.81);$

其中，θ：探测点的卫星天顶角；

$\mathrm{\θ}={\sin }^{-1}\left(\frac{\mathrm{Re}+H}{\mathrm{Re}}sin\mathrm{\σ}\right);$

其中，Re：地球半径；H：卫星高度；σ：探测点的卫星探测角。

(6)图像的投影

为便于计算，本发明采用的投影方式为等经纬度投影。

第二步：利用综合优化云检测算法，进行云覆盖检测；

本发明云检测方法在借鉴了NASA云检测可信度思想基础上，根据实际的资料获取情况来实现对研究地区的云检测。如图2所示，具体检测步骤如下：

(1)时间识别：白天云检测所限定的太阳高度角要小于85°，大于85°视为夜间。通过时间识别可以选择适当的云检测方法。

(2)海陆识别：区分出海洋、陆地与内陆湖泊。通过获取MOD03里的海陆标识数据集Land/Sea Mask实现海陆识别。

(3)雪识别：由于无法实时获取每日的冰雪覆盖数据，但为了提高云检测的可信度，在进行云检测之前，先进行积雪检测。目前基于反射特性的归一化雪盖指数(NDSI)具有普遍的实际操作意义，精度高，分类合理，是提取积雪信息的最佳技术手段。归一化雪盖指数(NDSI)是基于雪对可见光和短波红外波段的反射特性和反射差的相对大小的一种测量方法。

(4)洋面太阳耀斑区判别。当太阳反射角度θ_r位于0°～36°时，作太阳耀斑区处理；太阳反射角按下式确定：

其中，θ为太阳天顶角，θ₀为视角，为方位角。

(5)针对不同的下垫面，进行单一像元的分组检测。将检测项目分为高云、中云、低云三组检测：

第一组为高云检测：BT13.9，BT6.7，R1.38，BT11，BT3.7-BT12。

第二组为中云检测：R0.66，R0.87，R0.87/R0.66，BT11-BT7.3。

第三组为低云检测：BT11-BT3.9，R0.936/R0.865，B T11。

第三步：对于晴空区，分别选取海表温度和550nm反射率数据为示踪物，采用MMC反演得到表面海流场；

(1)示踪物选取

卫星遥感具有大范围系统反复的观测能力，可以被用来开发观测海洋表面环流，这种表面环流的估计方法一般都要观测示踪物及它们在一段时间内的移动，荧光素、叶绿素、悬浮泥沙、表面温度、及盐度都是潜在的示踪物，本方法主要以海表面温度和叶绿素浓度来作为估计海流的示踪物。

(2)示踪物追踪方法

如图3所示，在特征量跟踪法中用于估计海流的第一幅图称为“模板窗口(template window)”，第二幅图称为“搜索窗口(search window)”。利用最大相关系数法来跟踪一幅图到下一幅图的海表特征的位移。在“搜索窗口”中找出与“模板窗口”大小相同且相关性最好的那个子窗口，假设“搜索窗口”子窗口的中心坐标为(Xmcc，Ymcc)，“模板窗口”的中心坐标为(X，Y)，则有(Xmcc-X，Ymcc-Y)既可假设为流动的矢量。对许多“模板窗口”进行重复的操作，将得到完整的表面海流区域。

第四步：对于缺测的海表流数据，利用DINEOF方法进行插值补缺；

本方法采用DINEOF(Data Interpolating Empirical Orthogonal Functions)方法对缺值流场进行补缺。

DINEOF是一种自适应的EOF分解方法，通过经验正交函数来推测出缺值点的数据，而不需要通过计算最小的误差协方差来求取插值的结果，其基本原理如下：

(1)构造观测数据矩阵；

(X₀)_m×n原始观测矩阵为m×n的数据矩阵，m为数据的空间维数，n为数据时间上的维数，缺值点用0表示；(X₀)_m×n矩阵的每一行表示当前时刻流场在某一方向上的流矢量，流矢量个数为m；每一列代表有n个连续相邻时次的流场结果用于流场重构；

(2)对X₀进行SVD分解；

选择EOF的最优数n，第(4)步中说明选择原则。X₀的SVD分解可表示如下X₀＝UDV^T，U的维数为m×r，表示X₀的空间分解量，V的维数为r×n，表示X₀的时间分解量。D为对角特征值矩阵，维数为r×r。在得到U,D,V的值后，可按照X₀的重构算法(X_a)_ij＝(U_ND_NV_N ^T)_ij，来得到。

(3)缺值点数值由新的值代替；

X_a＝UDV^T，(X_a)_ij＝(U_ND_NV_N ^T)_ij直到收敛为止。

(4)EOF的最优选择数n的确定。

n的选择依赖于Mate Caro方法，亦称为交叉检验技术。如果n选取太小，重构的数据不足以反映数据内在的物理变化特征；n选太大，所引入的误差信息加大，计算量也会加大。采用MC客观分析的交叉验证方法以决定最优EOF阶数。

第五步：进行变分同化优化调整获得短时平均表面流；

变分同化实际上是偏微分方程的最优控制问题，涉及到偏微分方程理论、线性与非线性泛函分析、资料最优处理及数值计算等学科。

本专利提出的变分同化方法并结合正则化思想的表面海流场全局优化具体如下：

设特征量跟踪法得到的表面流场为(u^b,v^b)，我们要得到的是经过变分同化的质量控制得到流场(u,v)；

寻找表面流场的分析值(u,v)使得(u,v)满足散度场：

$\frac{\∂u}{\∂x}+\frac{\∂v}{\∂y}=D(x,y)---\left(1\right)$

(u^b,v^b)为推导得到的表面流场，使得分析值(u,v)与表面流场充分接近，满足式(1)时，同时满足下式：

$J\[u,v\]=\underset{\mathrm{\Ω}}{\∫}\[{(u-u^b)}^2+{(v-v^b)}^2\]dxdy=min!---\left(2\right)$

利用Lagrange乘子法，把式(1)作为约束条件，即条件变分法，引进新的泛函如下：

$J_1\[u,v\]=\underset{\mathrm{\Ω}}{\∫}\[{(u-u^b)}^2+{(v-v^b)}^2-2\mathrm{\λ}(x,y)(\frac{\∂u}{\∂x}+\frac{\∂v}{\∂y}-D(x,y\left)\right)\]dxdy=min---\left(3\right)$

利用欧拉方程得：

$0=\frac{{\mathrm{\δJ}}_1}{2}=\underset{\mathrm{\Ω}}{\∫}\[(u-u^b)\mathrm{\δ}u+(v-v^b)\mathrm{\δ}v-\mathrm{\λ}(x,y)(\frac{\∂\left(\mathrm{\δ}u\right)}{\∂x}+\frac{\∂\left(\mathrm{\δ}v\right)}{\∂y})\]dxdy---\left(4\right)$

$\underset{\mathrm{\Ω}}{\∫}\[(u-u^b+\frac{\∂\mathrm{\λ}}{\∂x})\mathrm{\δ}u+(\[(v-v^b+\frac{\∂\mathrm{\λ}}{\∂y})\mathrm{\δ}v\]dxdy-\underset{\∂\mathrm{\Ω}}{\∫}\mathrm{\λ}(\mathrm{\δ}u,\mathrm{\δ}v)nds=0---\left(5\right)$

得出：

$\left\{\begin{array}{c}u-u^b+\frac{\∂\mathrm{\λ}}{\∂x}=0\\ v-v^b+\frac{\∂\mathrm{\λ}}{\∂y}=0\end{array}\right.---\left(6\right)$

和边界条件λ|_Γ＝0，即：

$\left\{\begin{array}{c}u=u^b-\frac{\∂\mathrm{\λ}}{\∂x}\\ v=v^b-\frac{\∂\mathrm{\λ}}{\∂y}\end{array}\right.---\left(7\right)$

另一方面：(u,v)满足：

$\frac{\∂u}{\∂x}+\frac{\∂v}{\∂y}=\frac{\∂u^b}{\∂x}+\frac{\∂v^b}{\∂y}=D(x,y)---\left(8\right)$

将(7)代入(8)，

$\frac{\∂u^b}{\∂x}+\frac{\∂v^b}{\∂y}-\mathrm{\Δ}\mathrm{\λ}=D(x,y)---\left(9\right)$

得出：

$\left\{\begin{array}{c}\mathrm{\Δ}\mathrm{\λ}=D^b(x,y)-D(x,y)\\ \mathrm{\λ}{|}_{\mathrm{\Γ}}=0\end{array}\right.---\left(10\right)$

解出λ，即可得到(u,v)。

如图4所示，利用松驰法(Relaxation method)得到λ，从而得到分析值u、v；

令f(x,y)＝D^b(x,y)-D(x,y)，将(9)式离散得：

λ_i+1,j+λ_i-1,j+λ_i,j+1+λ_i,j-1-4λ_i,j＝f_i,j (11)

每一格点上给定λ的初始猜值一般来说，它不满足方程(10)，经过m次叠代后，余量重新求出(i,j)点值，使得余量等于或接近于0，从而得到λ_i,j的最后叠代值。将得到的λ(x,y)代入(7)式，就可以得到表面海流的分析值(u,v)。

本发明不局限于上述最佳实施方式，任何人在本发明的启示下都可得出其他各种形式的产品，但不论在其形状或结构上作任何变化，凡是具有与本申请相同或相近似的技术方案，均落在本发明的保护范围之内。

Claims (6)

1.一种MCC表面海流反演方法，其特征在于：所述方法包括

对MODIS卫星进行预处理；

利用综合优化云检测算法，进行云覆盖检测；

对于晴空区，分别选取海表温度和550nm反射率数据为示踪物，采用MMC反演得到表面海流场；

对于缺测的海表流数据，利用DINEOF方法进行插值补缺；

进行变分同化优化调整获得短时平均表面流。

2.根据权利要求1所述的一种MCC表面海流反演方法，其特征在于：所述对MODIS卫星进行预处理包括

(1.1)进行辐射值计算：将卫星观测到所储存的计数值转换为可用的物理值；

其中，辐亮度L和反射率数据定标公式如下：

R_B,T,FS＝reflectance_scales_B(SI_B,T,FS-reflectance_offsets_B)；

式中，R_B,T,FS为反射通道的反射率值，SI_B,T,FS是数据集中16位计数值，reflectance_scales_B和reflectance_offsets_B是对应波段的缩放比和截距参数；

L_B,T,FS＝radiance_scales_B(SI_B,T,FS-radiance_offsets_B)；

式中，L_B,T,FS为辐射通道的辐亮度值，SI_B,T,FS是数据集中16位计数值，radiance_scales_B和radiance_offsets_B是对应波段的缩放比和截距参数；

由普朗克公式：

$L=\frac{2{\mathrm{hc}}^2}{{\mathrm{\λ}}^5(e^{\frac{hc}{\mathrm{\λ}kBT}}-1)};$

可以将辐亮度L_B.T.F.S转换为亮温BT：

$BT=\frac{hc}{\mathrm{\λ}kln(\frac{2{\mathrm{hc}}^2}{{\mathrm{L\λ}}^5}+1)};$

其中，h＝6.626×10^-34J·s为普朗克常数，c＝2.9979×10⁸m·s^-1为真空光速，k＝1.38×10^-23J·K^-1为玻尔兹曼常数，λ代表波段中心波长；

(1.2)进行条带噪声去除：定位噪声带行数据，对其数据进行滤波平滑处理，对噪声带采取空间域的超限邻域平均-选择式掩膜平滑的方法来减弱噪声带的影响；

(1.3)进行边缘数据重叠的纠正：采用Kriging插值方法和分形插值算法对MODIS资料进行几何纠正；

(1.4)进行太阳天顶角纠正：将不同时刻太阳天顶角下的探测数据换算成相当于太阳处于天顶时的探测值；

其订正公式为：

R_c＝R_B.T.F.S/cos(θ)；

式中，R_c为订正后的反射率，R_B.T.F.S为定标后的反射率，θ为太阳天顶角；

(1.5)进行临边效应的订正：将各探测值换算成相当于卫星在天顶探测时的测值；

其订正方法为：

T＝T_B+ΔT；

其中，T：订正后的温度；T_B：Planck亮度温度；

$T_B=C_2\mathrm{\υ}/In(\frac{c_1{\mathrm{\υ}}^3}{E}+1);$

式中，υ：探测波段的中心波数；E：定标后的辐射率；C₁,C₂为常数，C₁＝1.1910659·10^-5,C₂＝1.438833；ΔT：温度订正值；

$\mathrm{\Δ}T=(e^{0.00012{\mathrm{\θ}}^2}-1)\·(001072\·T_B-26.81);$

其中，θ：探测点的卫星天顶角；

$\mathrm{\θ}={\sin }^{-1}\left(\frac{\mathrm{Re}+H}{\mathrm{Re}}sin\mathrm{\σ}\right);$

其中，Re：地球半径；H：卫星高度；σ：探测点的卫星探测角；

(1.6)进行图像的投影：采用的投影方式为等经纬度投影。

3.根据权利要求1所述的一种MCC表面海流反演方法，其特征在于：所述利用综合优化云检测算法，进行云覆盖检测包括

(2.1)进行时间识别：白天云检测所限定的太阳高度角要小于85°，大于85°视为夜间；

(2.2)进行海陆识别：区分出海洋、陆地与内陆湖泊，通过获取MOD03里的海陆标识数据集Land/Sea Mask实现海陆识别；

(2.3)进行雪识别：采用归一化雪盖指数测量雪对可见光和短波红外波段的反射特性和反射差的相对大小；

(2.4)判别洋面太阳耀斑区；

当太阳反射角度θ_r位于0°～36°时，作太阳耀斑区处理；太阳反射角按下式确定：

其中，θ为太阳天顶角，θ₀为视角，为方位角；

(2.5)针对不同的下垫面，进行单一像元的分组检测，将检测项目分为高云、中云、低云三组检测。

4.根据权利要求1所述的一种MCC表面海流反演方法，其特征在于：所述利用最大相关系数方法反演得到表面海流场包括

在特征量跟踪法中用于估计海流的第一幅图设为模板窗口，第二幅图设为搜索窗口，利用最大相关系数法来跟踪一幅图到下一幅图的海表特征的位移，在搜索窗口中找出与模板窗口大小相同且相关性的子窗口，对模板窗口进行重复的操作，得到完整的表面海流区域；

其中，示踪物选取：以海表面温度和叶绿素浓度来作为估计海流的示踪物，所述示踪物包括荧光素、叶绿素、悬浮泥沙、表面温度及盐度。

5.根据权利要求1所述的一种MCC表面海流反演方法，其特征在于：所述对于缺测的海表流数据，利用DINEOF方法进行插值补缺包括通过经验正交函数来推测出缺值点的数据，其基本原理包括：

(5.1)构造观测数据矩阵；

(X₀)_m×n原始观测矩阵为m×n的数据矩阵，m为数据的空间维数，n为数据时间上的维数，缺值点用0表示；(X₀)_m×n矩阵的每一行表示当前时刻流场在某一方向上的流矢量，流矢量个数为m；每一列代表有n个连续相邻时次的流场结果用于流场重构；

(5.2)对X₀进行SVD分解；

选择EOF的最优数n，X₀的SVD分解可表示如下X₀＝UDV^T，U的维数为m×r，表示X₀的空间分解量，V的维数为r×n，表示X₀的时间分解量，D为对角特征值矩阵，维数为r×r，在得到U,D,V的值后，可按照X₀的重构算法(X_a)_ij＝(U_ND_NV_N ^T)_ij，来得到；

(5.3)缺值点数值由新的值代替；

X_a＝UDV^T，(X_a)_ij＝(U_ND_NV_N ^T)_ij直到收敛为止；

(5.4)EOF的最优选择数n的确定；

通过交叉验证方法选择n，采用MC客观分析的交叉验证方法以决定最优EOF阶数。

6.根据权利要求1所述的一种MCC表面海流反演方法，其特征在于：所述进行优化调整包括

采用变分同化方法并结合正则化思想对表面海流场全局优化：

设特征量跟踪法得到的表面流场为(u^b,v^b)，经过变分同化的质量控制得到流场(u,v)；

寻找表面流场的分析值(u,v)使得(u,v)满足散度场：

$\frac{\∂u}{\∂x}+\frac{\∂v}{\∂y}=D(x,y)---\left(1\right)$

(u^b,v^b)为推导得到的表面流场，使得分析值(u,v)与表面流场充分接近，满足式(1)时，同时满足下式：

$J\[u,v\]=\underset{\mathrm{\Ω}}{\∫}\[{(u-u^b)}^2+{(v-v^b)}^2\]dxdy=min!---\left(2\right)$

利用Lagrange乘子法，把式(1)作为约束条件，即条件变分法，引进新的泛函如下：

$J_1\[u,v\]=\underset{\mathrm{\Ω}}{\∫}\[{(u-u^b)}^2+{(v-v^b)}^2-2\mathrm{\λ}(x,y)(\frac{\∂u}{\∂x}+\frac{\∂v}{\∂y}-D(x,y\left)\right)\]dxdy=min---\left(3\right)$

利用欧拉方程得：

$0=\frac{{\mathrm{\δJ}}_1}{2}=\underset{\mathrm{\Ω}}{\∫}\[(u-u^b)\mathrm{\δ}u+(v-v^b)\mathrm{\δ}v-\mathrm{\λ}(x,y)(\frac{\∂\left(\mathrm{\δ}u\right)}{\∂x}+\frac{\∂\left(\mathrm{\δ}v\right)}{\∂y})\]dxdy---\left(4\right)$

$\underset{\mathrm{\Ω}}{\∫}\[(u-u^b+\frac{\∂\mathrm{\λ}}{\∂x})\mathrm{\δ}u+(\[(v-v^b+\frac{\∂\mathrm{\λ}}{\∂y})\mathrm{\δ}v\]dxdy-\underset{\∂\mathrm{\Ω}}{\∫}\mathrm{\λ}(\mathrm{\δ}u,\mathrm{\δ}v)nds=0---\left(5\right)$

得出：

$\left\{\begin{array}{c}u-u^b+\frac{\∂\mathrm{\λ}}{\∂x}=0\\ v-v^b+\frac{\∂\mathrm{\λ}}{\∂y}=0\end{array}\right.---\left(6\right)$

和边界条件λ|_Γ＝0，即：

$\left\{\begin{array}{c}u=u^b-\frac{\∂\mathrm{\λ}}{\∂x}\\ v=v^b-\frac{\∂\mathrm{\λ}}{\∂y}\end{array}\right.---\left(7\right)$

另一方面：(u,v)满足：

$\frac{\∂u}{\∂x}+\frac{\∂v}{\∂y}=\frac{\∂u^b}{\∂x}+\frac{\∂v^b}{\∂y}=D(x,y)---\left(8\right)$

将(7)代入(8)，

$\frac{\∂u^b}{\∂x}+\frac{\∂v^b}{\∂y}-\mathrm{\Δ}\mathrm{\λ}=D(x,y)---\left(9\right)$

得出：

$\left\{\begin{array}{c}\mathrm{\Δ}\mathrm{\λ}=D^b(x,y)-D(x,y)\\ \mathrm{\λ}{|}_{\mathrm{\Γ}}=0\end{array}\right.---\left(10\right)$

解出λ，即可得到(u,v)；

利用松驰法得到λ，从而得到分析值u、v；

令f(x,y)＝D^b(x,y)-D(x,y)，将(9)式离散得：

λ_i+1,j+λ_i-1,j+λ_i,j+1+λ_i,j-1-4λ_i,j＝f_i,j (11)

其中，每一格点上给定λ的初始猜值它不满足方程(10)，经过m次叠代后，余量重新求出(i,j)点值，使得余量等于或接近于0，从而得到λ_i,j的最后叠代值；将得到的λ(x,y)代入(7)式，就可得到表面海流的分析值(u,v)。