基于137Cs本底值計算模型和地理加權回歸克里格對137Cs本底值空間分布的預測

李 豪， 文安邦， 劉 濤， 李 婷

(1.四川農(nóng)業(yè)大學資源學院， 611130，成都；2.中國科學院水利部成都山地災害與環(huán)境研究所，610041，成都)

土壤侵蝕造成土地資源退化，生產(chǎn)力下降，是危害人類生存與發(fā)展的世界性環(huán)境問題之一。深入了解土壤侵蝕規(guī)律、獲取可靠的土壤侵蝕信息是開展水土保持工作的基礎。核示蹤技術是土壤侵蝕研究的一種重要方法，具有成本低、高效可靠等優(yōu)點，其中又以137Cs示蹤技術應用最為廣泛和成功[1]。自20世紀60年代以來，在國際原子能機構(IAEA，International Atomic Energy Agency)等國際組織的大力推動下，其技術體系逐漸完善，測定的土壤侵蝕數(shù)據(jù)的準確性得到廣泛驗證和認可。

在一定區(qū)域內(nèi)，未受擾動、無侵蝕亦無沉積發(fā)生的土壤條件下137Cs的面積濃度代表該區(qū)域內(nèi)每m2上137Cs大氣沉降的累積量，稱為該區(qū)域的137Cs本底值(137Cs reference inventory，Bq/m2)。137Cs示蹤法研究土壤侵蝕的原理是分析研究區(qū)某一土壤剖面的137Cs含量與137Cs本底值之間的差異，選擇合適的計算模型，定量計算得到該處的土壤侵蝕量或沉積量，因此準確的137Cs本底值是運用137Cs示蹤技術研究土壤侵蝕的關鍵因素。目前，獲取137Cs本底值數(shù)據(jù)的主要手段是野外選擇符合沒有人為擾動、無侵蝕和沉積發(fā)生等條件的樣地采集土壤樣品，然后在實驗室測試得到；但由于近幾十年來人類活動強度大，在許多區(qū)域難以找到完全符合標準的采樣地，該方法的實現(xiàn)存在一定難度。有學者指出，由于該方法在采樣過程中普遍存在的一系列等問題，導致獲取的137Cs本底值數(shù)據(jù)的可靠性存在疑問[2-3]。

以實測數(shù)據(jù)為基礎、基于GIS的空間插值技術已廣泛用于獲取土壤、環(huán)境和氣象等要素的高分辨率空間分布數(shù)據(jù)。137Cs具有與上述要素類似的空間分布特性，空間插值技術亦適用于137Cs本底值研究；但由于現(xiàn)階段我國137Cs本底值的實測樣點數(shù)量偏少(近30年來，全國137Cs本底值研究樣點僅有100余個[4])，且空間分布不均勻，并不適合通過137Cs本底值實測數(shù)據(jù)直接進行空間插值。采用GIS的空間插值技術獲取區(qū)域137Cs本底值資料，須另辟途徑。

筆者以長江上游的四川省為研究對象，嘗試采用137Cs本底值理論計算模型、主成分分析與地理加權回歸克里格插值相結合的技術方法，對137Cs本底值空間分布進行研究，探討通過氣象站點相對豐富的實測數(shù)據(jù)、基于本底值理論計算模型和GIS空間插值技術獲取區(qū)域137Cs本底值數(shù)據(jù)的合適方法，以期獲取更高精度的區(qū)域137Cs本底值資料，為在該區(qū)開展土壤侵蝕的核示蹤研究、進行水土流失治理等工作提供有效的數(shù)據(jù)支持。

1 研究區(qū)概況

四川省位于我國西南，地處長江上游，E 92°21′～108°12′，N 26°03′～34°19′，區(qū)域面積48.5萬km2。本區(qū)處于第1級青藏高原和第2級長江中下游平原的過渡帶，區(qū)內(nèi)海拔高差懸殊，地貌復雜多樣，區(qū)域以山地地貌為主，全省可分為東部四川盆地、川西北高原和川西南山地等3大部分。區(qū)內(nèi)季風氣候明顯，氣候類型多樣，差異顯著：東部為亞熱帶濕潤、半濕潤氣候，雨熱同季，水熱條件好，年降水量1 000～1 200 mm，但多云霧、日照時間較短；西部為高山高原高寒氣候，氣候垂直變化明顯，總體上以寒溫帶氣候為主，冬寒夏涼，水熱不足，年降水量500～900 mm。四川省土壤類型豐富，區(qū)域分異明顯，據(jù)第2次土壤普查，全省土壤類型共有25個土類、66個亞類，以水稻土、紫色土、黃壤等為主，本區(qū)是我國紫色土分布最集中的地區(qū)。

四川省是我國土壤侵蝕較嚴重的省區(qū)之一，土壤侵蝕分布廣，面積大，危害嚴重。本區(qū)土壤侵蝕面積和數(shù)量約占整個長江上游總量的70%，長江上游的泥沙主要來源于地處金沙江和嘉陵江中下游的四川盆地丘陵地區(qū)。根據(jù)第1次全國水利普查公報數(shù)據(jù)，全省水蝕面積達11.4萬km2，占幅員面積的23.54%，其中中度以上土壤侵蝕面積占比達到了57.63%[5]。

圖1 研究區(qū)域及氣象站點分布圖Fig.1 Spatial distribution of meteorological stations in the studied area

2 數(shù)據(jù)來源與研究方法

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

筆者采用的四川省147個氣象站點的經(jīng)緯度、年均降水量(觀測精度為0.1 mm)等數(shù)據(jù)來自國家氣象信息中心提供的“中國基本、基準和一般地面氣象觀測站1981—2010年累年值年值數(shù)據(jù)集”。各氣象站點的空間分布如圖1所示。雖然該數(shù)據(jù)集與本文的研究對象在時間尺度上并不完全相同;但是，該數(shù)據(jù)集具有較長的時間序列(30年)，均為地面實測數(shù)據(jù)，可信度高，且站點數(shù)量多，能保證插值的精度，滿足研究要求。在參照類似研究的技術方法[4]后，我們采用該數(shù)據(jù)集作為本研究的降水數(shù)據(jù)。

2.2 數(shù)據(jù)處理

D.E.Walling等[6]建立的137Cs本底值全球空間分布模型(以下簡稱W&H模型)是基于核爆產(chǎn)物大氣沉降數(shù)據(jù)、建立在認識137Cs沉降過程的基礎上的理論模型，物理意義明確，目前應用最為廣泛。本文各氣象站點的137Cs本底值模型計算值(以下簡稱137Cs本底計算值)通過該模型計算得到，本文涉及的137Cs本底值數(shù)據(jù)均衰變校正至2016年1月1日。

W&H模型[6]按照10°緯度帶劃分計算帶，對各帶分別計算137Cs本底值數(shù)據(jù)，而四川省地處N 26°～34°之間，跨越了N 20°～30°和30°～40° 2個計算帶。為提高預測結果的精度，筆者對四川省進行了細分，以N30°為界劃分為2個研究區(qū)域(圖1)。其中，低于N30°區(qū)域(以下簡稱L30區(qū)域)包含63個站點，高于N 30°區(qū)域(以下簡稱H30區(qū)域)包含84個站點，對這2個研究區(qū)域分別預測其137Cs本底計算值的空間分布。

降水、地理位置及進行大氣核試驗的場所等是影響我國137Cs本底值分布的主要因素[7]。而我國的大氣核試驗均在新疆羅布泊地區(qū)進行，對本研究區(qū)域的137Cs本底值分布基本沒有影響；因此，筆者選擇年降水量(代表降水)、經(jīng)度和緯度(代表地理位置)等3個影響因素作為輔助變量，用于回歸分析和空間插值。筆者采用的技術方法主要包括主成分分析和地理加權回歸克里格插值。

2.3 主成分分析

主成分分析(principal component analysis, PCA)是一種應用廣泛的統(tǒng)計分析方法，主要用于降維和消除變量間的相關性。其通過正交變換將1組存在相關性的變量轉換為1組數(shù)量較少的、線性不相關的變量(稱之為主成分)，在盡可能多地保留原來變量的信息的基礎上，可減少變量的數(shù)量，有效降低分析問題的復雜性[8]。

2.4 地理加權回歸克里格插值

克里格插值法是一種基于地統(tǒng)計學的經(jīng)典空間插值方法，目前已在多個領域得到廣泛應用。有學者提出，通過引入與目標變量具有較好相關性的影響因素變量、構建回歸克里格模型(regression Kriging, RK)，可以有效提高克里格法的插值精度[9-11]。RK法的基本思想是建立目標變量和影響因素之間的回歸方程，分離出趨勢項，然后對回歸殘差(即隨機因素)采用傳統(tǒng)的普通克里格插值法(ordinary Kriging, OK)進行插值，最后對回歸預測的趨勢項和回歸殘差的估計值進行求和，從而得到目標變量的預測值。RK法可用以下式[12]表示：

(1)

式中:y(pi)為位置pi處的目標變量；βj,(pi)為回歸系數(shù)；xj,(pi)為位置pi處的第j個自變量；ε(pi)為位置pi處的回歸殘差；n為自變量的個數(shù)。

傳統(tǒng)的RK法中，回歸分析部分大都采用基于最小二乘法(ordinary least squares, OLS)的全局回歸(global regression, GR)，各回歸系數(shù)為常數(shù)，即為全局回歸克里格法(global regression Kriging, GRK)。由于沒有考慮變量的空間非平穩(wěn)性，因此在處理土壤、降水等具有較強空間非平穩(wěn)性的變量時，該方法的插值精度仍受到一定的限制。

地理加權回歸克里格法(geographically weighted regression Kriging, GWRK)法對GRK法進行了改進，是將GRK法中基于OLS的全局回歸替換為地理加權回歸的局部加權回歸(local regression, LR)，其余步驟、方法保持不變的一種空間插值方法。GWRK法的核心是地理加權回歸(geographically weighted regression, GWR)。GWR法對全局回歸模型進行了改進，通過構建空間權重函數(shù)，對空間非平穩(wěn)性進行了量化。在GWR中，回歸系數(shù)不再是常數(shù)，而是空間位置的函數(shù)，位置(ui,vi)處的回歸系數(shù)β通過下式[12]計算得到：

(ui,vi)=[XTW(ui,vi)X]-1XTW(ui,vi)Y。

(2)

式中:β為回歸系數(shù)矩陣；X為自變量設計矩陣(XT為其轉置矩陣)；W(ui，vi)為空間權重矩陣，由空間權重函數(shù)W(i)求得，其作用是定量衡量領域內(nèi)不同空間位置k(k=1，2，…，n)樣點的觀測值對于回歸點(ui，vi)回歸系數(shù)估計的影響程度[12](n為領域樣點數(shù))；Y為因變量矩陣。

已有研究表明，GWRK的插值結果能揭示被空間非平穩(wěn)性所掩蓋的一些局部變化，更真實地反映目標變量的空間變異情況，使插值精度得到進一步的提高[13]。基于上述分析，筆者采用GWRK法進行四川省137Cs本底計算值的插值研究。

2.5 模型精度的評價

從147個站點數(shù)據(jù)中隨機選取30個站點構建驗證數(shù)據(jù)集。采用其余117個站點數(shù)據(jù)建立137Cs本底值插值模型，并求取30個驗證站點的137Cs本底值預測值，通過比較驗證站點的理論計算值與插值預測值之間的平均誤差(mean error, ME)、平均絕對誤差(mean absolute error，MAE)、均方根誤差(root mean square error, RMSE)[13]、平均絕對百分誤差(mean absolute relative error, MARE)和均方根相對誤差(root mean square relative error, RMSRE)等指標評估不同插值模型的插值精度。其中，MARE和RMSRE定義如下：

(3)

(4)

式中:ya,i為位置i處的137Cs本底計算值，Bq/m2；ye,i為位置i處的137Cs本底預測值，Bq/m2；n為驗證數(shù)據(jù)集的站點個數(shù)。

2.6 研究軟件平臺

常規(guī)統(tǒng)計分析和PCA采用SPSS 23完成，基于OLS的GR和GWR回歸分析由GWR4完成，半方差函數(shù)計算通過GS+9進行，克里格插值和各空間分布圖的繪制由ArcGIS 10.4完成。各變量空間分布圖的分辨率均為1 km×1 km。

3 結果與分析

3.1 137Cs本底值描述性統(tǒng)計分析

對研究區(qū)域各站點的137Cs本底計算值指標進行了描述性統(tǒng)計分析(表1)。分析結果表明，在2個研究區(qū)域各站點的137Cs本底計算值在數(shù)量上存在一定的差異，H30區(qū)域和L30區(qū)域的平均值分別為1 167.42和616.39 Bq/m2，而137Cs本底計算值的變異系數(shù)分別為8.65%和7.48%，在分異程度上基本一致。2個區(qū)域的137Cs本底計算值數(shù)據(jù)均通過了J-B(Jarque-Bera)檢驗(P>0.05)，符合正態(tài)分布。

表1 研究區(qū)域站點的137Cs本底計算值描述性統(tǒng)計分析和正態(tài)性檢驗

3.2 影響因素主成分分析

表2為計算的研究區(qū)域137Cs本底計算值與各影響因素間的相關系數(shù)，結果顯示：137Cs本底計算值與各影響因素均存在較高的相關，與經(jīng)度、降水量呈極顯著正相關，與緯度呈極顯著負相關。另一方面，各影響因素間也存在著不同程度的相互作用，如降水量與經(jīng)度間達到了極顯著正相關。因此，為了消除變量間的信息冗余性和多重共線性，更好地反映各變量間的真實關系，首先對降水量、經(jīng)度和緯度等3個影響因素進行了主成分分析。

對各影響因素變量進行主成分分析后，根據(jù)方差累計貢獻率≥85%的原則提取了主成分(表3)，結果表明：對于2個研究區(qū)域，前2項主成分(F1、F2)的方差累計貢獻率均超過90%，分別達到了95.39%和94.66%，表明前2項主成分已經(jīng)表達了原始數(shù)據(jù)中絕大部分信息，采用前2項主成分即可較好地替代原數(shù)據(jù)信息；因此，后續(xù)的分析中將降水量、經(jīng)度和緯度等3個原始影響因素指標轉換為2個主成分變量F1和F2，2個主成分的空間分布如圖2所示。

表2 137Cs本底值與各影響因素的相關系數(shù)

注：同一變量行上方數(shù)據(jù)表示H30區(qū)域的結果，下方數(shù)據(jù)表示L30區(qū)域的結果；*表示在0.05水平上具有顯著性，** 表示在0.01水平上具有顯著性,下同. Notes: The data in the top of the same row refer to the results of H30, the data in the bottom of the same row refer to the results of L30. * refers to significance at 0.05 level, ** refers to the significance at 0.01 level. The same below.

表3 主成分的特征值與方差貢獻率

Notes:F1refers to for the first principal component,F2refers to for the second principal component, andF3refers to for the third principal component. The same below.

(a) H30區(qū)域F1變量 F1 in H30；(b) H30區(qū)域F2變量 F2 in H30；(c) L30區(qū)域F1變量 F1 in L30；(d) L30區(qū)域F2變量 F2 in L30圖2 研究區(qū)域主成分變量空間分布圖Fig.2 Spatial distribution of principal components in the studied area

3.3 回歸分析

以各站點的2個主成分值為自變量，137Cs本底計算值為因變量，分別采用基于OLS的GR和GWR2種模型進行了回歸分析。2種回歸分析的結果分別如表4和表5所示。

表4 基于OLS的GR回歸計算結果

Notes: GR stands for global regression. OLS stands for ordinary least square. VIF stands for variance inflation factor. The same below.

表5 GWR回歸計算結果

Notes:GWR stands for geographically weighted regression. The same below.

方差膨脹因子(variance inflation factor，VIF)是衡量自變量間是否存在多重共線性的主要指標。由表4的計算結果來看，對原始變量進行主成分分析后，回歸分析各自變量的VIF值遠小于10，這說明經(jīng)過PCA處理后，各自變量間基本不存在多重共線性問題，回歸分析的結果是可靠的。從GWR模型的計算結果來看，各解釋變量的回歸系數(shù)值不再是常數(shù)，如H30區(qū)域的F1回歸系數(shù)的取值范圍介于-57.19與-34.18間，變異系數(shù)達到了33.46%，這說明137Cs本底計算值與各解釋變量間的相互關系具有顯著的空間非平穩(wěn)性。

對2種模型的回歸效果進行了比較。表6的分析結果表明：GWR模型的各項診斷指標均明顯優(yōu)于基于OLS的GR模型。如在H30區(qū)域，GWR模型的RSS與GR模型相比，大幅減少了57.33%；AICc值也從GR模型的674.61顯著下降到GWR模型的639.37。這說明基于局部回歸的GWR模型能夠更好的探測和解釋變量間的空間非平穩(wěn)性，相對于GR模型，具有更好的回歸效果。

表6 GR和GWR模型的診斷指標

3.4 回歸殘差的克里格插值

正態(tài)性檢驗的結果顯示，GR模型在2個研究區(qū)域、GWR模型在L30區(qū)域的殘差均通過了J-B檢驗，為正態(tài)分布，而GWR模型在H30區(qū)域的回歸殘差不符合正態(tài)分布；因此，首先對該殘差數(shù)據(jù)進行5/7次方轉換，使其滿足進行克里格插值的條件。

在GS+9軟件中，對2種模型的回歸殘差(包括H30區(qū)域轉換后的回歸殘差)進行半方差函數(shù)擬合。從表7的擬合結果可以看出，各半方差函數(shù)的R2均達到0.80，取得了較好的擬合效果。此外，各殘差的半方差函數(shù)的塊基比(即C0/(C0+C))大都小于25%，說明各殘差存在較強的空間自相關性，適合使用克里格法進行空間插值。

表7 GR和GWR模型回歸殘差的半方差函數(shù)擬合結果

根據(jù)上述半方差函數(shù)的擬合結果，采用OK法分別對2種模型的回歸殘差進行了空間插值。根據(jù)2個主成分變量及其回歸系數(shù)、殘差插值等數(shù)據(jù)的計算結果，分別獲得了基于GRK模型和GWRK模型的137Cs本底計算值空間分布圖(圖3)。

3.5 模型預測精度評價

作為對照，同時采用OK法對各站點的137Cs本底計算值進行了空間插值。通過構建的驗證數(shù)據(jù)集分別評價了OK、GRK和GWRK等模型的預測精度。各評價指標的值越接近0，說明相應模型的預測精度越高。由表8的分析結果可以看出，在各研究區(qū)域，模型預測精度均為GWRK>GRK>OK。綜合考慮了降水、空間位置等多個影響因素的GRK模型和GWRK模型對137Cs本底值的預測精度較OK法均有較大提高，而進一步探究了各影響因素空間非平穩(wěn)性的GWRK模型的預測效果是這3種模型中最好的。

(a) H30區(qū)域GWRK結果 Result of GWRK model in H30；(b) H30區(qū)域GRK結果 Result of GRK model in H30；(c) L30區(qū)域GWRK結果 Result of GWRK model in L30；(d) L30區(qū)域GRK結果 Result of GRK model in L30圖3 四川省137Cs本底計算值空間分布Fig.3 Spatial distribution of 137Cs reference inventory of Sichuan province

表8 不同預測方法的預測精度對比

Notes: ME stand for mean error, MAE stand for mean absolute error, RMSE stand for root mean square error, MARE stand for mean absolute relative error, and RMSRE stand for root mean square relative error.

同時收集了已發(fā)表文獻中137Cs本底值的實測數(shù)據(jù)，與GWRK模型的預測結果進行了對比分析。從表9的對比結果可以看出，GWRK模型的預測值普遍略低于實測值，差異幅度大致為10%?？紤]到測定137Cs活度的γ能譜儀的測試誤差普遍在5%～10%之間，我們?nèi)∑渖舷?0%作為實測本底值和預測值是否存在差異的標準。由此可見，雖然預測值略低于實測值，但差異不大，實測值與本文的預測結果是基本吻合的。

3.6 四川省137Cs本底值的空間分布

由圖3的插值結果可以看出，通過GRK和GWRK法得到的四川省137Cs本底計算值空間分布的總趨勢基本一致，均呈現(xiàn)北高南低、東高西低及高原低盆地高、盆周山地高盆中丘陵低的趨勢，該預測結果與前人對全球137Cs沉降規(guī)律的認識基本一致[7]。一些局部的、細節(jié)性的137Cs本底計算值空間分布特征在插值結果中也有所反映：在四川盆地西部邊緣出現(xiàn)了1個局部137Cs本底計算值峰值區(qū)，大致位于雅安市、樂山一線，137Cs本底計算值達到了1 300～1 400 Bq/m2。出現(xiàn)這一現(xiàn)象主要是由于該區(qū)的“華西雨屏”[19]現(xiàn)象造成降水量較高，而137Cs主要來源于降水帶來的濕沉降所致。此外，由于W&H計算模型的局限——模型中按照10°的緯度帶劃分計算區(qū)域，造成北緯30°兩側區(qū)域的137Cs本底計算值存在“斷崖式”突變。在資料精度允許的條件下，今后可以通過細分基本區(qū)域，改進模型，提高模擬結果的精度和可信度[7]。

表9 137Cs本底值實測值與GWRK法預測值比較

4 結論

以地處長江上游、土壤侵蝕特征具有典型性的四川省為研究對象，基于147個氣象站點的實測數(shù)據(jù)以及W&H模型的137Cs本底值分布模型，采用GWRK法進行空間插值，獲得四川省137Cs本底計算值的空間分布數(shù)據(jù)。初步研究結果表明：

1)應用PCA法將經(jīng)緯度、降水量等3個影響因素轉換為2個主成分變量，主成分變量能表達原始數(shù)據(jù)中90%以上的信息，同時有效的消除了變量間的多重共線性，為后續(xù)的空間插值分析奠定了基礎。

2)根據(jù)PCA的分析結果，通過基于OLS的GRK和GWRK模型對四川省137Cs本底計算值空間分布進行了預測，并評價了不同模型的預測精度?；貧w分析的結果表明，基于局部回歸的GWRK模型比GRK模型能更好地解釋137Cs本底計算值與各影響因素間相互關系的空間非平穩(wěn)性，具有更好的擬合效果；從最終的插值結果來看，GWRK模型綜合考慮了降水、空間位置等影響因素，并進一步探究了各影響因素對137Cs本底值影響的空間非平穩(wěn)性，能更好地揭示影響因素和137Cs本底值間的相互關系，其預測精度最高，且與實測結果較吻合。

3)以實測數(shù)據(jù)為基礎、基于GIS空間插值技術獲取的四川省137Cs本底計算值數(shù)據(jù)具有較高空間分辨率，并能比較準確地表達局部區(qū)域137Cs本底值發(fā)生突變的細節(jié)信息。筆者提出的PCA+GWRK插值法可以有效的揭示137Cs本底值的空間分布規(guī)律及不同因素對其的影響作用，對運用137Cs示蹤技術開展侵蝕泥沙研究具有一定的借鑒作用，可以為水土流失治理、土壤退化防治及水土保持工程等工作提供有效的技術支持。