杭州湾三维悬浮泥沙输运模型初始场 的伴随法反演研究*

杜云飞 张继才 王道胜

(1. 浙江大学海洋学院 舟山 316000； 2. 中国地质大学(武汉)海洋学院 海洋地质资源湖北省重点实验室 武汉 430074； 3. 南方海洋科学与工程广东省实验室(广州) 广州 511458)

由于径流输入和局部泥沙再悬浮， 近岸、河口海域通常具有高悬浮泥沙浓度的特征(Yanget al，2016)。杭州湾作为典型的强潮河口湾， 具有潮流急， 潮差大， 水体含沙量高的显著特征。大量悬浮泥沙的存在和大范围含沙量时间、空间变化的存在， 使杭州湾悬浮泥沙分布比较复杂(中国海湾志编纂委员会， 1992)。因此， 准确了解杭州湾海域悬浮泥沙的时空分布和输运规律， 不仅具有重要的科学价值， 还对杭州湾海域及近岸地区水质、地貌、生态环境、海岸工程、港口建设以及滩涂开发等具有重要意义。数值模拟具有成本低、预报性强、可移植性强、定量的特点， 可有效获取悬浮泥沙的时间和空间分布(Wanget al，2018b)， 从而弥补传统现场观测和卫星遥感观测的不足。因此， 悬沙输运模型已经成为研究泥沙输运的一个有力工具(Papanicolaouet al， 2008)。

随着对泥沙输运机制认知的深入以及科学计算能力的不断进步， 悬沙输运模型从一维发展到了三维。众多学者利用三维悬沙输运模型开展悬浮泥沙的研究工作(Normant， 2000； Dufoiset al， 2014； Chenet al， 2015)。三维悬沙输运模型对初始场比较敏感(Wanget al，2018a)， 即 悬 浮 泥 沙 浓 度(suspended sediment concentration， SSC)在数值模拟初始时刻的空间分布直接影响数值模拟的精度， 初始场的空间分布与实际SSC 的空间分布越接近， 相应的数值模拟精度越高， 反之则越低。因此， 对于短期的悬浮泥沙数值模拟而言， 初始场的准确给定至关重要(Leeet al，2007； Yanget al， 2016)。到目前为止， 在悬浮泥沙输运的数值模拟中， 初始场的确定方法主要有四种： (1)假定在整个模型计算区域内为常数分布， 以模型运行稳定时的结果作为最终的模拟结果(Huet al，2009； Gourgueet al， 2013)； (2)根据现场实测水文资料插值得到(曹振轶等，2002； Duet al， 2010)； (3)由经验、半经验的公式确定(Woolnoughet al， 1995； Schwabet al， 2000)； (4)用卫星遥感观测资料获得(Chenet al， 2010； Ramakrishnanet al， 2012； Yanget al， 2016)。前三种方法与实际SSC 的空间分布偏差较大， 因而数值模拟的精度较差； 随着卫星遥感观测手段的不断发展， 初始场的确定更多依赖于卫星遥感数据及与之相关的反演方法和经验公式。卫星遥感反演方法虽然能够获取海洋表层大面积、连续的SSC， 但由于受到云覆盖或恶劣天气条件的限制， 卫星遥感数据存在不同程度的缺测， 这在一定程度上限制了卫星遥感方法的应用范围和准确度。

伴随同化方法是建立在严格的数学基础上的一种有效的四维变分同化技术， 它将所要解决的实际问题作为最小值问题来解决， 流体力学方程组及其初边值条件作为约束条件， 使得根据具体问题设计的代价函数达到最小(Sasaki， 1970)。伴随同化方法能够充分利用已有的观测资料， 实现观测资料和数值模型的有机结合， 得到更加符合实际的模拟结果， 并能对模型重要参数如初始场、边界条件等进行优化反演。目前为止， 伴随同化方法已广泛应用于海表温度数值模型(马继瑞等， 2002)、风暴潮预报模型(Penget al， 2006)、海洋生态系统动力学模型(Wanget al， 2013)、海洋环流模型(Songet al， 2016)以及悬沙输运模型(Wanget al， 2018b)等海洋数值模型中初始场的优化反演， 并取得了良好的应用效果。基于此背景， 本文利用Wang 等(2018a)构建的一个三维悬沙输运伴随同化模型， 针对模型初始场的优化反演进行一系列孪生实验和实际实验， 孪生实验验证伴随同化方法优化反演模型初始场的有效性并定量评估该模型数据同化的能力； 实际实验则根据孪生实验的结论同化杭州湾海域典型的小潮时期和大潮时期的GOCI卫星遥感资料所得表层SSC 数据， 获得模型初始场的最优估计。

1 模型介绍

1.1 正向模型

三维悬沙输运模型的控制方程(Wanget al，2018a)：

其中，C是悬浮泥沙浓度，t是时间，x、y是水平坐标，σ是垂向坐标，H是包含静水深和海面起伏的总水深，u、v和w分别是x、y和σ方向的流速，ws为悬浮泥沙的沉降速度，KH和KV分别为水平扩散系数和垂向扩散系数。

初始条件：

其中，C0是初始场。

边界条件：

其中，是垂直于边界的外法向量，Cobc是流入开边界处的悬浮泥沙浓度，B1代表陆地边界，B2代表海洋开边界，B3代表流入开边界，B4代表海面边界，B5代表底边界，E和D分别是侵蚀率和沉降率， 二者的计算公式如下(Partheniades， 1965)：

其中，M0是再悬浮率，τb是底部剪切应力值，τce是侵蚀的临界剪切应力值，τcd是沉降的临界剪切应力值，C1是近海底处的悬浮泥沙浓度值。

正向模型的离散形式和求解方式详见Wang 等(2018a)， 此处不再赘述。

1.2 伴随模型

根据拉格朗日乘子法， 引入伴随变量λ， 令拉格朗日函数关于各个变量和模型参数的梯度值为0， 可以推导出伴随模型以及代价函数关于初始场的梯度表达式， 具体的推导过程和差分格式详见Wang 等(2018a)， 此处不再赘述。

1.3 优化算法

由伴随模型求得代价函数关于初始场的梯度后， 可以利用优化算法通过迭代获得最优初始场， 迭代公式可以表示为：

本文中，αn是正值， 一般取为常数， 但调整步长过大或过小都会对反演结果产生影响。调整步长过大时， 在极小值附近会造成约化代价函数“过调整”的问题， 产生“锯齿效应”； 调整步长过小时， 会造成收敛速度慢， 计算效率低的问题。与以往调整步长取为常数不同， 本文中调整步长依赖于初始场的空间分布， 取为空间分布的形式， 即将αn取为初始场的第n个迭代步迭代结果的3%(公式13)， 这样能够保证迭代收敛的稳定性， 且不需要人为地凭经验给定调整步长， 可以避免由调整步长过大或者过小引起的上述问题， 同时避免需要反复调试才能获得最佳的调整步长， 提高了计算效率。

确定dn有许多可供选择的优化算法， 本文选择并比较了其中的三种， 分别是最速下降法(the gradient descent algorithm， GD)、共轭梯度法(the conjugate gradient algorithm， CG)和有限记忆 BFGS 法(the limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm， L-BFGS)。

在GD 中， 搜索方向为代价函数关于模型初始场的负梯度方向：

其中，dn为搜索方向，gn为代价函数关于模型初始场的梯度方向， 根据βn的不同求法， CG 又可以分为许多种类， 其中PRP、HS、LS 方法的数值表现较好， 这三种方法中βn的计算公式如下：

(1) PRP 方法(CG-PRP， Polaket al，1969； Polyak， 1969)：

在初始场优化反演的迭代过程中， 迭代收敛标准可以有很多， 当达到某一收敛标准时迭代终止。收敛标准可以是最后两个迭代步的代价函数值充分接近或小于某个值， 也可以是最后两个迭代步的代价函数关于初始场的梯度值充分接近或小于某值等。本文的孪生实验中， 迭代收敛标准是最后两个迭代步的约化代价函数值之差的绝对值小于3.0×10-5， 且最大迭代步设置为200， 即：

2 研究区域、模型设置及数据来源

本文的研究区域为杭州湾海域(图1)， 经纬度范围为29.902°—30.904°N、120.602°—121.806°E。背景流场数据由FVCOM(Chenet al，2003)计算得到， 与 Wang 等(2018b)相同。模型的水平分辨率为0.006°×0.006°； 垂向分为四层， 与输出背景流场数据的数值模型FVCOM 一致。正向模型的时间步长为180s， 总的模拟时间是31h， 小潮时期数值模拟实验的初始时刻是2011 年6 月26 日0：28， 结束时刻是2011 年6 月27 日7：28， 大潮时期数值模拟实验的初始时刻是2015 年3 月24 日0：28， 结束时刻是2015年3 月25 日7：28。河流流入开边界条件在小潮时期取为1.5×10-1kg/m3(中华人民共和国水利部， 2011)， 在大潮时期取为1.0×10-1kg/m3(中华人民共和国水利部， 2015)。根据唐建华(2007)的研究， 一个潮周期内杭州湾海域絮凝的平均沉速的量级为10-4m/s， 因此沉降速度取为1.0×10-4m/s。参照Hu 等(2009)在模拟杭州湾悬浮泥沙输运时再悬浮率的取值， 再悬浮率取为3.0×10-6m/s， 沉降和侵蚀的临界切应力均取为0.4N/m2。基于Wang 等(2018a)中的尺度分析， 水平扩散系数取为80m2/s， 垂向扩散系数在小潮时期取为3.0×10-3m2/s， 在大潮时期取为3.0×10-2m2/s。以上模型参数的具体描述和取值如表1 所示。

本文采用由静止水色卫星GOCI 的观测数据反演所得SSC 作为同化数据， 由He 等(2013)提供。GOCI 卫星能够对潮汐、赤潮、河流羽流、泥沙输运等短期的区域性海洋现象实现高时空分辨率的监测， 在水色遥感探测中具有很大的发展潜力(Heet al，2013)。GOCI 数据的空间分辨率为500m， 时间分辨率为 1h， 每天可以获得八幅观测图像(Ryuet al， 2011)。He 等(2013)构建了一个大气校正算法和经验关系模型， 可以反演得到小潮时期和大潮时期长江口和杭州湾海域表层SSC， 该数据在Wang 等(2018b)中有具体描述。

表1 主要模型参数 Tab.1 Main parameters of the model

数据同化的有效性需要用未被同化的独立观测进行检验(Elbernet al，2007)。因此， 在本文的孪生实验和实际实验中， 为了检验伴随同化方法优化反演模型初始场的有效性， 将随机选取10%的观测数据不参与同化而作为独立观测来检验数据同化的结果， 将这些观测数据称为“检验观测”(check observations， COs)； 其余90%的观测数据参与同化， 称其为“被同化观测”(assimilated observations， AOs)。孪生实验和实际实验中观测数据的唯一区别是， 孪生实验中的观测数据是给定模型参数值运行正向模型生成的表层“观测”， 实际实验中的观测数据则是实际GOCI 卫星遥感数据反演所得SSC。

3 孪生实验及结果分析

3.1 初始场的误差对模拟结果的影响

初始场的不确定性是模型误差的主要来源之一。基于三维悬沙输运伴随同化模型， 以小潮时期为例， 利用初始时刻实测的GOCI 卫星遥感观测数据经过水平空 间插值获得表层SSC 的初始场(图2a)， 插值方法采用MATLAB 软件自带的自然邻近插值法(natural)， 其余各层初始场由表层初始场线性插值得到， 插值公式与Wang 等(2018b)中相同， 该插值公式是由SSC 原位观测数据( 旸杨 等， 2008)线性拟合得到， 较好地体现了杭州湾SSC 由表及底逐渐增大的特点。运行正向模型获得整个研究区域的SSC“观测值”， 每小时输出一次结果。

为了评估初始场的误差对模拟结果的影响， 对初始场分别添加最大不超过10%、20%、30%的随机相对误差， 计算公式如下：

其中，C0+%为添加随机误差后的初始场，C0为给定初始场，p为最大随机相对误差百分比，r为介于-1—1之间的随机数。

图2 初始时刻实测的GOCI 卫星遥感观测数据经过水平空间插值后所得表层悬浮泥沙浓度的初始场 Fig.2 Initial field of the surface SSC obtained by horizontal spatial interpolation of GOCI satellite remote sensing observation data measured at initial time during neap tide

将垂向四层的SSC 模拟值与观测值之间的相对误差分别进行水平和垂向空间平均后得到整个研究区域的SSC 的平均相对误差(mean relative error， MRE)随模拟时间的变化情况(图3)。可以看出， SSC 模拟结果的MRE 均随初始场误差的增大而增大； 随着模拟时间的增加， 模拟结果的MRE 没有明显下降。这说明初始场对于三维悬沙输运模型是敏感参数， 初始场的给定越准确， 模拟结果的MRE 越小， 反之， 则越大。因此， 对于三维悬沙输运模型而言， 尤其是短期的SSC 数值模拟， 初始场的给定至关重要， 将直接影响SSC 的模拟精度。

3.2 不同优化算法对初始场优化反演的影响

图3 整个研究区域的SSC 模拟值与观测值之间的平均相对误差随模拟时间的变化情况 Fig.3 The change of mean relative errors of SSC between the simulation results and observations for the total study area with the model time increasing

图4 TE11—TE15 中反演后的表层初始场(a—e)、反演后的表层初始场与给定初始场的误差(f—j)、反演后的表层初始场与给定初始场的相对误差(k—o)的空间分布 Fig.4 The spatial distributions of the surface initial field after inversion in TE11—TE15 (a—e)； the spatial distributions of errors between the surface initial field after inversion and the given initial field in TE11—TE15 (f—j)； and the spatial distributions of relative errors between the surface initial field after inversion and the given initial field in TE11—TE15 (k—o)

图5 实验TE11—TE15 中约化代价函数(log10(J/J0)) (a)、代价函数关于初始场的梯度的L1 范数(Gradient of C0) (b)、SSC 模拟结果与“被同化观测”(AOs)之间的平均绝对误差(MAEs at AOs) (c)、SSC 模拟结果与“检验观测”(COs)之间的平均绝对误差(MAEs at COs) (d)、初始场的反演结果与给定分布之间的平均相对误差(MREs of C0) (e)随迭代步的变化 Fig.5 The values of log(J/J0) versus the iteration steps (a)； the variations of Gradient of C0 (b)； the variations of MAEs at AOs (c)； the variations of MAEs at COs (d)； and the variations of MREs of C0 (e) in TE11—TE15

表2 实验TE11—TE15 中， 数据同化前后的约化代价函数(log10(J/J0))、代价函数关于初始场的梯度的L1 范数(Gradient of C0)、SSC 模拟结果与AOs 之间的平均绝对误差(MAEs at AOs)、SSC 模拟结果与COs 之间的平均绝对误差(MAEs at COs)和初始场的反演结果与给定分布之间的平均相对误差(MREs of C0) Tab.2 log10(J/J0)， gradient of C0， MAEs at AOs， MAEs at Cos， and MREs of C0 before and after data assimilation in TE11—TE15

本节探讨 GD、CG-PRP、CG-HS、CG-LS 和L-BFGS 五种算法在三维悬沙输运伴随同化模型中优化反演初始场的表现。设计孪生实验TE11—E15， 以图2a 所示的初始场作为给定初始场， 初始场的初始猜测值取为给定初始场的水平空间平均值的一半， 即常数0.03325kg/m3， 其余模型参数均为默认值。从反演结果和误差分布(图4)来看， 五种算法的反演结果与给定分布较为吻合， 利用GD、CG-PRP 和CG-LS算法反演所得初始场的误差主要集中分布在杭州湾东侧湾口的海洋开边界处， 而在其余的大部分区域误差较小， 均取得了良好的反演效果； 然而利用CG-HS 和L-BFGS 算法反演所得初始场的误差在湾顶、湾口和南部的大部分区域都很大， 反演效果较差， 且CG-HS 算法要比L-BFGS 算法更差。从统计结果(图5、表2)来看， 五组实验中约化代价函数、SSC 模拟结果与“观测数据”之间的平均绝对误差(mean absolute error， MAE)均有不同程度的下降， 可以定量说明数据同化的效果(Jinet al，2015； Liuet al， 2017； Zonget al， 2018)。以约化代价函数为例， 实验TE15的约化代价函数下降速度最快， 但只下降了2 个数量级； 实验TE11、TE12和TE14的约化代价函数下降速度相当， 均下降了3 个数量级； 实验TE13的约化代价函数下降速度最慢， 在第3 步已达到收敛准则， 迭代停止。此外， 代价函数关于初始场的梯度的L1范数、反演结果与给定分布之间的MRE 两个统计量可以直观地说明初始场的反演效果。实验TE11、TE12和TE14中两个统计量分别经过171、181、181 步均大幅度下降， 表明初始场的反演效果较好； 而实验TE13和TE15中反演结果的误差较大， 反演效果较差。综合以上分析， GD、CG-PRP 和CG-LS 三种算法反演初始场的效果相当， 但GD 算法具有易于实现、迭代平稳易于控制、良好的可操作性等优点， 在本模型中具有很大的优势； L-BFGS 算法的反演速率最快， 但最终的反演结果较差； CG-HS 算法无论是反演速率还是反演结果都是最差的。以上数值实验的结论与Chen 等(2013)和Zhang 等(2015)的研究结论相似。因此， 本节之后的孪生实验和实际实验选用GD 算法用于模型初始场的优化反演。

3.3 初始猜测值对初始场优化反演的影响

本节设计孪生实验TE21—TE23研究初始场的初始猜测值对反演结果的影响。初始场的初始猜测值分别取为给定初始场的水平空间平均值的1/2、给定初始场的水平空间平均值和给定初始场的水平空间平均值的3/2， 即常数0.03325、0.0665 和0.09975kg/m3。从反演结果和误差分布(图6)来看， 三组实验的反演结果相差不大， 反演结果与给定分布吻合良好， 相对误差在大部分区域均小于20%。从统计结果(图7、表3)来看， 三组实验中相关统计量均大幅度下降。其中， 约化代价函数均下降了至少3 个数量级， 代价函数关于初始场的梯度的1L范数均至少下降了2 个数量级， SSC 模拟结果与观测之间的MAE 均至少下降了1 个数量级， 初始场的反演结果与给定分布之间的MRE 均小于16%， 表明初始场的反演结果充分收敛到了“真实值”附近。本节数值实验结果表明， 初始场的反演效果受 初始猜测值的影响较小， 初始场的“真实值”均能得到精确反演。因此， 本节之后的敏感性理想实验中初始猜测值均取为0.03325kg/m3。

图6 TE21—TE23 中反演后的表层初始场(a—c)、反演后的表层初始场与给定初始场的误差(d—f)、反演后的表层初始场与给定初始场的相对误差(g—i)的空间分布 Fig.6 The spatial distributions of the surface initial field after inversion in TE21—TE23 (a—c)； the spatial distributions of errors between the surface initial field after inversion and the given initial field in TE21—TE23 (d—f)； and the spatial distributions of relative errors between the surface initial field after inversion and the given initial field in TE21—TE23 (g—i)

图7 实验TE21—TE23 中log10(J/J0) (a)、Gradient of C0 (b)、MAEs at AOs(c)、MAEs at COs (d)、MREs of C0(e)随迭代步的变化 Fig.7 The values of log(J/J0) versus the iteration steps (a)； the variations of Gradient of C0 (b)； the variations of MAEs at AOs (c)； the variations of MAEs at COs (d)； and the variations of MREs of C0 (e) in TE21—TE23

表3 实验TE21—TE23 中， 数据同化前后的log10(J/J0)、Gradient of C0、MAEs at AOs、MAEs at COs 和MREs of C0 Tab.3 log10(J/J0)， gradient of C0， MAEs at AOs， MAEs at Cos， and MREs of C0 before and after data assimilation in TE21—TE23

3.4 卫星遥感数据空间分布数量对初始场优化反演的影响

在实际的卫星遥感观测中， 由于云覆盖或恶劣天气条件等因素的影响， 卫星遥感数据往往存在一定程度的缺失， 因此有必要讨论卫星遥感数据数量对初始场优化反演的影响。本节设计孪生实验TE31—TE34， 四组实验中“被同化观测”的数量依次随机减少1/10、1/4、1/2、3/4。从反演结果和误差分布(图8)可以看出， 随着卫星数据数量的减少， 反演精度相差不大， 反演结果在大部分区域与给定分布吻合良好， 误差较大的区域主要集中在湾口的海洋开边界处。从统计结果(图9、表4)可以看出， 反演结果与给定分布之间的MRE 分别为10.94%、10.82%、12.71%和 13.62%， 依然取得了较高的反演精度， 实验TE34相比于实验TE31初始场的反演结果与给定分布之间的MRE 只增大了2.68%， 相差不大。据此我们认为， 在目前的实验设置下， 卫星遥感数据数量对初始场的优化反演影响较小， 即便是分布稀疏或云覆盖现象较为严重， 依然能够取得良好的反演效果， 表明伴随同化方法是还原SSC 空间分布的有效手段。

图8 TE31—TE34 中反演后的表层初始场(a—d)、反演后的表层初始场与给定初始场的误差(e—h)、反演后的表层初始场与给定初始场的相对误差(i—l)的空间分布 Fig.8 The spatial distributions of the surface initial field after inversion in TE31—TE34 (a—d)； the spatial distributions of errors between the surface initial field after inversion and the given initial field in TE31—TE34 (e—h)； and the spatial distributions of relative errors between the surface initial field after inversion and the given initial field in TE31—TE34 (i—l)

图9 实验TE31—TE34 中log10(J/J0) (a)、Gradient of C0 (b)、MAEs at AOs (c)、MAEs at COs (d)、MREs of C0(e)随迭代步的变化 Fig.9 The values of log(J/J0) versus the iteration steps (a)； the variations of Gradient of C0 (b)； the variations of MAEs at AOs (c)； the variations of MAEs at COs (d)； and the variations of MREs of C0 (e) in TE31—TE34

表4 实验TE31—TE34 中， 数据同化前后的log10(J/J0)、Gradient of C0、MAEs at AOs、MAEs at COs 和MREs of C0 Tab.4 log10(J/J0)， gradient of C0， MAEs at AOs， MAEs at COs， and MREs of C0 before and after data assimilation in TE31—TE34

3.5 卫星遥感数据同化时间窗口的宽度对初始场优化反演的影响

本节设计孪生实验TE41—TE44研究卫星遥感数据同化时间窗口的宽度对初始场优化反演的影响。SSC 数值模拟的起始时刻是GOCI 卫星遥感观测数据同化时间窗口的起始时刻， 同化窗口依次增大， 分别为1、2、3、4h， 即分别覆盖2 次、3 次、4 次、5 次的GOCI 卫星遥感观测。从反演结果和误差分布(图10)可以看出， 随着同化窗口的宽度的增大， 初始场的反演效果越来越好， 当同化窗口的宽度增加到3h， 在杭州湾中部和东南部两个误差相对较大的区域得到明显改善。从统计结果(图11、表5)可以看出， 四组实验的约化代价函数、代价函数关于初始场的梯度的1L范数以及模拟结果和观测之间的MAE 分别经过200、200、194、191 步均显著下降， 四组实验初始场的反演结果与给定分布之间的MRE 分别为13.53%、10.03%、8.15%和8.73%， 即增大同化窗口的宽度可以有效改善反演结果。由此可见， 初始场的优化反演效果对同化窗口的宽度较为敏感， 当SSC 数值模拟的起始时刻是GOCI 卫星遥感观测数据同化窗口的起始时刻， 只需要至少3h 的同化窗口即可得到改善效果明显的反演结果。因此， 对于SSC 的数值模拟而言， 合理设置卫星遥感数据同化时间窗口的宽度可以有效提高初始场的优化反演效果。

3.6 背景水动力场的误差对初始场优化反演的的影响

图10 TE41—TE44 中反演后的表层初始场(a—d)、反演后的表层初始场与给定初始场的误差(e—h)、反演后的表层初始场与给定初始场的相对误差(i—l)的空间分布 Fig.10 The spatial distributions of the surface initial field after inversion in TE41—TE44 (a—d)； the spatial distributions of errors between the surface initial field after inversion and the given initial field in TE41—TE44 (e—h)； and the spatial distributions of relative errors between the surface initial field after inversion and the given initial field in TE41—TE44 (i—l)

图11 实验TE41—TE44 中log10(J/J0) (a)、Gradient of C0 (b)、MAEs at AOs (c)、MAEs at COs (d)、MREs of C0(e)随迭代步的变化 Fig.11 The values of log(J/J0) versus the iteration steps (a)； the variations of Gradient of C0 (b)； the variations of MAEs at AOs (c)； the variations of MAEs at COs (d)； and the variations of MREs of C0 (e) in TE41—TE44

表5 实验TE41—TE44 中， 数据同化前后的log10(J/J0)、Gradient of C0、MAEs at AOs、MAEs at COs 和MREs of C0 Tab.5 log10(J/J0)， gradient of C0， MAEs at AOs， MAEs at COs， and MREs of C0 before and after data assimilation in TE41—TE44

本文中背景水动力场是由FVCOM 模拟生成， 用于驱动三维悬沙输运模型。由于背景流场的模拟结果与实际流场不可避免地存在误差， 因而三维悬沙输运模型的模拟精度会很大程度上依赖于背景水动力场的模拟精度(Mehtaet al， 1989， Xieet al， 2009， Wanget al，2018b)， 因此有必要讨论背景水动力场的误差对初始场优化反演的影响。李丹(2016)利用FVCOM 对长 江口杭州湾海域的水动力场进行模拟， 将模拟结果与实测资料进行对比， 潮流流速误差的最大值在15%以内。Wang 等(2018b)利用FVCOM 对杭州湾海域水动力场的模拟结果显示， 模拟的潮流流速的相对误差大多不超过20%。因此， 本节设计孪生实验TE51—TE53， 三组实验中分别对背景水动力场中流速的三个分量u、v、w均添加最大不超过10%、20%和30%的随机误差， 误差添加方式与3.1 节中相同。从反演结果和误差分布(图12)来看， 三组实验中初始场的反演结果相差不大， 反演结果均与给定分布吻合良好。从统计结果(图13、表6)来看， 三组实验相关统计量分别经过179、179、180 步都下降到了相同的数量级， 三组实验中反演结果与给定分布之间的MRE 均小于12%。因此，背景水动力场的误差对初始场优化反演的影响较小， 本文所采用的FVCOM 模拟结果足以满足研究要求。

图12 TE51—TE53 中反演后的表层初始场(a—c)、反演后的表层初始场与给定初始场的误差(d—f)、反演后的表层初始场与给定初始场的相对误差(g—i)的空间分布 Fig.12 The spatial distributions of the surface initial field after inversion in TE51—TE53 (a—c)； the spatial distributions of errors between the surface initial field after inversion and the given initial field in TE51—TE53 (d—f)； and the spatial distributions of relative errors between the surface initial field after inversion and the given initial field in TE51—TE53 (g—i)

图13 实验TE51—TE53 中log10(J/J0) (a)、Gradient of C0(b)、MAEs at AOs (c)、MAEs at COs (d)、MREs of C0(e)随迭代步的变化 Fig.13 The values of log(J/J0) versus the iteration steps (a)； the variations of Gradient of C0 (b)； the variations of MAEs at AOs (c)； the variations of MAEs at COs (d)； and the variations of MREs of C0 (e) in TE51—TE53

表6 实验TE51—TE53 中， 数据同化前后的log10(J/J0)、Gradient of C0、MAEs of AOs、MAEs of COs 和MREs of C0 Tab.6 log10(J/J0)， gradient of C0， MAEs at AOs， MAEs at COs， and MREs of C0 before and after data assimilation in TE51—TE53

3.7 伴随同化方法和插值方法重构初始场的能力比较

为了比较伴随同化方法和插值方法重构模型初始场的能力， 我们构建一个新的初始场(图14a)， 该初始场由图 2a 所示初始场运行正向模型所得与实际GOCI 卫星遥感观测时刻对应的“观测数据”也即3.1节中的“观测数据”经过时间平均后获得。利用新的初始场运行正向模型， 获得与实际观测数据一一对应的受云覆盖影响的“被同化观测”和“检验观测”用于本节的数据同化和初始场重构。在伴随同化方法中， 初始场的初始猜测值取为新的初始场的水平空间平均值的一半， 即常数0.0535kg/m3， 卫星遥感数据同化时间窗口的宽度设置为3h， 其余模型设置不变。在插值方法中， 分别利用双三次插值法、自然邻近插值法和最近邻点插值法插值数值模拟初始时刻的“观测数据”， 获得重构的初始场。从初始场的重构结果(图14b—e)和误差分布(图14f—m)可以看出， 伴随同化方法重构的初始场与给定分布最接近， 在绝大部分区域二者吻合良好， 很好地重现了SSC 变化剧烈的区域， 且保持了很好的平滑性， 重构结果与给定分布之间的MRE 只有11.53%； 双三次插值法和自然邻近插值法重构初始场的效果较差， 虽然总体保持了与给定分布一致的空间变化规律， 但在杭州湾中西部、南部陆地边界处和中东部存在三个相对误差局部高值区， 插值结果对三处SSC 的剧烈变化刻画不够细致； 最近邻点插值法重构的初始场与给定分布偏差最大， 在上述三个局部区域SSC 均远高于给定初始场， 重构结果与给定分布之间的MRE 高达46.95%， 且插值结果的平滑性较差， 与实际偏差较大。因此， 与传统的插值方法相比， 伴随同化方法是重构模型初始场更有效的手段。

4 实际实验及结果分析

图14 给定的(a)和重构的(b—e)表层初始场的空间分布； 重构的表层初始场与给定初始场的误差(f—i)以及相对误差(j—m)的空间分布 Fig.14 The spatial distribution of the given surface initial field (a)； the spatial distributions of initial field of the surface layer reconstructed by adjoint assimilation method， bicubic interpolation method， natural neighbor interpolation method， and nearest neighbor interpolation method (b—e)； the spatial distributions of errors between the surface initial field of the reconstructed surface initial field and the given initial field (f—i)； and the spatial distributions of relative errors of the reconstructed surface initial field and the given initial field (j—m)

实际实验中， 分别同化典型的小潮时期和大潮时期GOCI 卫星遥感资料所得杭州湾海域表层SSC 数据， 利用伴随同化方法对模型初始场进行优化反演。根据孪生实验的结论， 优化算法选用GD 算法， 初始场的初始猜测值取为初始时刻实测的GOCI 卫星遥感观测数据经过空间插值获得的SSC(图2)， 卫星遥感数据同化时间窗口的宽度设置为3h， 其余模型设置不变。从图15 可以看出， 实际初始场缺测严重， 在缺测区域优化反演后的初始场(图16)与反演前的插值初始场(图2)有很大不同。在小潮时期， 反演前的初始场的SSC 高值区出现在杭州湾的西部、中部和东南部， 而反演后的初始场的SSC 高值区出现在杭州湾的东南部及南部沿岸的大部分地区。在大潮时期， 反演前的初始场的SSC 高值区出现在杭州湾的中南部、东部和东南部的大部分区域， 而反演后的初始场的SSC高值区出现在杭州湾的东南部及南部沿岸的大部分地区。无论是小潮时期还是大潮时期， 反演后的初始场的SSC 均呈现出由南到北逐渐递减的趋势， 这与基于一个潮汐周期得到的杭州湾GOCI 卫星遥感观测中SSC“北低南高”的分布特征相似(Heet al， 2013； 江彬彬等， 2015； Liuet al， 2018)， 显然优化反演后的初始场更符合实际， 证明了伴随同化方法优化模型初始场的有效性。从图17 可以看出， 优化模型初始场可以很大程度上提高SSC 的模拟精度。小潮(大潮)时期， 经过260(102)步迭代， 满足收敛标准， 约化代价函数、代价函数关于初始场的梯度的L1范数均大幅下降并趋于平稳， “被同化观测”和模拟结果之间的MAE 从8.72×10-2kg/m3(3.66kg/m3)下降到了2.72×10-2kg/m3(7.25×10-1kg/m3)， 下降了68.81%(80.19%)； “检验观测”和模拟结果之间的 MAE 从 9.02×10-2kg/m3(3.63kg/m3)下降到了2.88×10-2kg/m3(7.39×10-1kg/m3)， 下降了68.07%(79.64%)。可见， 反演后SSC的模拟结果和观测之间的误差大幅度降低， 模拟结果与观测更接近。以上结果表明， 初始场对SSC 的模拟结果有重要的影响， 利用伴随同化方法可以获得模型的最优初始场， 从而很大程度上改善SSC 的模拟结果。

图15 小潮(a)和大潮(b)时期实际初始场的空间分布 Fig.15 The spatial distribution of actual initial field during neap tide (a) and spring tide(b)

图16 实际实验中， 小潮(a)和大潮(b)时期初始场的反演结果 Fig.16 The inversion result of initial field during neap tide (a) and spring tide (b) in practical experiments

图17 实际实验小潮时期的log10(J/J0)(a)、Gradient of C0(b)、MAEs at AOs(c)、MAEs at COs (d) 随迭代步的变化； 实际实验大潮时期的log10(J/J0) (e)、Gradient of C0 (f)、MAEs at AOs (g)、MAEs at COs (h) 随迭代步的变化 Fig.17 The values of log(J/J0) (a， e)， the variations of Gradient of C0 (b， f)， the variations of MAEs at AOs (c， g)， and the variations of MAEs at COs (d， h) during neap tide and spring tide in practical experiments

根据以上分析， 实际初始场缺测严重时， 利用插值方法获得的初始场(图2)不能体现杭州湾南部沿岸高SSC 的特征， 且插值结果的平滑性较差。而利用伴随同化方法反演所得初始场在有观测区域与实际初始场较为一致， 且在初始场的观测数据缺测区域， 伴随同化方法的反演效果要更好， 很好地反演出了SSC“北低南高”的分布特征， 且总体上保持了很好的 平滑性， 更符合实际。这是因为插值方法仅考虑了SSC 的空间分布， 未考虑时间变化， 不能充分地利用观测资料， 因而不能较好地重现SSC 的时空变化特征， 特别是在云覆盖较为严重的区域， 传统的插值方法往往是失效的， 模拟所得SSC 与实际偏差较大； 而伴随同化方法既考虑SSC 的空间分布， 又考虑其时间变化， 可以获得与实际观测更加接近的初始场， 因而模拟结果更接近观测值， 且更符合实际。

5 结论

本文基于一个三维悬沙输运伴随同化模型， 通过一系列孪生实验和实际实验， 对模型初始场进行优化反演研究。

孪生实验中， 首先验证了三维悬沙输运模型对于初始场的相对敏感性， 引出了利用伴随同化方法优化初始场的必要性。其次， 探讨了初始场对于不同影响因素的相对敏感性， 包括优化算法、初始猜测值、卫星遥感数据数量、卫星遥感观测数据同化时间窗口的宽度和背景流场误差。结果表明， 在五种优化算法中， GD 算法是最优的； 对本模型设置而言， 初始猜测值对反演效果影响较小， 初始场的“真实值”均能得到精确反演； 初始场对卫星遥感数据数量较不敏感， 即便是分布稀疏或云覆盖现象较为严重， 依然能够取得良好的反演效果； 初始场对同化窗口的宽度较为敏感， 合理设置同化窗口的宽度可以有效改善初始场的优化反演效果； 初始场对背景流场的误差较不敏感， 即便是30%的随机误差， 反演结果与给定分布之间的MRE 依然小于12%。最后， 比较了伴随同化方法和插值方法重构模型初始场的能力。结果表明， 与传统的插值方法相比， 伴随同化方法是重构模型初始场更有效的手段。

实际实验中， 利用孪生实验的相关结论， 在杭州湾海域同化典型的小潮时期和大潮时期的GOCI 卫星遥感资料所得表层SSC 数据， 优化反演了模型初始场。结果表明， 无论是小潮时期还是大潮时期， 利用伴随同化方法均得到了更加符合实际的最优初始场， 使得SSC 的数值模拟结果和实测数据更加接近。

该研究进一步证实了伴随同化方法优化模型初始场的有效性， 并为改进悬沙输运模型及其他海洋数值模型的初始化方案提供了一种有效的方法。未来我们将尝试同时同化更多、更长时间尺度的观测资料， 为悬沙输运的数值模拟提供更为精确的初始场。