1 项目介绍

见[IOCAS 科创计划记录] 目录

2 项目进展记录

2.1 第一阶段：2019.8.5——2019.9.5

还在享受着快乐的暑假时光，突然接到导师在群里传的第一阶段任务：

第一阶段五人分工：
（1）历史温盐剖面数据搜集和质量控制：
搜集历史温盐数据（1980年-至今）并对数据进行质量控制。数据来源主要是https://www.nodc.noaa.gov/OC5/WOD/pr_wod.html。质量控制方法需查找相关文献。
（2）数据插值方法：查找已有的数据插值方法，选择最适合的插值方法，将（1）处理的温盐剖面数据插值成1/8º×1/8º空间分辨率1-12个月月平均网格化数据。
（3）地转流计算方法一：查找学习地转流计算方法，利用（2）插值后温盐数据计算地转流。
（4）地转流计算方法二：查找学习地转流计算方法，利用（2）插值后温盐数据计算地转流。
（5）温跃层的判定：查找相关判定方法，利用（2）插值后温盐数据计算逐月平均的温跃层分布。
参考文献：
Qiu, B., S. Chen, D. L. Rudnick, and Y. Kashino, 2015: A new paradigm for the North Pacific subthermocline low-latitude western boundary current system. Journal of Physical Oceanography, 45(9), 2407-2423.
给大家提供一篇参考文献，请同学们自己搜集更多的文献。

导师给了五个任务，我们一人一个，我选了(2) ，进行数据插值。
首先到WOD（World Ocean Database）网站上下载数据，我们研究的空间范围为$-10^{\circ}S-40^{\circ}N$−10∘S−40∘N，$115^{\circ}E-180^{\circ}E$115∘E−180∘E

使用WOD select选取所需时空范围内的CTD数据。数据可以下载为三种类型：ASCII格式、csv格式和netCDF格式。

我首先尝试了netCDF，发现可能是为了节省空间，温盐数据都按经纬度、深度被存为了一维数组，并且缺测的数据直接被略过，我也没找到缺测点的位置信息，导致无法与经纬度、深度一一对应，读取失败。
随后我就尝试了ASCII格式，并利用WOD提供的wodtodepthmatrix.exe将其转为矩阵，实现经纬度、深度、年月日、测量数据、qcflag的一一对应，至此原始数据（raw data）准备完毕。

下面进行了质量控制，即剔除不合理的数据。其实也就是根据数据自带的qcflag以及经纬度范围进行了筛选，其余的首先是我们担心将局部细节剔除以及受限于编程能力就没有进行。

现在轮到插值环节了。
再来回顾一下任务

（2）数据插值方法：查找已有的数据插值方法，选择最适合的插值方法，将（1）处理的温盐剖面数据插值成1/8º×1/8º空间分辨率1-12个月月平均网格化数据。

我首先查了一下文献，发现包括老师提供的文献,好几篇文献都提到克里金插值是比较好的方法^{1 ,2 ,3}。
$\left\{\begin{array}{l}{v_{e}=\frac{\sum_{i=1}^{n} w_{i} v_{i}}{\sum_{i=1}^{n} w_{i}}} \\ {w_{i}=e^{\left(\frac{\Delta x^{2}}{2 L_{x}^{2}}-\frac{\Delta y^{2}}{2 L_{y}^{2}}\right)}}\end{array}\right.$⎩⎨⎧​ve​=∑i=1n​wi​∑i=1n​wi​vi​​wi​=e(2Lx2​Δx2​−2Ly2​Δy2​)​
其中，v_e为所求插值格点上的数值v_i为散点数据的值,w_i为权重，并且按照文献[3]，L_x ,L_y都取0.5$^{\circ}$∘，插值后的图像结果看起来并不理想。

可以明显看出受反距离权重以及部分区域数据量不足的影响，会产生区域性斑块状的分布。用月平均的插值数据与WOA气候态数据（0.25$^{\circ}$∘×0.25$^{\circ}$∘）进行了对比，得到各月平均绝对值误差如下图所示。

误差不超过1$^{\circ}$∘C，但图像并不令人满意，在老师的指引下，我又看了WOA的插值方法^{4 ,5}。
该方法先将区域按纬度分为宽度为1$^{\circ}$∘的条带，将这个条带内格点上的数据赋值为该条带内所有数据的平均值作为初始数据场(first-guess field)，以1月海表为例，如下图所示。
然后对各个格点上的数据进行修正，修正方法为：
$\left\{\begin{array}{l}{C_{i, j}=\frac{\sum_{s=1}^{n} W_{s} Q_{s}}{\sum_{s}^{n} W_{s}}} \\ {W_{s}=e^{-\frac{E r^{2}}{R^{2}}}(r \leq R)} \\ {G_{i, j}=F_{i, j}+C_{i, j}}\end{array}\right.$⎩⎪⎨⎪⎧​Ci,j​=∑sn​Ws​∑s=1n​Ws​Qs​​Ws​=e−R2Er2​(r≤R)Gi,j​=Fi,j​+Ci,j​​
其中，i, j为待插值格点的坐标，C_i,j为该点修正值，n为影响半径内数据点的个数，Q_s为first-guess field与数据点的值之差，W_s为权重函数，r为数据点与格点间的距离，R为影响半径，这里取R²=100000m²，E=4，G_i,j为最终格点上的值。
我按照上述方法进行了插值，以1月海表为例，结果如下：

与克里金插值相比，使用WOA的插值方法得到的图像的平滑程度已经得到了比较大的提升，但仍然有部分区域有温度的斑状分布。为了解决该问题，我联想到了图像处理中的模糊处理方法。将温度看作像素，进行模糊处理。在众多滤波算子中，我选择了对细节保留较好的运动模糊算子，用自带的fspecial和imfilter函数对数据矩阵进行处理，以1月海表为例，结果如下：

可以明显看出，相较于克里金插值和WOA插值方法，模糊处理后的SST平滑程度有非常明显的提高，下图为模糊处理后得到的各月表层温度：

fspecial函数需要两个参数，我将与WOA的年平均误差作为目标函数对参数进行寻优：
$\begin{array}{ll}{\min } & {\text { error }=\sum_{i=1}^{n}\left|T_{i}-T_{i}^{*}(a, b)\right|} \\ {\text { s.t. }} & {a>0, b>0}\end{array}$min s.t. ​ error =∑i=1n​∣Ti​−Ti∗​(a,b)∣a>0,b>0​
用MATLAB的fmincon函数经多次赋初始值得到局部最优解$a=65.4019,b=485.9608$a=65.4019,b=485.9608
使用运动模糊处理的方法，不仅能让平滑程度显著提高，误差也明显降低了，月平均及年平均误差分别如下图所示：


至此，我的任务告一段落，将数据了交给下面要计算地转流和温跃层的同学。

标签：IOCAS,2019.8,2019.9,转流,插值,温盐,方法,circ,数据
来源： https://blog.csdn.net/cyzzym000/article/details/100887279