空間多觀測樣本的地理加權回歸模型?

栗春曉, 李芙蓉

(中國海洋大學數學科學學院, 山東 青島 266100)

隨著衛星遙感、地理信息系統和計算機技術不斷發展,空間數據被廣泛應用于生態學、林業、醫學和海洋科學等各領域。空間數據具有空間異質性,這意味著隨著地理位置的變化,變量間的關系或結構會發生變化,這種空間異質性普遍存在于空間數據中。例如,日常生活中房價及其影響因素的關系在不同區域具有明顯的差異;海洋科學研究發現溫度和鹽度的關系隨著地理位置和深度而變化[1]。在這種情況下,普通最小二乘回歸(Ordinary least squares regression,OLSR)不能準確、全面地描述空間數據的真實特征,因為該模型是全局模型,沒有考慮地理位置的影響,得到的結果只是所研究區域的一種平均值[2-3]。

為了進一步研究空間數據的特性,Brunsdon等[2]基于變系數回歸模型提出了地理加權回歸(Geographically weighted regression,GWR)模型。該模型在估計某個空間點上的回歸系數時,可以利用其鄰近空間點上的觀測樣本信息(每一空間點僅有一個觀測樣本),并按照距離地遠近施以不同的權重來進行估計,這克服了傳統回歸模型的回歸系數為常數的缺點[3]。在過去20年中,GWR成為探索空間非平穩關系的重要工具,被廣泛用于生態學、林業、房地產等領域[4-9]。然而,GWR模型也存在很多局限性,例如,缺乏統一的方法來估計局部系數,離群值的強烈影響,弱數據問題等[10-11]。

為了克服GWR模型的局限性,學者們也提出了多種改進方法。Lesage[10,12]借助于貝葉斯方法提出了貝葉斯地理加權回歸模型(Bayesian geographically weighted regression,BGWR),Charlton[13]提出了同時包含全局變量和局部變量的混合地理加權回歸模型(Mixed geographically weighted regression,MGWR),Huang等[14]提出了時空地理加權回歸模型(Geographi-cally and temporally weighted regression,GTWR),假設回歸系數是變系數回歸模型中地理位置和觀察時間的函數。研究發現,這些關于GWR模型以及GWR改進模型的研究都是基于每個空間點上僅有一個觀測樣本的情況。然而,在實際應用中,不同的空間位置上往往有多個觀測樣本。例如,在研究海表面熱通量與中尺度渦旋所攜帶的海表面溫度異常關系的季節變化時,同一季節不同年份的數據就可以視作多個觀測樣本;在研究不同普查區域家庭收入與學齡兒童體重指數(BMI)的關系時,同一普查區域不同兒童的數據同樣可以視作多個觀測樣本。

當每個空間點上的觀測樣本數量相同時,一種可能的策略是分別對每組觀測樣本(這里定義所有空間點上的一次觀測為一組觀測)利用GWR模型對回歸系數進行估計,再對基于每組觀測樣本得到的回歸系數求平均數或中位數,從而得到最終回歸系數的估計。然而,該想法忽略了每組觀測樣本所得回歸系數都應相同這一事實。此外,當每個空間點上的樣本數量不一樣時,甚至無法定義一組觀測樣本,使得已有的GWR模型根本無法應用。另一種可能的策略是對每個點上的所有觀測樣本分別進行普通最小二乘回歸,但這未能利用回歸關系在鄰近的空間點上具有相近值這一特點,導致估計效果不理想。特別的,當每個空間上的樣本數量較少時,會導致嚴重的估計誤差。

鑒于此,本文對經典GWR模型進行拓展,構建了一種可以處理空間點上任意樣本數量的地理加權回歸模型(Multiple-replicate geographically weighted regression,MRGWR)。MRGWR模型在估計回歸系數時,充分利用了回歸關系在鄰近點上具有相似性的特點,并對鄰近空間點的觀測樣本施加不同權重。本文進一步基于數值仿真實驗與實證研究,對模型的有效性與適用性進行了分析。

1 模型

1.1 GWR模型

假設(ui,vi)是第i個空間點位置,{(yi,xi),i=1, 2, …,n}是觀測樣本在觀測位置(u1,v1)、 (u2,v2), …, (un,vn)處的觀測值,xi=[xi1,xi2, …,xip]T表示觀測位置(ui,vi)處的p維解釋變量,yi為對應的響應變量。將經典OLSR模型推廣到空間上,則GWR模型可寫成

(1)

式中:βik表示第k個解釋變量在(ui,vi)處所對應的回歸系數;εi為獨立同分布的白噪聲過程均值為0,方差為σ2。對于GWR模型,βi=[βi1,βi2, …,βip]T的估計值為

(2)

其中,

Wi=W(ui,vi)為n×n空間權重矩陣,其對角線元素wij表示觀測點i與觀測點j之間的核函數,wij與點i到點j的空間距離相關。

1.2 MRGWR模型

模型(1)假定在每一個空間點(ui,vi)上,僅有一個觀測樣本(yi,xi),然而,在許多空間數據集中,不同的空間位置上往往有多個觀測樣本。考慮下面這種情況,空間點(ui,vi)上具有hi個觀測樣本,樣本觀測值為

(yi1;xi11,xi12, …,xi1p), (yi2;xi21,xi22, …,xi2p), …, (yihi;xihi1,xihi2, …,xihip)。

其中,(yil;xil1,xil2, …,xilp),l=1, 2, …,hi是觀測點在(ui,vi)處響應變量和解釋變量的第l個觀測值,i=1, 2, …,n。

此時,模型(1)需改寫為

(3)

MRGWR模型作為GWR模型的推廣,沿用了GWR模型中最基本的思想,即估計(ui,vi)處的回歸系數時,可以利用其鄰近空間點上的觀測樣本信息,并按照距離地遠近施以不同的權重來進行估計。在MRGWR模型中,回歸系數的估計值應使得以下加權殘差平方和取得最小值

(4)

其中,wij為估計(ui,vi)處的回歸系數時,(uj,vj)處的觀測樣本所被賦予的權重。wij與(ui,vi)和(uj,vj)之間的距離有關,其具體表達式將在1.3節詳細說明。需要指出的是,wij與樣本本身無關,同一點上的不同樣本被賦予同樣權重。

公式(4)可以寫成以下矩陣形式:

其中,

Wij=diag[wij,…,wij]hj×hj,

X=[X1,X2,…,Xj,…,Xn]T,

y=[y11, …,y1h1, …,yn1, …,ynhn]T,βi=[βi1,βi2, …,βip]T,Wij是(ui,vi)處觀測樣本與(uj,vj)上hj個觀測樣本之間地理距離的空間權重矩陣。

由此解出空間點(uj,vj)上的回歸系數為

(5)

圖1 GWR模型(a)和MRGWR模型(b)的估計原理示意圖Fig.1 Schematic of GWR model (a) and MRGWR model (b) estimation principles

1.3 核函數和帶寬

符合以上要求的空間核函數類型較豐富。常用的函數類型包括Gauss核函數、Exponential核函數、Bi-square核函數、Tri-cube核函數等[15]。在本文中,采用的是Exponential核函數,其表達式為

wij=exp(-‖dij‖/b)。

(6)

其中,b為帶寬參數,控制著權重隨著距離增加的衰減速度,其大小決定了擬合曲面的光滑性。

與核函數的選擇相比,回歸系數對于帶寬的選擇表現得更敏感,帶寬過大回歸系數估計的偏差過大,帶寬過小又會導致回歸系數估計的方差過大。因此,選取合適的帶寬參數對于MRGWR模型估計的準確性至關重要。目前,常用的最優帶寬選取準則有交叉驗證方法(CV)、廣義交叉驗證方法(GCV)以及修正的赤池信息準則(AICC)[16-19]。本文采用CV方法,如果有需要的話,也可以直接采用其他準則。

2 數值仿真實驗

本節通過數值仿真實驗來評估普通最小二乘回歸,GWR和MRGWR模型在不同情況下的估計性能。

2.1 實驗設計

模擬實驗中,選取平方域[0, 1]×[0, 1]構成地理區域,從中隨機生成500個空間位置,每個位置上的響應由下式產生:

yil=βi1xil1+βi2xil2+βi3+εil。

(7)

在以往的數值實驗模擬中,解釋變量的值通常從白噪聲過程中產生[20-21]。在海洋科學的許多應用中,用于擬合回歸模型的解釋變量具有明顯的空間結構[22]。因此,為了模擬實際情況,在數值實驗中基于空間過程生成解釋變量。具體來說,先定義{z1(ui,vi)}和{z2(ui,vi)},i=1, 2, …, 500是服從空間高斯過程的兩個獨立觀測,且其均值為零,變量協方差矩陣定義為

在數值實驗中,回歸系數被設計成具有空間上連續變化的形式,并由一個空間高斯過程生成。具體而言,βi=[βi1,βi2,βi3]服從空間高斯分布過程,且均值為0,協方差矩陣由一個各向異性的指數函數定義

Cov(βik,βjk)=

基于以上模型,設計以下三組實驗來評估所建立的MRGWR、GWR和OLSR模型系數估計的性能。

實驗1:解釋變量具有不同的空間自相關性。在研究中設定σ=1,h=10,φ=[0.1, 0.3, 1]。

實驗2:空間點上具有不同的觀測樣本數量。在研究中設定σ=1,φ=0.3,h=[3, 10, 30]。

實驗3:回歸模型具有不同的信噪比。在研究中設定h=10,φ=0.3,σ=[0.3, 1, 3]。

每組實驗均開展100次數值模擬來評估不同模型的估計效果。為了量化每種模型估計的性能,引入估計的均方誤差(MSE),并定義為

(8)

2.2 實驗結果分析

2.2.1 解釋變量具有不同的空間自相關性 在實驗1中,通過設定φ=[0.1, 0.3, 1]來改變解釋變量的空間自相關性,進而評估在不同φ值下回歸模型的估計性能。圖2展示了回歸系數的真值以及當φ=0.3時,OLSR、GWR和MRGWR三種模型的估計結果。其中,GWR模型的回歸系數估計值是基于每組觀測樣本得到的回歸系數求平均得到的(定義所有空間點上的一次觀測為一組觀測)。

((a)—(c)回歸系數β1、β2和β3的空間結構;(d)—(f)OLSR模型的估計結果;(g)—(i)GWR模型的估計結果;(j)—(l)MRGWR模型的估計結果。該模擬中,空間自相關性φ=0.3, 觀測樣本數量h=10, 標準差σ=1,右側圖示表示回歸系數的大小。(a)—(c) denote the spatial structure of the regression coefficients β1、β2 and β3; (d)—(f) denote the estimation results of OLSR model; (g)—(i) denote the estimation results of GWR model; (j)—(l) denote the estimation results of MRGWR model. In this simulation, spatial autocorrelation φ=0.3, number of replicates h=10, standard deviation σ=1, right-hand icons indicate the magnitude of the regression coefficients.)圖2 回歸系數估計值對比Fig.2 Comparison of regression coefficient estimates

從圖2可以看出,MRGWR估計的回歸系數與真實值更加一致, GWR和OLSR模型估計的回歸系數與真實值差異較大。一方面,OLSR模型未能利用回歸關系在鄰近的空間點上具有相近值這一特點,導致估計效率果不理想。另一方面,盡管GWR模型利用了空間相鄰點的觀測信息,但它忽略了對于每組觀測樣本回歸系數都是相同的這一事實。相反,MRGWR模型估計結果更準確,因為它充分利用了回歸關系在鄰近點上具有相似性的特點,并對鄰近空間點的觀測樣本施加不同的權重。

表1展示了MRGWR、GWR以及OLSR模型在解釋變量具有不同強度的空間自相關性之下回歸系數的估計誤差。由表1可知,在解釋變量具有三種強度的空間自相關性下,GWR和MRGWR模型的估計性能均顯著優于OLSR模型,且MRGWR模型的估計性能相比于GWR模型更優。隨著空間自相關性增強(即φ不斷增大),OLSR和GWR的估計誤差逐漸變大,雖然MRGWR的估計誤差也在變大,但幅度很小,其值遠低于OLSR和GWR的估計誤差。這說明當解釋變量存在強空間自相關性時,MRGWR相對于GWR估計效率上的優勢更加突出。例如當φ=0.1時,MRGWR模型的MSE約為GWR模型的1/2;但當φ=1時,MRGWR模型的MSE僅為GWR模型的1/4。

表1 不同解釋變量間自相關性下各模型100次模擬的平均均方誤差Table 1 Average mean squared error for 100 simulations of each model under different spatial autocorrelation in the explanatory variables

2.2.2 不同的觀測樣本數量 在實驗2中,我們通過設定h=[3, 10, 30]來評估不同模型在空間點上不同觀測樣本數量下的性能,表2比較了三種模型在100組模擬下的平均MSE。

由表2可知,對于空間點上各種觀測樣本數量,MRGWR均保持較低的誤差,且都優于OLSR和GWR模型。此外,隨著h值的減少,三種模型的MSE都變大,估計性能隨之降低,這種退化對OLSR來說是最嚴重的,但對GWR和MRGWR來說MSE的變動幅度要小很多,因為它們利用的是鄰近的觀測樣本。在h=3或10時,雖然GWR要顯著優于OLSR,但當h=30時,OLSR卻更優于GWR,而MRGWR模型的MSE都小于OLSR模型。由此可見,相比于OLSR和GWR模型,MRGWR模型表現更優。

表2 不同的觀測樣本數量下各模型100次模擬的平均均方誤差Table 2 Average mean squared error for 100 simulations of each model under different number of replicates

2.2.3 不同的信噪比 在實驗3中,我們通過設定標準差σ=[0.3, 1, 3]來改變回歸模型的信噪比,進而評估三種不同信噪比設置下回歸模型的性能。表3比較了三種模型在100組實驗下的平均MSE。從表3中可以看出,對于模型的各種信噪比,MRGWR模型的估計性能總是顯著優于OLSR和GWR模型。此外,隨著σ值從0.3變為1或3,模型信噪比從較高水平變為中等水平或較低水平,OLSR的MSE大幅增加,且波動幅度較大,而MRGWR和GWR模型的MSE增幅較小。特別的,當σ=0.3時,GWR模型的MSE大于OLSR,但當σ=1或3時,結果相反。相比之下,在所有三種情況下,MRGWR估計的MSE總是小于OLSR估計的MSE,且MRGWR的分布較為集中。

表3 不同的標準差下各模型100次模擬的平均均方誤差Table 3 Average mean squared error for 100 simulations of each model under different standard deviation

2.3 數值仿真實驗總結

為了檢驗模型在不同情況下的估計性能,本文分別從解釋變量的空間自相關性、空間點上觀測樣本數量以及模型信噪比三個方面展開具有空間結構的數值仿真實驗。結果表明,無論解釋變量的空間自相關性、觀測樣本數量以及模型信噪比如何變化,MRGWR模型的估計性能總是優于普通最小二乘回歸和GWR模型。特別是當解釋變量在空間上高度相關、觀測樣本數量較少、模型信噪比較低時,普通最小二乘回歸和GWR模型的估計性能會存在不同程度上惡化。相比之下,MRGWR模型在這種情況下仍能產生合理的估計,表明其具有較強的穩健性。

3 實際應用

3.1 問題背景

中尺度渦旋占據了海洋中90%的動能,是海洋上層最顯著的運動形式,對海洋中熱量、溶解氣體以及營養物質的輸送起著重要的作用[23]。傳統的中尺度渦旋研究主要是從海洋自身的動力過程出發,缺少對渦旋與大氣間相互作用的考慮[24-26]。然而,近十幾年的觀測和數值模擬結果表明海洋中尺度渦旋與大氣之間存在著強烈的相互作用,反過來可以給海洋中尺度渦旋帶來反饋[27-30]。因此,中尺度渦旋與大氣相互作用已經成為當前物理海洋和氣候研究的一個前沿領域。

中尺度渦旋與大氣之間存在多種形式的相互作用。特別是中尺度渦旋攜帶的海表面溫度(Sea surface temperature,SST)異常可以對海氣界面的熱通量產生強烈的影響[1,29,31-32]。也就是說,暖(冷)異常會促進(抑制)熱量從海洋進入大氣。除了極大地影響大氣邊界層的動態過程外[33-34],這種中尺度渦旋引起的海面熱通量異常反過來又抑制了中尺度SST異常,成為渦旋勢能耗散的重要途徑[29],這一過程被稱為渦旋熱反饋[35]。渦旋熱反饋的強度可以通過中尺度海面熱通量異常對中尺度SST異常的回歸系數來衡量。然而,由于這種回歸關系受到背景大氣的動力和熱力學結構的調控,存在明顯的季節性和空間差異,因此難以通過經典的OLSR模型進行估計。

鑒于此,本文通過MRGWR模型來估計中尺度SST異常與海面凈熱通量異常間(本文中的熱通量定義向上為正)的回歸關系,進而揭示中尺度海面凈熱通量異常對中尺度SST異常的響應關系隨季節和空間的變化特征。

3.2 數據來源

本節數據來自第三代日本海洋通量數據集(J-OFURO3),包含1988—2017年(30年)的SST和海面凈熱通量的月平均場,水平空間覆蓋的范圍為全球海域(90°S—90°N,180°E—180°W),可通過以下網址下載:https://j-ofuro.isee.nagoya-u.ac.jp/en/dataset/entry-114.html。該數據的空間分辨率為0.25°×0.25°,能夠分辨出大部分中尺度渦旋引起的SST異常。為了提取中尺度SST和海面凈熱通量異常,我們首先對原始數據進行4°×4°的空間滑動平均,并剔除該滑動平均的信號,將剩余的異常場定義為中尺度異常。

為了能夠較好表征中尺度渦旋熱反饋,本文選取北太平洋(10°N—65°N,105°E—90°W)為研究區域,見圖3紅色矩形區域。

(紅色矩形為研究區域,白色為海域。Red rectangle is the study area, white is the sea area.)圖3 實驗海域Fig.3 Experimental sea area

3.3 模型應用

由于本文關心的是中尺度海面凈熱通量異常對中尺度SST異常的響應關系隨空間和季節的變化,所以對于不同月份(1月、2月……12月)的數據分別利用MRGWR模型進行估計,并將不同年份同一月份的數據視作多次觀測,即每個月每個空間點上都有30次觀測樣本。在此假定下,對于每個月,可建立如下回歸模型

Qil=βi1Til+βi0+εil。

式中:Qil表示在位置點(ui,vi)上第l年的中尺度海面凈熱通量異常;Til表示在位置點(ui,vi)上第l年的中尺度SST異常;βi1是這兩個變量之間的空間回歸系數,用于描述中尺度海面凈熱通量異常對中尺度SST異常的響應關系;βi0是空間回歸模型的截距;εil表示隨機誤差項。為了降低計算負擔,在北太平洋海域非零的觀測值地理位置中,隨機選取10 000個點進行估計。

圖4展示了由MRGWR模型估計得到的中尺度海面凈熱通量異常與中尺度SST異常之間的空間回歸系數β1的季節和空間的變化特征。從空間上來看,中尺度海面凈熱通量異常對中尺度SST異常的響應存在顯著的南北向差異,在副熱帶海區β1可達60～70 W/(m2·K),而在副極地海區β1則只有40 W/(m2·K)左右。此外,β1在黑潮延伸體區域呈現出局地增強特征。這些空間特征背后的物理機制值得相關領域的學者深入分析。在季節變化方面,β1整體呈現冬強夏弱的特征,這種季節變化在中高緯度要比低緯度更明顯。由于海氣間的熱量交換顯著依賴于海表面風速,因此中高緯度冬季的β1大概率是由冬季風暴過程所引起的較強風場導致的,但是否還有其他控制因素還需要進一步分析。綜上所述,中尺度海面凈熱通量異常對中尺度SST異常的響應關系存在顯著的季節和空間變化特征,MRGWR模型提供了一種可靠的提取該特征的統計方法。

(右側圖示表示β1的大小,單位為W/(m2·K)。陸地被灰色掩蓋。Right-hand icons indicate the magnitude of β1 in W/(m2·K).Land is masked by grey.)圖4 不同月份下β1的空間分布Fig.4 Spatial distribution of β1 for different months

4 結語

本文針對現有的GWR模型以及GWR改進模型無法處理空間點上多個樣本的情形,對GWR模型進行了拓展,構建了一種可以處理空間點上任意樣本數量的地理加權回歸模型——MRGWR。MRGWR模型在估計回歸系數時,充分利用了回歸關系在鄰近點上具有相似性的特點,并對鄰近空間點的觀測樣本施加不同權重。通過開展數值仿真實驗,評估了在不同的解釋變量空間自相關性、不同的觀測樣本數量以及不同的模型信噪比下,MRGWR模型的估計性能。數值實驗結果表明,與普通最小二乘回歸和GWR模型相比,MRGWR模型能夠有效地處理空間點上有多個觀測樣本的情形,且估計性能總是優于普通最小二乘回歸和GWR模型。本文進一步將MRGWR模型用于探究北太平洋中尺度海面凈熱通量異常對中尺度SST異常的響應關系,研究發現,中尺度海面凈熱通量異常對中尺度海面溫度異常的響應關系存在顯著的正相關以及明顯的空間和季節變化。