基于回歸克里格法的土壤鹽分采樣點(diǎn)布局優(yōu)化

徐 英 謝若禹 沈麗佳 馮紹元

(揚(yáng)州大學(xué)水利科學(xué)與工程學(xué)院， 揚(yáng)州 225009)

0 引言

土壤鹽漬化是干旱、半干旱區(qū)農(nóng)業(yè)和生態(tài)環(huán)境可持續(xù)發(fā)展的主要制約因素之一。掌握土壤鹽漬化時(shí)空分布規(guī)律是防治土壤鹽漬化的前提條件，而繪制土壤特性分布圖是表達(dá)土壤特性空間分布最有效的途徑之一[1-2]。目前，獲取土壤鹽分空間分布圖的主要途徑有3種：基于大量野外采樣數(shù)據(jù)，借助空間插值方法繪制[3]；通過遙感技術(shù)反演鹽分分布圖[4]；利用遙感數(shù)據(jù)以及過往資料與影像等作為輔助變量(亦稱為協(xié)變量)，利用空間插值方法繪制[5-6]。無論哪種途徑，都需要一定量的土壤鹽分實(shí)測數(shù)據(jù)用以建?；蝌?yàn)證。有研究表明，土壤特性分布圖的全部誤差由采樣誤差和數(shù)據(jù)分析誤差組成，其中采樣誤差占全部誤差的90%[7]，而80%～85%的采樣誤差源于樣本的代表性不足[8]?？梢?，優(yōu)化采樣設(shè)計(jì)對(duì)提高土壤鹽分分布圖繪制精度或降低野外采樣和室內(nèi)分析成本具有重要意義。

采樣設(shè)計(jì)包括確定采樣數(shù)量和采樣點(diǎn)分布。采樣數(shù)量通常由預(yù)算經(jīng)費(fèi)、允許采樣時(shí)間以及可用勞力決定，為保證一定的繪圖精度，優(yōu)化采樣點(diǎn)空間布局顯得尤為重要[7]。根據(jù)目的不同，采樣設(shè)計(jì)還要考慮不可采樣區(qū)域(如本研究區(qū)的居民區(qū))、如何利用以往采樣數(shù)據(jù)以及其他輔助信息等約束條件，在繪圖精度和采樣成本之間取得平衡[9]。在具有土壤實(shí)測數(shù)據(jù)或先驗(yàn)方差的情況下，研究者們利用克里格估計(jì)方差只依賴于采樣點(diǎn)空間布局、與實(shí)際取值無關(guān)的特點(diǎn)，以平均普通克里格估計(jì)方差(Ordinary Kriging estimated variance，OKV)最小(采樣數(shù)固定)為準(zhǔn)則優(yōu)化土壤采樣布局[10-12]。普通克里格法(Ordinary Kriging，OK)采樣設(shè)計(jì)能很好地表達(dá)出研究變量的地理空間變異[11]，變量均質(zhì)性越強(qiáng)，采樣也越趨于規(guī)則性空間分布，有利于表達(dá)變程范圍內(nèi)主變量的變異，但因不能很好地覆蓋受環(huán)境因素(包括人類活動(dòng)相關(guān)的因素)約束或影響的屬性空間，故對(duì)特異值極其不敏感[12]。與此相反，考慮土壤特性與各種協(xié)變量(包括環(huán)境變量、地理坐標(biāo)、遙感數(shù)據(jù)等)之間的線性、非線性關(guān)系，以平均估計(jì)方差最小為準(zhǔn)則優(yōu)化采樣點(diǎn)布局，可得到更趨于聚集性空間分布的樣本[12-13]。這樣的采樣點(diǎn)布局雖然對(duì)地理空間覆蓋性較差，但適合表征變量在屬性空間(由各協(xié)變量所代表的屬性空間)的變異，有助于包含更多極值點(diǎn)，進(jìn)而降低外延誤差[14-15]，對(duì)于應(yīng)用遙感數(shù)據(jù)優(yōu)化采樣過程非常有效[16]。因土壤特性同時(shí)受成土過程和環(huán)境因素的影響，以繪制土壤特性分布圖為目的的采樣布局，不僅需要良好的覆蓋地理空間，同時(shí)需要在屬性空間有良好的分布。為此，研究者們引入各種環(huán)境因子作為協(xié)變量，用協(xié)變量與土壤特性之間的回歸關(guān)系表達(dá)屬性空間分布，用殘差表達(dá)地理空間分布，試圖在二者之間尋找平衡。BRUS等[13]將環(huán)境因子作為協(xié)變量，以平均回歸克里格估計(jì)方差(Regression Kriging estimated variance，RKV)最小為準(zhǔn)則，通過空間模擬退火(Spatial simulated annealing，SSA)算法優(yōu)化了荷蘭Noord-Brabant省某區(qū)域地下水位動(dòng)態(tài)監(jiān)測網(wǎng)，結(jié)果表明回歸克里格法(Regression Kriging，RK)得到的監(jiān)測點(diǎn)布局能在屬性空間和地理空間之間取得很好的平衡，估值精度也比OK法和多元線性回歸法更高。SZATMARI等[9]以地面高程和土地利用類型為輔助變量，在匈牙利中部Mez?f?ld地區(qū)優(yōu)化了土壤有機(jī)質(zhì)采樣點(diǎn)布局，結(jié)果也表明RK法優(yōu)化得到的采樣點(diǎn)布局能很好地覆蓋屬性空間和地理空間，同時(shí)SSA允許優(yōu)化中考慮不可采樣區(qū)域等約束條件。盡管RK法被廣泛地用于優(yōu)化各種土壤特性的采樣點(diǎn)布局，但針對(duì)土壤鹽分采樣點(diǎn)布局優(yōu)化的研究很少。

受降雨、蒸發(fā)、地下水位、高程、土質(zhì)等自然因素與灌溉、排水、荒地分布、作物布局等人類活動(dòng)的共同影響，內(nèi)蒙古河套灌區(qū)土壤鹽分時(shí)空分布格局在不斷發(fā)生變化。為了更好地進(jìn)行鹽漬化防治和鹽漬化土壤管理，建立相對(duì)穩(wěn)定的土壤鹽分監(jiān)測網(wǎng)絡(luò)非常必要。為滿足各種插值方法的需求，土壤鹽分采樣在保證插值精度、降低成本的前提下，其樣點(diǎn)分布應(yīng)在屬性空間和地理空間均有較好的代表性。為此，本文擬引入?yún)f(xié)變量，用RK法(即平均RKV最小化為優(yōu)化準(zhǔn)則，以SSA為尋優(yōu)方法)對(duì)永濟(jì)灌域一研究區(qū)進(jìn)行土壤鹽分采樣布局優(yōu)化，以期為河套灌區(qū)土壤鹽漬化監(jiān)測提供參考。

1 材料與方法

1.1 研究區(qū)概況與數(shù)據(jù)采集

(1)研究區(qū)概況

內(nèi)蒙古河套灌區(qū)永濟(jì)灌域?qū)儆跍貛Т箨懶詺夂?，干旱少雨，蒸發(fā)強(qiáng)烈，因大量引黃灌溉造成地下水位埋深較淺，土地鹽堿化嚴(yán)重。

研究區(qū)(東經(jīng)108°49′～108°52′，北緯40°48′～40°52′)位于內(nèi)蒙古河套灌區(qū)永濟(jì)灌域合濟(jì)渠控制范圍內(nèi)，如圖1所示，研究區(qū)南北約7 km，東西約4 km，面積約為1 574 hm2。研究區(qū)平均地下水埋深為1.7 m；高程在1 040～1 041 m之間且地勢平緩(自然坡降平均為1/5 000)，區(qū)內(nèi)土壤以粉砂質(zhì)壤土為主，部分區(qū)域?yàn)樯百|(zhì)壤土；研究區(qū)多年平均降水量125.2 mm，多年蒸發(fā)量在1 999～2 346 mm之間；平均氣溫為8.3℃，年日照時(shí)數(shù)3 100 h左右；主要作物有小麥、玉米、葵花、瓜菜等。

圖1 研究區(qū)位置及樣點(diǎn)分布Fig.1 Distributions of study area and sampling point

(2)土樣采集

采用嵌套式采樣，即大尺度嵌套小尺度，大尺度采樣范圍(下稱大尺度區(qū))即為整個(gè)研究區(qū)，小尺度采樣區(qū)(下稱小尺度區(qū))位于研究區(qū)的中部偏南，面積約為20.25 hm2。

在谷歌地圖上分別勾勒出大尺度區(qū)和小尺度區(qū)，以400 m×400 m規(guī)則網(wǎng)格在大尺度區(qū)布置89個(gè)采樣點(diǎn)(圖1b)，以50 m×50 m規(guī)則網(wǎng)格在小尺度區(qū)布置100個(gè)點(diǎn)位(圖1c)，注意點(diǎn)位避開村莊和道路。所有采樣點(diǎn)采樣深度均為0～10 cm、10～20 cm、20～40 cm(每個(gè)采樣點(diǎn)重復(fù)2次)，以0～40 cm土層平均鹽分進(jìn)行研究。采樣時(shí)間為2020年9月1日(秋澆前)。

(3)數(shù)據(jù)測定

鹽分測定：將土樣自然風(fēng)干、研磨、過2 mm篩，按土水質(zhì)量比1∶5進(jìn)行土壤浸提、振蕩、過濾，土壤溶液的電導(dǎo)率測定使用DDSJ-308A型電導(dǎo)率儀，并以經(jīng)驗(yàn)公式換算成土壤全鹽量，公式為

M=2.882EC1∶5+0.183

式中M——全鹽量，g/kg

EC1∶5——土水質(zhì)量比1∶5電導(dǎo)率，mS/cm

與荒地距離測定：在谷歌地圖中勾勒出研究區(qū)荒地的邊界，將其導(dǎo)入ArcGIS 10.0軟件中計(jì)算各網(wǎng)格結(jié)點(diǎn)與荒地邊界的最短距離。

與溝道距離測定：在谷歌地圖中畫出研究區(qū)溝道，將其導(dǎo)入ArcGIS 10.0軟件中計(jì)算各網(wǎng)格結(jié)點(diǎn)與溝道的最短距離。

質(zhì)地測定：使用激光粒度分析儀(馬爾文Mastersizer 2000型)測定土壤顆粒組成，土粒分級(jí)采用美國制分級(jí)標(biāo)準(zhǔn)(粒徑小于0.002 mm為粘粒；粒徑[0.002 mm，0.05 mm]為粉粒；粒徑(0.05 mm，2 mm)為砂粒)。

1.2 研究方法

本文基于RK法計(jì)算研究區(qū)內(nèi)平均估計(jì)方差，并使用SSA尋找平均估計(jì)方差全局最小的采樣點(diǎn)布局，最后對(duì)尋找到的采樣點(diǎn)布局進(jìn)行評(píng)價(jià)。

1.2.1回歸克里格法(RK)

RK是一種混合插值方法，由線性回歸和OK法結(jié)合而成，即通過建立協(xié)變量(通常為容易獲得的環(huán)境變量)和主變量(所研究變量)之間的回歸方程，分離趨勢項(xiàng)，并對(duì)殘差進(jìn)行OK插值，最后將回歸預(yù)測的趨勢項(xiàng)和殘差的普通克里格估計(jì)值進(jìn)行空間疊加，從而得到未采樣點(diǎn)x0處主變量的估值。因此，RK法中趨勢項(xiàng)表達(dá)了協(xié)變量對(duì)主變量空間變異的影響，殘差估值反映了去除協(xié)變量影響后的主變量空間變異。此時(shí)，x0處的估計(jì)方差計(jì)算式為[17]

(1)

c(0)——?dú)埐畹姆讲?/p>

c0——采樣點(diǎn)和未采樣點(diǎn)殘差之間的協(xié)方差向量

C——采樣點(diǎn)殘差的方差-協(xié)方差矩陣

q——采樣點(diǎn)的協(xié)變量矩陣

q0——未采樣點(diǎn)的協(xié)變量向量

式(1)整合了回歸方差和殘差的普通克里格估計(jì)方差，因此，RKV既反映了協(xié)變量代表的屬性空間插值精度，又反映了地理空間的插值精度。采樣點(diǎn)布局優(yōu)化問題可以轉(zhuǎn)換為在采樣空間內(nèi)尋求平均RKV最小的采樣集。

為了評(píng)價(jià)RK法在優(yōu)化土壤鹽分采樣布局方面的表現(xiàn)，以O(shè)K法(即尋求OK估計(jì)方差平均值最小)作為對(duì)照方法。普通克里格估計(jì)方差計(jì)算見文獻(xiàn)[18]。

1.2.2空間模擬退火算法(SSA)

SSA是目前最常用的空間尋優(yōu)算法之一，其運(yùn)行原理是模擬固體退火過程：固體受熱熔解為液體后冷卻時(shí)，粒子運(yùn)動(dòng)隨溫度下降而逐漸有序，最后圍繞晶體格點(diǎn)做微小振動(dòng)。退火過程中，系統(tǒng)能量隨溫度降低趨于最小值。

本文將研究區(qū)按照50 m×50 m離散成6 352個(gè)網(wǎng)格節(jié)點(diǎn)，去除村莊內(nèi)的節(jié)點(diǎn)剩余5 187個(gè)節(jié)點(diǎn)，以這些節(jié)點(diǎn)為潛在采樣點(diǎn)進(jìn)行優(yōu)化。

使用SSA將從一個(gè)初始解(采樣方案)出發(fā)(初始解在已知采樣點(diǎn)中隨機(jī)抽取生成)，經(jīng)過大量擾動(dòng)和概率接受后，最終求得滿足目標(biāo)函數(shù)的全局最優(yōu)解。

(1)溫度參數(shù)

溫度參數(shù)決定了算法的循環(huán)次數(shù)，本研究中初始溫度T0=100℃，結(jié)束溫度Te=0.1℃，衰減系數(shù)α為0.995。

(2)目標(biāo)函數(shù)

目標(biāo)函數(shù)代表退火過程中的系統(tǒng)能量，SSA的目標(biāo)就是找到目標(biāo)函數(shù)全局最小的解。本研究目標(biāo)函數(shù)設(shè)為研究區(qū)內(nèi)平均RKV，其公式為

(2)

式中σ2——研究區(qū)內(nèi)潛在采樣點(diǎn)處的RKV

xe,j——研究區(qū)域內(nèi)采樣點(diǎn)的點(diǎn)位

S——當(dāng)前采樣點(diǎn)組合

ne——潛在采樣點(diǎn)數(shù)

以平均RKV最小為目標(biāo)函數(shù)可以同時(shí)在地理空間和特征空間內(nèi)對(duì)采樣點(diǎn)布局進(jìn)行優(yōu)化，同時(shí)RKV的計(jì)算獨(dú)立于采樣點(diǎn)的觀測值，可以在實(shí)際采樣前計(jì)算[17]。

(3)Metropoils準(zhǔn)則

Metropoils準(zhǔn)則定義了擾動(dòng)后新解被接受的概率Pc，計(jì)算見文獻(xiàn)[10]，本文中c初始值為0.025，衰減系數(shù)為0.995。Metropoils準(zhǔn)則使得SSA在一開始能以一定概率接受劣化解，隨著算法運(yùn)行，接受劣化解的概率逐漸較小，這樣避免優(yōu)化結(jié)果陷入局部最優(yōu)。

(4)擾動(dòng)

擾動(dòng)的作用是在當(dāng)前解的基礎(chǔ)上產(chǎn)生一個(gè)新解，擾動(dòng)的原理是在當(dāng)前解中任取一采樣點(diǎn)，以該點(diǎn)為圓心，擾動(dòng)范圍h為半徑畫圓，將原采樣點(diǎn)替換為圓內(nèi)隨機(jī)一個(gè)潛在采樣點(diǎn)。h初始值取研究區(qū)最大長度，本文h初始值為6 000 m，h隨溫度衰減而衰減，衰減系數(shù)為0.995。

1.2.3合理采樣數(shù)計(jì)算

根據(jù)經(jīng)典統(tǒng)計(jì)學(xué)，使用Cochran公式確定土壤鹽分的合理采樣數(shù)，計(jì)算公式為

(3)

式中n′——合理采樣數(shù)

ta——置信度為a時(shí)的概率度

σ——樣本標(biāo)準(zhǔn)差

K——采樣允許誤差，%

1.2.4采樣布局評(píng)價(jià)方法

為評(píng)價(jià)RK法在土壤鹽分采樣布局優(yōu)化中的表現(xiàn)，本文除了通過與OK法比較，還通過G函數(shù)、F函數(shù)和K-S檢驗(yàn)來評(píng)價(jià)優(yōu)化結(jié)果的代表性。

G函數(shù)即為最鄰近距離分布函數(shù)(Nearest neighbour distribution function)，F(xiàn)函數(shù)為間隔距離分布函數(shù)(Empty space distribution function)。G函數(shù)與F函數(shù)常一起分析樣點(diǎn)的空間分布情況，檢驗(yàn)它們對(duì)地理空間的覆蓋程度或代表性。G函數(shù)描述了從任意采樣點(diǎn)到其最近采樣點(diǎn)的距離分布。將優(yōu)化采樣集(本文指優(yōu)化后的采樣數(shù)和采樣位置)的G函數(shù)與其理論分布(本文在潛在采樣點(diǎn)中隨機(jī)抽取99組樣本近似計(jì)算G函數(shù)的95%置信限，樣本容量與優(yōu)化后采樣點(diǎn)數(shù)相同)比較，優(yōu)化采樣集的G函數(shù)值相對(duì)理論值越大，表明樣點(diǎn)越趨于聚集性空間分布，反之則樣點(diǎn)趨于規(guī)則(均勻)空間分布。F函數(shù)描述空間任意點(diǎn)(潛在采樣點(diǎn))到其最近采樣點(diǎn)的距離分布。仍然將優(yōu)化采樣集的F函數(shù)與相應(yīng)的理論分布(獲取方法與G函數(shù)相似)比較，優(yōu)化采樣集F函數(shù)值相對(duì)理論值越小，表明研究區(qū)域內(nèi)的間隔空間越少，樣點(diǎn)越趨于聚集性空間分布，反之則表明樣點(diǎn)趨于規(guī)則(均勻)空間分布。G函數(shù)和F函數(shù)計(jì)算式為

(4)

其中

I(ri≤r)=1

式中F′——G或F函數(shù)，%

r——預(yù)先設(shè)定的距離

I(·)——指示函數(shù)

n——采樣數(shù)(G函數(shù))或潛在采樣點(diǎn)數(shù)(F函數(shù))

ri——第i個(gè)采樣點(diǎn)到其他采樣點(diǎn)的最小距離(G函數(shù))或第i個(gè)潛在采樣點(diǎn)到采樣點(diǎn)的最小距離(F函數(shù))

科爾莫戈羅夫-斯米爾諾夫檢驗(yàn)(Kolmogorov-Smirnovtest，K-S)常用于檢驗(yàn)兩組數(shù)據(jù)是否來自同一個(gè)分布。本文用于檢驗(yàn)協(xié)變量樣本(由優(yōu)化后采樣點(diǎn)上的協(xié)變量取值組成)對(duì)協(xié)變量總體分布(潛在采樣點(diǎn)上協(xié)變量取值估計(jì))的代表性，即采用K-S法檢驗(yàn)優(yōu)化采樣集對(duì)屬性空間的覆蓋程度或代表性。

本文變異函數(shù)模型通過GS+軟件確定，其他計(jì)算均用Java語言編程完成；用Excel繪制變異函數(shù)圖，其他圖形在ArcGIS 10.0平臺(tái)上完成。

2 結(jié)果與分析

2.1 協(xié)變量選擇

2.1.1環(huán)境變量與土壤含鹽量的相關(guān)關(guān)系

協(xié)變量通常選擇與主變量(本研究為土壤含鹽量)關(guān)系密切且較容易獲取的環(huán)境變量。本文初步確定的備選變量為土質(zhì)(本研究區(qū)以粉砂質(zhì)壤土為主，土質(zhì)用砂粒含量表示；無采樣數(shù)據(jù)時(shí)，可用當(dāng)?shù)氐臄?shù)字化土壤地圖)、溝道分布(用與溝道的最短距離表征)和荒地分布(用與荒地的最短距離表征)。因研究區(qū)面積不大，且地形平坦，以往研究也表明研究區(qū)高程對(duì)土壤鹽分分布影響不明顯[19]，因此，備選變量未包括高程。不同采樣尺度上，3個(gè)備選變量與土壤含鹽量的皮爾遜相關(guān)系數(shù)見表1。由表1可以看出，大尺度(采樣分辨率為400 m×400 m)上土壤含鹽量與荒地分布、土質(zhì)均呈極顯著相關(guān)，且荒地分布與鹽分分布關(guān)系更密切，相關(guān)系數(shù)為-0.509；小尺度(采樣分辨率為50 m×50 m)上鹽分分布受荒地分布影響較小(相關(guān)系數(shù)僅為-0.057)，受土質(zhì)影響最大，其次為與溝道最短距離，且相關(guān)系數(shù)均達(dá)到極顯著水平。可見，不同尺度上影響土壤鹽分分布的主要環(huán)境變量不同，這表明協(xié)變量選擇受采樣尺度的影響。

表1 各環(huán)境變量與土壤含鹽量的皮爾遜相關(guān)系數(shù)Tab.1 Pearson correlation coefficient between environment variable and soil salt

2.1.2研究區(qū)協(xié)變量確定以及RK法中實(shí)測點(diǎn)處殘差計(jì)算

本文主要研究大尺度上采樣點(diǎn)的優(yōu)化布局，小尺度采樣僅是為了提高變異函數(shù)的計(jì)算精度，因此，選取在大尺度上與土壤鹽分相關(guān)系數(shù)最大的與荒地最短距離為協(xié)變量，不僅因?yàn)樗诖蟪叨壬吓c鹽分關(guān)系最密切，而且因?yàn)樵撟兞揩@取非常容易。此時(shí)，因?yàn)樾〕叨壬消}分與荒地分布關(guān)系不密切，故而，在RK法中，僅采用大尺度網(wǎng)格結(jié)點(diǎn)上的土壤含鹽量和與荒地最短距離確定回歸方程(y=-0.403 6x，x和y分別表示與荒地最短距離和土壤含鹽量)，分離趨勢項(xiàng)，計(jì)算所有點(diǎn)(包括大尺度和小尺度采樣點(diǎn))的殘差。

圖2為殘差的變異函數(shù)。圖2表明，所有點(diǎn)(189個(gè)點(diǎn))殘差的變異函數(shù)可較好地用指數(shù)模型擬合，小尺度上點(diǎn)的參與能更細(xì)致地刻畫較小滯后距上的變異特征，且圖形比較連續(xù)、無明顯突變；與僅用大尺度數(shù)據(jù)計(jì)算的殘差變異函數(shù)相比，小尺度數(shù)據(jù)加入后，塊金值縮小42.9%，塊金值/基臺(tái)值縮小32.2%(表2)?？梢?，用大尺度上的回歸方程分離趨勢項(xiàng)、計(jì)算殘差是合理的。

圖2 土壤鹽分殘差變異函數(shù)變化曲線Fig.2 Variogram of soil salt residuals

表2 變異函數(shù)模型參數(shù)Tab.2 Variogram model parameters

需說明的是：因?yàn)橥寥利}分符合對(duì)數(shù)正態(tài)分布(文中未列出有關(guān)計(jì)算)，因此，所有計(jì)算均是基于實(shí)測土壤含鹽量對(duì)數(shù)變換后的值進(jìn)行的。

2.2 優(yōu)化結(jié)果與評(píng)價(jià)

2.2.1合理采樣數(shù)確定

由式(3)可知，精度越高，對(duì)應(yīng)的采樣數(shù)越多。計(jì)算結(jié)果表明在95%置信度下相對(duì)誤差分別為5%、10%和15%時(shí)，合理采樣數(shù)分別為77、55和37。

初步選擇合理采樣數(shù)為55，但該采樣數(shù)僅是從統(tǒng)計(jì)學(xué)角度出發(fā)，為了保證該采樣數(shù)具有更好的合理性，繪制了不同數(shù)量樣點(diǎn)(從89個(gè)已知點(diǎn)中隨機(jī)抽取)用于回歸克里格插值時(shí)所有潛在采樣點(diǎn)平均估計(jì)方差的變化曲線(圖3)。由圖3可以看出，隨著樣點(diǎn)數(shù)的增加，平均估計(jì)方差逐漸減小并趨于平緩，樣點(diǎn)數(shù)為55時(shí)的平均估計(jì)方差相對(duì)于大尺度89個(gè)采樣點(diǎn)時(shí)的平均估計(jì)方差增大6.7%。因此認(rèn)為圖3選取55為優(yōu)化方案采樣數(shù)合理。

圖3 不同采樣數(shù)的平均回歸克里格估計(jì)方差變化曲線Fig.3 Mean RKV based on different sampling numbers

2.2.2SSA優(yōu)化結(jié)果

最終優(yōu)化后采樣點(diǎn)布局如圖4所示，分別以采樣方案的平均RKV最小和平均OKV最小為目標(biāo)函數(shù)，對(duì)采樣方案進(jìn)行優(yōu)化。潛在采樣點(diǎn)的平均RKV和平均OKV分別為0.080 7和0.114 3，RK法的估計(jì)方差估值精度相較OK法提高了29.39%。這表明以荒地分布作為協(xié)變量，采用RK法優(yōu)化采樣布局，可以有效提高土壤鹽分估值精度。

圖4 SSA優(yōu)化采樣點(diǎn)布局圖Fig.4 Optimized sampling point layout by SSA

2.2.3采樣布局評(píng)價(jià)

從估計(jì)方差的全局最優(yōu)性、采樣點(diǎn)對(duì)地理空間和屬性空間的代表性3方面，對(duì)基于RK法的優(yōu)化結(jié)果進(jìn)行評(píng)估。

(1)評(píng)估優(yōu)化結(jié)果的全局最優(yōu)性。本文跟蹤了優(yōu)化過程中目標(biāo)函數(shù)值的變化。由圖5可以看出，隨著溫度衰減次數(shù)增加，平均估計(jì)方差的總體趨勢下降，在下降過程中一直上下波動(dòng)，最終趨于平穩(wěn)。估計(jì)方差的上下波動(dòng)體現(xiàn)了SSA接受劣化解的能力，避免了優(yōu)化結(jié)果陷入局部最優(yōu)的狀況，隨著溫度衰減，SSA接受劣化解的概率逐漸下降，平均估計(jì)方差僅在小范圍內(nèi)波動(dòng)，最后不變。這表明優(yōu)化結(jié)果具有全局最優(yōu)性。

圖5 優(yōu)化過程中平均回歸克里格估計(jì)方差變化曲線Fig.5 Change of mean RKV during optimization

(2)評(píng)價(jià)優(yōu)化方案對(duì)地理空間的代表性。圖6為優(yōu)化后采樣點(diǎn)布局的G函數(shù)與F函數(shù)。由圖6可知，優(yōu)化后采樣點(diǎn)布局的G函數(shù)與F函數(shù)均位于隨機(jī)分布95%置信區(qū)間內(nèi)，即G函數(shù)與F函數(shù)均與其理論值無顯著差異，這表明優(yōu)化后的采樣點(diǎn)布局在地理空間內(nèi)既不是規(guī)則分布，也不是聚集型分布，而是隨機(jī)分布并較好覆蓋整個(gè)研究區(qū)。實(shí)際上，由圖4a也可以看出，雖然個(gè)別點(diǎn)相距較近，但并未出現(xiàn)大量點(diǎn)聚集在同一區(qū)域的現(xiàn)象。此外，由F函數(shù)可以看出任意潛在采樣點(diǎn)在600 m(鹽分變程)鄰域內(nèi)找到距其最近的采樣點(diǎn)的概率為96.5%，這為空間插值提供鄰域保證。因此，可以認(rèn)為優(yōu)化后的采樣點(diǎn)布局較好地覆蓋了地理空間，具有較好的地理空間代表性。

圖6 G函數(shù)和F函數(shù)變化曲線Fig.6 Changing curves of G-function and F-function

(3)評(píng)價(jià)優(yōu)化方案對(duì)屬性空間的代表性。如前所述，本文采用K-S法檢驗(yàn)優(yōu)化采樣集對(duì)協(xié)變量所表達(dá)的屬性空間代表性。由優(yōu)化后采樣點(diǎn)與荒地最短距離組成的樣本和由50 m×50 m離散結(jié)點(diǎn)與荒地最短距離組成的總體屬于同一分布。通過SPSS軟件分析，得到相伴概率P為0.275，遠(yuǎn)大于0.05，零假設(shè)成立，這表明優(yōu)化后的協(xié)變量樣本來自于協(xié)變量總體，優(yōu)化采樣集在荒地分布所表達(dá)的屬性空間上具有較好的代表性。

3 討論

本研究仍是基于地統(tǒng)計(jì)學(xué)理論對(duì)土壤采樣布局進(jìn)行優(yōu)化，這就要求研究變量的變異函數(shù)已知或具有估計(jì)的條件。當(dāng)研究變量的變異函數(shù)確定后，則可以計(jì)算出研究區(qū)平均估計(jì)方差，因此本文的研究方法可以適用于大多數(shù)土壤屬性的采樣點(diǎn)布局優(yōu)化，例如金屬離子、土壤有機(jī)質(zhì)等。對(duì)土壤鹽分而言，河套灌區(qū)的鹽分研究已經(jīng)積累了大量不同尺度、不同時(shí)期的鹽分采樣數(shù)據(jù)[20-24]，這可為鹽分變異函數(shù)模型的建立和采樣方案的優(yōu)化提供有力支撐；研究還表明，河套灌區(qū)某一區(qū)域周年內(nèi)鹽分空間分布格局相對(duì)穩(wěn)定[23]，這也有利于鹽分采樣布局的優(yōu)化，綜上本文的研究方法可以很好地適用于河套灌區(qū)的鹽分采樣點(diǎn)布局優(yōu)化。而當(dāng)無過往數(shù)據(jù)或先驗(yàn)方差可利用時(shí)，本文的方法則難以使用。為獲取精度較高的變異函數(shù)模型，可使用嵌套采樣，此時(shí)，如何在綜合考慮不同尺度環(huán)境協(xié)變量基礎(chǔ)上，設(shè)計(jì)嵌套采樣方案還需要進(jìn)一步探討。

本文研究方法適用于所有存在鹽漬化問題的區(qū)域，但由于自然因素和人類活動(dòng)的不同，不同地區(qū)影響土壤鹽分分布的因素也不同。如張芳等[25]研究表明，新疆維吾爾自治區(qū)奇臺(tái)綠洲地形因素對(duì)表層和中層土壤鹽分分布格局影響較大，其中高程與鹽分呈極顯著相關(guān)關(guān)系；劉洋等[26]研究認(rèn)為黃河三角洲地區(qū)土壤鹽分空間變異程度與距海洋遠(yuǎn)近以及人類活動(dòng)等因素有關(guān)；李亮等[27]研究結(jié)果表明，河套灌區(qū)內(nèi)荒地在作物生育期有排水積鹽作用，荒地對(duì)土壤含鹽量存在影響。本文選擇荒地分布作為協(xié)變量，主要因?yàn)樗潜姸嗯c土壤鹽分關(guān)系密切的環(huán)境變量中較容易獲取的。此外，本研究中大、小尺度上各環(huán)境變量和土壤鹽分的相關(guān)分析以及其他相關(guān)研究[28]表明，同一個(gè)環(huán)境變量在不同研究尺度上與主變量關(guān)系的密切程度可能不同，不同尺度上土壤鹽分分布的主控因子不同。因此，在實(shí)際使用本文研究方法中，應(yīng)根據(jù)地區(qū)特點(diǎn)和研究尺度選擇合適的協(xié)變量。

綜上所述，本研究以研究區(qū)內(nèi)平均估計(jì)方差最小為目標(biāo)，運(yùn)用SSA對(duì)土壤鹽分采樣點(diǎn)布局進(jìn)行優(yōu)化，優(yōu)化后的采樣點(diǎn)布局在地理空間和屬性空間均能較好地覆蓋。研究結(jié)果對(duì)提高土壤鹽分分布圖繪制精度或降低野外采樣和室內(nèi)分析成本具有重要意義。

4 結(jié)論

(1)在采樣數(shù)為55時(shí)，以荒地分布作為協(xié)變量、運(yùn)用RK法優(yōu)化土壤鹽分采樣點(diǎn)布局，可使50 m×50 m潛在采樣點(diǎn)上的平均RKV(0.080 7)相較OK法優(yōu)化結(jié)果的平均OKV(0.114 3)縮小29.39%；F函數(shù)圖和G函數(shù)圖表明，優(yōu)化后的采樣點(diǎn)布局在地理空間具有較好的代表性；K-S檢驗(yàn)則表明該采樣點(diǎn)布局在荒地分布所代表的屬性空間亦具有良好的代表性。

(2)因不同研究尺度上，影響土壤鹽分分布的主要環(huán)境因素不同，所以協(xié)變量具有尺度依賴性。本研究區(qū)在大尺度(采樣網(wǎng)格400 m×400 m，范圍1 574 hm2)上以荒地分布為協(xié)變量是合理的，但小尺度上或多尺度嵌套采樣優(yōu)化時(shí)，協(xié)變量需要另外擇優(yōu)選用。