全國近岸海域水質(zhì)空間插值算法精度分析

李俊龍，丁 頁，劉喜惠，劉 方，李 曌，陳 平，劉海江，李銘煊，高鋒亮

中國環(huán)境監(jiān)測總站，國家環(huán)境保護環(huán)境監(jiān)測質(zhì)量控制重點實驗室，北京 100012

近岸海域環(huán)境要素的空間分布和表征技術(shù)研究是開展環(huán)境監(jiān)測和評價，構(gòu)建中國近岸海域環(huán)境監(jiān)測技術(shù)體系的重要組成部分［1］?？臻g插值技術(shù)形成于20世紀(jì)50年代，理論體系不斷完善［2-4］，尤其在大尺度環(huán)境要素的分析和表征方面，空間插值技術(shù)具有明顯的優(yōu)勢和價值［5-7］。近岸海域環(huán)境要素具有空間相關(guān)性、依賴性和異質(zhì)性，因此，在現(xiàn)階段人力、財力和技術(shù)裝備能力有限的情況下，空間插值技術(shù)對于評估采樣點周圍空白區(qū)域的環(huán)境屬性信息十分重要［8-10］。隨著地理統(tǒng)計科學(xué)的迅速發(fā)展，空間插值算法已有很多，但各種算法與近岸海域環(huán)境監(jiān)測站位的匹配性、準(zhǔn)確性和適用性分析成為難題，相關(guān)研究也較少見，制約了空間插值技術(shù)在環(huán)境監(jiān)測領(lǐng)域的應(yīng)用。為提高空間插值算法應(yīng)用的科學(xué)性，更準(zhǔn)確地評價和表征近岸海域環(huán)境質(zhì)量，基于實測歷史數(shù)據(jù)，通過對3種插值算法13種模型模擬值與實測值的比較和分析，優(yōu)化和篩選出實際偏離程度最小的插值模型，從而為空間插值技術(shù)在近岸海域環(huán)境監(jiān)測領(lǐng)域中的應(yīng)用提供參考和借鑒。

1 實驗部分

1.1 數(shù)據(jù)來源

研究區(qū)域為全國近岸海域，根據(jù)國家環(huán)境保護部《全國近岸海域環(huán)境質(zhì)量監(jiān)測實施方案》，按照網(wǎng)格布點法，確定299個國控站位，面積約28萬平方公里。

海水水質(zhì)數(shù)據(jù)來源于全國近岸海域環(huán)境監(jiān)測網(wǎng)2009年監(jiān)測數(shù)據(jù)，監(jiān)測項目包括水溫、懸浮物、鹽度、pH、溶解氧、化學(xué)需氧量、石油類、活性磷酸鹽、無機氮、汞、銅、鉛、鎘、非離子氨和六價鉻等。

基礎(chǔ)地理信息數(shù)據(jù)根據(jù)國家測繪局1∶25萬標(biāo)準(zhǔn)分幅地圖加工而成。

1.2 水質(zhì)評價方法

各站位的污染物濃度采用年度平均值，根據(jù)《海水水質(zhì)標(biāo)準(zhǔn)》(GB 3097—1997)［11］和《近岸海域環(huán)境監(jiān)測規(guī)范》(HJ 442—2008)［12］，按照單因子判別法進行水質(zhì)評價，即某一測點海水中任何評價指標(biāo)超過一類海水標(biāo)準(zhǔn)的，該測點水質(zhì)即為二類，超過二類海水標(biāo)準(zhǔn)的即為三類，依次類推［13］。

1.3 空間插值方法

1.3.1 反距離加權(quán)(IDW)

IDW法是一種基于相近相似原理的確定性插值方法，假設(shè)每個采樣點有一個局部影響，與未采樣點距離最近的若干個點對未采樣點值的貢獻最大，其貢獻與距離成反比，表達(dá)式為

式中:Z是評估的海水水質(zhì);Zi是站位i(i=1，2…n)的海水水質(zhì);d是評估站位和采樣站位之間的距離;p是分析定義的冪參數(shù);n代表用于評估計算用的站位數(shù)量。

影響IDW法準(zhǔn)確度的主要因子是p值［14］，當(dāng)p等于1時，是線性距離衰減插值;當(dāng)p大于1時，是非線性距離衰減插值。對較大的冪，較近的采樣站位被給定較高的權(quán)重份額;對較小的冪，權(quán)重比較均勻地分配給各采樣站位。該法的優(yōu)點是可通過權(quán)重調(diào)整空間插值等值線的結(jié)構(gòu)，缺點是如不了解研究區(qū)域內(nèi)插屬性的結(jié)構(gòu)分布特征，不合理的加權(quán)會導(dǎo)致較大的偏差，而且易受極值的影響。研究采用常用的幾個p參數(shù)(p=1，2，3，4，5，6［2，15］)的模擬計算結(jié)果進行比較。另外，臨近點的大小和數(shù)量同樣影響結(jié)果的準(zhǔn)確性，最近采樣點值取12，范圍設(shè)置 Maximum distance為1.7。

1.3.2 樣條函數(shù)法(Spline)

Spline法是用數(shù)學(xué)函數(shù)來逼近曲面的一種方法，就是按一定的光滑性要求對接起來的分段多項式。該方法計算量小，速度快，保留了局部的變化特征，有較好的視覺效果，適合于平滑的表面，一般要求有連續(xù)的一階和二階導(dǎo)數(shù);其缺點是難以對函數(shù)內(nèi)插的誤差進行估算，點稀疏時效果不好，對有限的監(jiān)測數(shù)據(jù)模擬，Spline法難以滿足內(nèi)插的精度，表達(dá)式為

式中:Z為估算值;n為站位數(shù);λ為一系列線性方程解的系數(shù);ri是估算點到i站位的距離;R(ri)和T(x，y)的表達(dá)式為

式中:τ為權(quán)重系數(shù);k0為改正后的貝塞爾函數(shù);c為常數(shù);a為線性方程的系數(shù)。分別采用Regularized和Tension 2種類型的Spline法進行插值驗證，權(quán)重系數(shù)取0.1，模型樣點個數(shù)取12。

1.3.3 克里格方法(Kriging)

Kriging法是建立在地統(tǒng)計學(xué)基礎(chǔ)上的一種插值方法，該算法以區(qū)域化變量理論為基礎(chǔ)，認(rèn)為任何在空間連續(xù)性變化的屬性是不規(guī)則的，不能用簡單的平滑數(shù)學(xué)函數(shù)進行模擬，可用隨機表面給予較恰當(dāng)?shù)拿枋?。Kriging法著重于權(quán)重系數(shù)的確定，從而對空白站位的變量值提供最好的線性無偏估計，公式為

式中:λi是站位i海水水質(zhì)Z(xi)的權(quán)重系數(shù)，表示各站位海水水質(zhì)Z(xi)對估計值Z(x0)的貢獻。為達(dá)到線性無偏估計，使估計方差最小，權(quán)重系數(shù)由方程組求得:

式中:Z(xi)是站位i的海水水質(zhì);Z(xi+h)是距離i站位h的點的水質(zhì);n代表距離為h的成對站位的數(shù)目。根據(jù)實驗半變異函數(shù)，確定合理的變異函數(shù)理論模型。分別采用 Spherical、Circular、Exponential、Guassian、Linear 5 個變異函數(shù)模型進行驗證，最近采樣點取值12，范圍設(shè)置 Maximum distance為1.7。

1.4 插值結(jié)果準(zhǔn)確度檢驗

利用交叉檢驗法［16-17］驗證各插值模型的模擬效果。其原理為每個監(jiān)測站位的海水水質(zhì)都從數(shù)據(jù)列中去掉一次，其海水水質(zhì)由周圍站位的實測值通過插值來估算，然后計算站位模擬值與實測值之間的誤差。采用模擬值與實測值的pearson相關(guān)系數(shù)(Corr)、平均誤差(ME)、平均絕對誤差(MAE)和均方差(RMSE)作為評估插值方法的標(biāo)準(zhǔn)。Corr反映驗證值與模擬值之間的吻合程度，ME反映總體模擬值誤差的大小，MAE反映模擬值可能的誤差范圍，RMSE反映模擬值的靈敏度和極值效應(yīng)，公式為

式中:C(xi，yi)為測站樣本點之間的協(xié)方差;C(xi，x0)為測站樣本點與插值點之間的協(xié)方差，μ為拉格朗日乘子。

插值數(shù)據(jù)的空間結(jié)構(gòu)特性用半變異函數(shù)γ(h)表示:

1.5 交叉檢驗的實現(xiàn)

利用 Visual Studio．NET 2005開發(fā)平臺，C++語言和ArcGIS Engine組件編寫程序，根據(jù)坐標(biāo)文件生成矢量圖層，對圖層中各站位實現(xiàn)循環(huán)，分別去除一個站位，利用剩余站位進行插值計算，生成柵格圖層，再從該圖層中提取去除點的插值結(jié)果，進行誤差統(tǒng)計。

1.6 數(shù)據(jù)統(tǒng)計、插值的實現(xiàn)

數(shù)據(jù)的統(tǒng)計、相關(guān)分析采用SPSS 13.0。插值計算采用ArcGIS Desktop 9.2，Albers雙標(biāo)緯等積圓錐投影，Krasovsky橢球體，D_Krasovsky_1940大地基準(zhǔn)面，中央經(jīng)線105°，標(biāo)準(zhǔn)緯線25°和47°。

最優(yōu)插值模型與其他模型用RI進行比較:

2 結(jié)果與討論

2.1 模擬海水水質(zhì)分布

將各插值模型的模擬結(jié)果進行柵格化處理，形成全國近岸海域水質(zhì)分布圖(圖1)。

由圖1可見，所有插值模型模擬的海水水質(zhì)污染程度都由近岸向遠(yuǎn)海遞減。這與中國近岸海域環(huán)境問題的實際情況相符，近海水體中污染物濃度主要受陸源徑流的影響，總體呈現(xiàn)近岸高、遠(yuǎn)岸低的趨勢特征;主要超標(biāo)污染物是無機氮和活性磷酸鹽，長江口、黃河口、珠江口、海河口和閩江口等及鄰近海域受陸源輸入影響，富營養(yǎng)化程度較嚴(yán)重［18］。相比較，Spline-Regularized法在廣西北部灣海域的模擬結(jié)果為污染程度嚴(yán)重，相對該海域水質(zhì)的實際水平有所高估，說明其在局部海域水質(zhì)的插值模擬效果較差。

圖1 不同插值模型的近岸海域水質(zhì)模擬結(jié)果

2.2 模擬結(jié)果誤差分析

交叉檢驗統(tǒng)計各插值模型的誤差結(jié)果(表1)。

表1 不同插值模型的交叉檢驗誤差結(jié)果

圖2 不同插值模型的交叉檢驗結(jié)果散點圖

由表1可見，13種插值模型的ME值范圍為－0.081 7～0.057 2，從小到大次序為 Spline-Regularized＜Spline-Tension＜IDW(p6)＜IDW(p5)＜IDW(p4)＜IDW(p3)＜IDW(p2)＜IDW(p1)＜Kriging-Linear＜Kriging-Spherical，Kriging-Circular＜Kriging-Guassian＜Kriging-Exponential;MAE值范圍為0.625 1～0.853 5，從小到大次序為IDW(p4)＜IDW(p5)＜IDW(p6)＜IDW(p3)＜IDW(p2)＜Spline-Tension＜IDW(p1)＜Kriging-Spherical，Kriging-Circular＜ Kriging-Guassian ＜Kriging-Linear＜Kriging-Exponential＜Spline-Regularized;RMSE值范圍為0.934 4～1.216 0，從小到大的次序為IDW(p2)＜Spline-Tension＜IDW(p3)＜IDW(p4)＜IDW(p1)＜IDW(p5)＜IDW(p6)＜ Kriging-Guassian ＜ Kriging-Spherical，Kriging-Circular＜Kriging-Linear＜Kriging-Exponential＜Spline-Regularized;Corr范圍為0.643～0.750，從大到小次序為 Spline-Tension＜IDW(p2)＜IDW(p3)＜IDW(p4)＜IDW(p1)＜IDW(p5)＜IDW(p6)＜Spline-Regularized＜Kriging-Guassian ＜ Kriging-Spherical， Kriging-Circular＜Kriging-Linear＜Kriging-Exponential，各插值模型的模擬值和實測值在置信度0.01水平上顯著性相關(guān)，說明模擬結(jié)果與實測結(jié)果很吻合，基本能反映海水水質(zhì)的平面分布規(guī)律。

比較各插值模型的模擬值與實測值(圖2)。

從線性擬合與最佳擬合效果(圖中虛線，1∶1)的偏離度來看，IDW(p=2，3，4，5，6)法和 Spline法的效果優(yōu)于其他方法，但在Spline法的結(jié)果中，個別站位海水水質(zhì)出現(xiàn)負(fù)值，偏離實際情況，說明該模型的適用性相對較差。綜合考慮各插值模型的誤差及與實際情況的吻合程度，推薦最佳插值模型是IDW(p2)，按RI值偏離最佳插值模型的次序為IDW(p2)＜Spline-Tension＜IDW(p3)＜IDW(p4)＜IDW(p1)＜IDW(p5)＜IDW(p6)＜Kriging-Guassion ＜ Kriging-Spherical，Kriging-Circular＜Kriging-Linear＜Kriging-Exponential＜Spline-Regularized。

3 結(jié)論

基于2009年全國近岸海域海水水質(zhì)監(jiān)測數(shù)據(jù)，利用交叉檢驗法，分析了常用的3種插值算法13種模型的模擬精度，驗證了其在近岸海域水質(zhì)平面分布表征中的適用性，并根據(jù)模擬效果提出最優(yōu)插值模型。

1)插值模型的模擬結(jié)果與實測結(jié)果具有顯著相關(guān)性，說明用插值模型能較好地表征海水水質(zhì)空間分布，其模擬結(jié)果可以代表實際監(jiān)測站位的水質(zhì)平面分布特征，13種空間插值模型適用于全國近岸海域海水水質(zhì)的分析和評價。

2)隨研究區(qū)域和站位設(shè)置的不同，適宜的插值模型也不相同。3種插值算法13種模型中，IDW(p2)插值模型最適宜模擬近岸海域水質(zhì)分布。但該插值模型不能代表不同區(qū)域、不同要素的適用性，今后有待于進一步研究不同污染物插值模型的精度。

［1］王業(yè)耀，李俊龍，劉方．中國近岸海域環(huán)境監(jiān)測技術(shù)路線研究［J］．中國環(huán)境監(jiān)測，2013，29(5):118-123．

［2］BurroughP A， McDonnellR A． Principlesof Geographical Information Systems．Spatial Information Systems and Geostatistics［M］．New York:Oxford University Press，1998:333．

［3］Ramirez-Herrera M，Urrutia-Fucugauchi J．Morphotectonic zones along the coast of the Pacific continentalmargin［J］．Southern Mexico Geomorphology，1999，28:237-250．

［4］Javier B P，Jose L H．Comparing the performance of two spatial interpolation methods for creating a digital bathymetric model of the Yucatan submerged platform［J］．Pan-American Journal of Aquatic Sciences，2007，2(3):247-254．

［5］向運榮，張軼男，劉俊．粵港珠江三角洲區(qū)域空氣質(zhì)量分級技術(shù)及信息發(fā)布方法研究［J］．中國環(huán)境監(jiān)測，2009，25(1):7-10．

［6］曾武濤，劉紅衛(wèi)．珠江三角洲地區(qū)空氣自動監(jiān)控網(wǎng)絡(luò)數(shù)據(jù)空間模擬適用方法［J］．中國環(huán)境監(jiān)測，2008，24(1):51-54．

［7］陳彥軍，李偉鏗，張寶春，等．基于GIS的珠三角區(qū)域空氣質(zhì)量時空演化分析模型研究［J］．中國環(huán)境監(jiān)測，2012，28(5):136-141．

［8］Klenke M，Schenke H W．A new bathymetric model for the central Fram strait［J］．Marine Geophysical Researches，2002，23:367-378．

［9］Wright D J，Bloomer S H．Bathymetry of the tonga trench and forearc:amap series［J］．Marine Geophysical Researches，2000，21:489-511．

［10］Merino M． Upwelling on the Yucatan Shelf:hydrographic evidence ［J］． JournalofMarine Systems，1997，13:101-121．

［11］GB 3097—1997 海水水質(zhì)標(biāo)準(zhǔn)［S］．

［12］HJ 442—2008 近岸海域環(huán)境監(jiān)測規(guī)范［S］．

［13］中華人民共和國環(huán)境保護部．中國近岸海域環(huán)境質(zhì)量公報［R］．北京:環(huán)境保護部，2003—2009．

［14］張慧智，史學(xué)正．中國土壤溫度的空間插值方法比較［J］．地理研究，2008，27(6):1 307-1 329．

［15］Kravchenco A，Bullock D G．A comparative study of interpolation methods for mapping soil properties［J］．Agronomy Journal，1999，13:101-121．

［16］Isaak E H．An Introduction to Applied Geostatistics［M］．New York:Oxford University Press，1989:561．

［17］Webster R，Oliver M A．Geostatistics for Environmental Science［M］．Toronto:John Wiley and Sons，2001:271．

［18］霍云龍．我國近海主要環(huán)境問題研究［D］．廈門:國家海洋局第三海洋研究所，2011．