A modified cressman interpolation method applied to process total nitrogen data in Bohai surface waters
-
摘要:
根据渤海常规监测站的站位分布特征，提出了一种用来处理渤海表层总氮数据的改进的Cressman插值（MCI）方法。由渤海常规监测总氮数据反应的污染状况设定污染物浓度的两种理想分布，通过理想实验发现利用MCI方法插值得到的误差比传统的Cressman方法、Kriging方法和Local Polynomial（LP）方法更小，验证了MCI方法的有效性。实际实验利用交叉验证方法来计算插值的误差，并用MCI方法插值出2009年和2010年共4个月份的表层总氮空间分布，为渤海污染状况分析提供较为准确的数据支持。
Abstract:A modified Cressman interpolation (MCI) method was proposed to study the total nitrogen (TN) data in Bohai surface waters, based on the distribution of Bohai routine monitoring stations.The pollution concentration station data features two idealized distributions.The twin experiments suggest reduced interpolation biases using MCI method, compared to the traditional Cressman, Kriging and Local Polynomial (LP) method.In practical experiments, cross validation scheme was applied to evaluate the interpolation biases.Bohai Surface TN concentration was recalculated following MCI method for four quarters in 2009 and 2010, providing an improved dataset for Bohai pollution investigation.
-
Keywords:
- interpolation radius /
- cross validation /
- Bohai Sea /
- total nitrogen
-
对于环境污染状况评估而言，精确充足的数据是必不可少的。尽管可以通过一些监测手段获取到重要的污染物数据，但是此类数据往往是离散的。这些离散数据必须要进行整合(比如数据拟合、插值等)得到连续的再分析数据才能为我们进一步的研究提供依据。
目前已有许多插值方法，例如反距离插值法，Kriging方法和LP方法等应用于诸多领域[1-5]。Kriging方法通过对测量数据做统计分析来描述空间结构特征，但是这种方法较为复杂，并且计算量较大[6]；作为反距离插值方法的一类，Cressman[7]插值方法则比较灵活，且易于实现。LP方法是一种局部加权最小二乘方法，既考虑全部数据点的趋势性变化，又能较好地反映局部特征[8]，但在数据稀少时易造成矩阵奇异性。在样本点距离太远或者严重不足的情形下需要对插值方法进行改进才能得到更好的插值效果[9-15]。本文针对渤海常规污染物监测站位的布设特点和表层总氮浓度的分布特征，提出了一种改进的Cressman方法来处理渤海表层总氮数据，为渤海污染状况分析提供较为准确的数据支持。
1 材料与方法
本文收集到的数据是由北海环境监测中心提供的渤海常规监测总氮数据。2009年和2010年4个月份(2009年5月、8月、10月和2010年5月)的渤海常规监测站点的表层总氮浓度分布如图 1所以，其中圆点的大小跟总氮浓度成正比。从图上我们可以看出，近岸处监测站位分布比较密集，且浓度较高；而在渤海中部区域，站点分布较为稀疏，浓度较小。这种特殊的监测站位分布有其合理性，能够利用较少的观测站位信息来初步获取渤海海域的污染物空间分布特征。要得到其它非监测站位的总氮浓度信息，需要通过其他有效方法来实现，比如插值方法。本文利用Cressman方法、Kriging方法和LP方法这3种插值方法来计算渤海表层总氮的空间分布。
![]()
图 1 渤海常规监测站点的表层总氮浓度分布Fig. 1 Distribution of surface TN concentration of Bohai routine monitoring stationsCressman方法和Kriging方法可以统一用如下公式表示：
(1) 其中：z(xi)代表目标网格点x0的数据；z(xi)代表采样网格点xi的数据；λi代表权重系数；n表示采样点的个数[16]。Cressman插值方法的权重系数λi的计算公式为
(2) 其中：
(3) 其中：R表示插值半径；ri表示是采样点和目标点之间的距离。如果用传统的Cressman插值方法对图 1中的数据进行插值，插值半径必须足够大才能保证所有区域都能被插值到，这就会导致插值结果出现较大的误差。同时注意到，渤海中部区域站点污染物浓度很小并且差距不大，为此对传统的Cressman插值进行如下改进：
1) 将渤海中部站点的表层总氮浓度较小(＜0.01 mg/L)的数据做平均，作为整个浓度场的背景值；
2) 将所有的观测数据减去背景值，得到一组新的数据；
3) 对处理后的上述数据做传统的Cressman插值，得到一组插值后的表层总氮浓度场；
4) 将第一步得到的背景值加到上述方法得到的浓度场中，即可插值出整个渤海表层总氮浓度分布。
利用等值线绘图软件Surfer对数据进行Kriging插值和LP插值，该软件具有的强大插值功能和绘制图件能力，能迅速地将离散的数据通过插值转换为连续的三维曲面。Kriging方法通过建立变异函数理论模型，通过采样数据值的初始设置和变异函数的空间特性来对每个插值点做最优无偏估计[17]。Kriging的核心是建立变异函数，其定义为：
(4) 其中：N(h)是距离为h的点对数；Z为特征值；r(h)为变异函数值。其中，公式(1) 中的λi满足如下条件：
(5) 其中：μ表示Lagrange系数；γ表示观测点之间的距离。
LP方法是一种局部加权最小二乘拟合方法，在搜索范围内能对所有数据点选择适当特定阶数的多项式[18]。设某一阶的多项式函数
(6) 根据数据点的分布确定搜索范围，利用最小二乘原理:
(7) 其中：Wi为权值；Z(Xi，Yi)是数据点(Xi，Yi)的值。若多项式函数确定，则目标网格点的值也确定。
2 结果与讨论
2.1 理想实验
为验证MCI方法在渤海常规监测站位的分布特征下插值的有效性，根据渤海表层总氮浓度初步反应的分布特征设置污染物浓度的两种理想分布，如图 2所示，其中子图(a)和(b)描述的分布分别记为分布1和分布2。从理想分布中挑选出跟监测站点位置相同的污染物浓度值进行插值，来比较MCI方法和其他传统插值方法的优劣。其中，网格分辨率设置为1′×1′。
![]()
图 2 设定渤海表层污染物浓度的两种理想空间分布Fig. 2 Two prescribed distributions of Bohai Surface pollution concentration因部分监测站位在渤海中部区域的站位相距较远，用传统的Cressman插值时，为了保证至少有两个站位包含在插值半径之内，插值半径必须足够大。插值半径的求法如下：首先将插值半径设置为1.5°，如果在此插值半径内包含两个或者两个以上的站位，那么不需要再对插值半径进行调整；否则，将插值半径在原来的基础上增加0.1°，直至找到符合要求的半径为止。
MCI方法与传统的Cressman插值相比，插值半径的选择较为灵活，可以避免插值半径过大造成误差较大，本文将插值半径设置为0.1°到0.9°。利用MCI得到的插值结果和给定分布的误差如图 3所示。从图 3可以看出，插值半径对误差的影响明显，因此需要对插值半径进行多次尝试以得到最优插值结果，最优插值结果对应下的插值半径定义为最优插值半径。不用分布对应的最优插值半径也不同，分布1的最优插值半径为0.4°，而分布2则为0.3°；同一插值半径下分布1的误差要大于分布2，这可能是因为分布1的空间结构变化更为剧烈。
![]()
图 3 不同插值半径下利用MCI方法的插值结果跟给定分布的误差Fig. 3 Errors between interpolation results with MCIM and the prescribed distributions among different interpolation radii表 1表示最优插值半径下的给定分布和插值结果的误差统计结果。从表 1中能够清晰的看出，利用MCI方法插值得到的误差要小于传统的Cressman，Kriging方法和LP方法，说明了MCI方法的有效性。而从插值之后的污染物分布和给定分布对比来看(图 4~7)，用MCI方法插值得到的分布也非常接近给定分布，只在局部出现了微小差异，局部误差最大值不超过0.25 mg/L。用Cressman方法插值得到的分布总体来看比较光滑，但是与给定分布相比差距很大，尤其是在渤海中部区域，浓度值大于0.3 mg/L，而给定分布在渤海中部浓度值小于0.1 mg/L。利用Kriging方法得到的分布尽管看上去与给定分布的变化趋势较为一致，但是数值偏小且误差较大，局部地区的插值误差甚至达到1.0 mg/L。用LP方法插值得到的分布与给定分布相比差距极为明显，误差最大。
表 1 插值之后的结果跟给定分布的误差(单位：mg·L-1)Tab. 1 Errors of interpolation results and prescribed distributions (unit:mg·L-1)
2.2 实际实验
从2.1部分的理想试验可以看出，MCI方法能够利用渤海常规监测站点的位置信息成功的插值出给定的污染物分布，并且误差要小于传统的Cressman方法，Kriging方法和LP方法，从而MCI方法可应用在实际情形中。
实际实验中利用交叉验证方法来统计插值误差。首先，从渤海常规监测总氮数据中选取约10%(共选取10次)的数据作为验证数据，剩余数据用来插值。在每个插值半径下，验证数据和插值结果的误差都要统计10次，将这10次的平均误差作为插值的误差。不同插值半径下利用MCI方法得到的插值结果跟验证数据的误差如表 2所示。图 8表示最优插值半径下利用MCI方法得到的插值结果。
表 2 利用MCI方法得到的插值结果跟验证数据的误差(单位：mg·L-1)Tab. 2 Errors between interpolation results and data with MCI method (unit: mg·L-1)
(1) 不同月份下，MCI方法在插值时的最优插值半径及误差也会不同。例如2009年8月份的最优插值半径是0.2°，误差为0.115 mg/L, 而在2009年10月份则为0.3°，误差为0.082 mg/L。
(2) 从整体上看，渤海表层总氮浓度分布呈现三湾较高而渤海中部及渤海海峡附近较低的分布趋势。三湾中莱州湾的浓度最高(最高浓度值>4.5 mg/L)，这一原因可能是莱州湾海流较弱致使海水交换能力不强，加之黄河等河流的输入，导致浓度较高；渤海湾次之，该海域海流虽较强，但有海河等大型河流的输入[19]；辽东湾的表层总氮浓度最低。渤海中部因受陆源影响小导致浓度较低，渤海海峡附近与外界海水交换能力较强导致浓度较低。
(3) 渤海表层总氮浓度季节性变化明显。夏秋季节渤海表层总氮的浓度较高，可达4.5 mg/L以上；而在春季浓度较低，浓度最高不到1.8 mg/L。夏秋季节降水较多，通过河流注入等方式进入到渤海的污水较多；加之农作物施肥、喷洒农药等原因，更加剧了渤海的污染程度。
3 结论
(1) 本文收集到2009年5月、8月、10月和2010年5月共4个月份的渤海常规监测站的表层总氮数据，根据站位分布特征提出了一种处理渤海表层总氮数据的MCI方法。MCI方法通过引入“背景值”避免插值半径选择多大，从而减小了插值误差。为验证MCI方法的有效性，首先根据渤海表层总氮监测数据初步反应的污染情形设置污染物浓度的两种理想分布，实验发现MCI方法插值得到的误差要小于传统的Cressman方法、Kriging方法和LP方法，从而对实际实验做出一定参考。
(2) 实际实验利用交叉验证方法来计算MCI方法插值时的误差，并分别插值出2009年和2010年共4个月份的表层总氮空间分布，对渤海总氮浓度的时空分布规律做出分析，从而为渤海污染状况分析提供较为准确的数据支持。
-
![]()
图 1 渤海常规监测站点的表层总氮浓度分布
Fig. 1. Distribution of surface TN concentration of Bohai routine monitoring stations
![]()
图 2 设定渤海表层污染物浓度的两种理想空间分布
Fig. 2. Two prescribed distributions of Bohai Surface pollution concentration
![]()
图 3 不同插值半径下利用MCI方法的插值结果跟给定分布的误差
Fig. 3. Errors between interpolation results with MCIM and the prescribed distributions among different interpolation radii
表 1 插值之后的结果跟给定分布的误差(单位：mg·L-1)
Tab. 1 Errors of interpolation results and prescribed distributions (unit:mg·L-1)
下载: 导出CSV
表 2 利用MCI方法得到的插值结果跟验证数据的误差(单位：mg·L-1)
Tab. 2 Errors between interpolation results and data with MCI method (unit: mg·L-1)
下载: 导出CSV
-
[1] LARGUECHE F Z B.Estimating soil contamination with kriging interpolation method[J].American Journal of Applied Sciences, 2006, 3(6):1894-1898. doi: 10.3844/ajassp.2006.1894.1898
[2] 阳文锐, 王如松, 黄锦楼, 等.反距离加权插值法在污染场地评价中的应用[J].应用生态学报, 2007, 18(9):2013-2018. http://www.cnki.com.cn/Article/CJFDTOTAL-YYSB200709016.htm [3] BARGAOUI Z K, CHEBBI A.Comparison of two kriging interpolation methods applied to spatiotemporal rainfall[J].Journal of Hydrology, 2009, 365(1/2):56-73.
[4] CHIAO L Y, CHEN Y N, GUNG Y.Constructing empirical resolution diagnostics forkriging and minimum curvature gridding[J].Journal of Geophysical Research, 2014, 119(5):3939-3954.
[5] VERDIN A, FUNK C, RAJAGOPALAN B, et al.Kriging and local polynomial methods for blending satellite-derived and gauge precipitation estimates to support hydrologic early warning systems[J].IEEE Transactions on Geoscience and Remote Sensing, 2016, 54(5):2552-2562. doi: 10.1109/TGRS.2015.2502956
[6] KRAVCHENKO A N.Influence of spatial structure on accuracy of interpolation methods[J].Soil Science Society of America Journal, 2003, 67(5):1564-1571. doi: 10.2136/sssaj2003.1564
[7] CRESSMAN G P.An operational objective analysis system[J].Monthly Weather Review, 1959, 87(10):367-374. doi: 10.1175/1520-0493(1959)087<0367:AOOAS>2.0.CO;2
[8] 汪俊, 高金耀, 吴招才, 等.局部多项式插值方法在多源海底沉积厚度数据融合中的应用[J].海洋科学, 2009, 33(4):25-28. http://www.cnki.com.cn/Article/CJFDTOTAL-HYKX200904008.htm [9] GU L.Moving kriging interpolation and element-free Galerkin method[J].International Journal for Numerical Methods in Engineering, 2003, 56(1):1-11. doi: 10.1002/(ISSN)1097-0207
[10] TONGSUK P, KANOK-NUKULCHAI W.Further investigation of element-free Galerkin method using moving kriging interpolation[J].International Journal of Computational Methods, 2004, 1(2):345-365. doi: 10.1142/S0219876204000162
[11] 张静, 李旭祥, 蔡启闽, 等.非参数局部多项式法在大气环境数据分析中的应用[J].环境工程, 2010, 28(S):343-346. http://www.cnki.com.cn/Article/CJFDTOTAL-HJGC2010S1091.htm [12] 周永道, 马洪, 吕王勇, 等.基于多元局部多项式方法的混沌时间序列预测[J].物理学报, 2007, 56(12):6809-6814. doi: 10.3321/j.issn:1000-3290.2007.12.006 [13] 严华雯, 吴健平.加权最小二乘法改进遗传克里金插值方法研究[J].计算机技术与发展, 2012, 22(3):92-95. http://www.cnki.com.cn/Article/CJFDTOTAL-WJFZ201203025.htm [14] MONTEYNE G, UGRYUMOVA D, VANDERSTEEN G.Extension of local Polynomial method for periodic excitations[J].IFAC Proceedings Volumes, 2012, 45(16):61-65. doi: 10.3182/20120711-3-BE-2027.00015
[15] HUANG C C, ZHENG X G, TAIT A, et al.On using smoothing spline and residual correction to fuse rain gauge observations and remote sensing data[J].Journal of Hydrology, 2014, 508:410-417. doi: 10.1016/j.jhydrol.2013.11.022
[16] WEBSTER R, OLIVER M A.Geostatistics for environmental scientists[M].2nd ed.Chichester, England:John Wiley & Sons Ltd, 2007.
[17] 彭兆璇, 袁峰, 周涛发, 等.土壤中元素空间分布的体视化方法研究[J].计算机技术与发展, 2009, 19(5):195-197. http://www.cnki.com.cn/Article/CJFDTOTAL-WJFZ200905054.htm [18] 陈欢欢, 李星, 丁文秀.Surfer 8.0等值线绘制中的十二种插值方法[J].工程地球物理学报, 2007, 4(1):52-57. http://www.cnki.com.cn/Article/CJFDTOTAL-GCDQ200701010.htm [19] SHEN Y L, WANG C H, WANG Y G, et al.Inversion study on pollutant discharges in the Bohai Sea with the adjoint method[J].Journal of Ocean University of China, 2015, 14(6):941-950. doi: 10.1007/s11802-015-2501-8
-
期刊类型引用(0)
其他类型引用(1)
计量
- 文章访问数: 2766
- HTML全文浏览量: 1885
- PDF下载量: 59
- 被引次数: 1
