付昊天, 李芙蓉

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

空間與時間是一切人類活動、社會事件和物理現(xiàn)象的兩個基本維度。隨著觀測、實驗和數(shù)值模擬技術的發(fā)展,時空數(shù)據(jù)的來源已經十分豐富,處理具有時空依賴性的觀測數(shù)據(jù)已成為地球科學、環(huán)境科學、生物學和醫(yī)學等研究領域的一個關注點。時空數(shù)據(jù)的一個重要特征是不同變量間的關系在空間與時間上可能存在著顯著的變化。經典的回歸模型由于假定回歸系數(shù)為常數(shù),所以無法捕捉響應變量與協(xié)變量之間的回歸關系在時空域上復雜的變化特征。

變系數(shù)回歸模型為刻畫時空變量的關系提供了可能。早期的變系數(shù)回歸模型主要是用于估計響應變量與協(xié)變量之間的回歸關系在空間上的變化特征。例如,Brunsdon等[1]提出的地理加權回歸(Geographically weighted regression, GWR)模型是一種用于解決空間非平穩(wěn)性的空間統(tǒng)計方法,允許響應變量和協(xié)變量之間的關系隨著地理位置的變化而變化。GWR模型基于局部回歸的思想,在估計某個空間點上的回歸系數(shù)時,利用其鄰近空間點上的觀測信息,并按照距離遠近施以不同的權重來進行估計。在過去20年中,GWR模型被廣泛用于地球科學、農業(yè)和社會經濟等領域。Huang B等和Fotheringham等[2-3]通過將權重函數(shù)從空間推廣至時空域,構建了時空地理加權回歸模型(Geographically and temporally weighted regression, GTWR),實現(xiàn)了對隨時空變化的回歸系數(shù)的估計。Gelfand等[4]提出的空間變系數(shù)(Spatially varying coefficient, SVC)模型將回歸系數(shù)看作是一個多元空間高斯過程的實現(xiàn),并利用貝葉斯的方法進行估計。Song等[5]對SVC模型進行了拓展,建立了時空變系數(shù)(Spatiotemporally varying coefficients, STVC)模型。

GTWR和STVC模型都假定回歸系數(shù)隨時空連續(xù)變化。然而,在很多實際應用中回歸系數(shù)會呈現(xiàn)出團狀結構,即回歸系數(shù)可能會在時空域的某些子區(qū)域內保持相對穩(wěn)定,而在不同子區(qū)域間的界面上產生突變。例如,社會經濟活動中,社區(qū)房產價格和教育水平間的關系可能會受到政府宏觀政策的影響在時間上發(fā)生突變;海洋學研究中,海水物理化學性質間的關系在同一水團內部保持相對均勻,但在不同水團的邊界上產生顯著變化。

針對回歸系數(shù)在空間上的團狀結構,Li和Sang[6]提出了空間團狀系數(shù)回歸(Spatially clustered coefficient, SCC)模型。SCC模型結合了統(tǒng)計學中的變量選擇理論和圖論中的最小生成樹理論,通過對最小生成樹所連接的空間點上的回歸系數(shù)的差異實施懲罰來捕捉回歸系數(shù)在空間上的團狀結構。本文將SCC模型拓展至時空域,提出一種可以估計在時空域上呈團狀分布的回歸系數(shù)模型,即時空團狀系數(shù)回歸(Spatio-temporally clustered coefficient, STCC)模型。STCC模型通過對時空域內相鄰點上的回歸系數(shù)進行懲罰,來促進時空鄰近回歸系數(shù)之間的同質性。

1 STCC模型介紹

假設一組時空觀測數(shù)據(jù){(x(si,tj),y(si,tj)),i=1,…,n,j=1,…,m},其中s1,…,sn∈Rd(d≥2)代表空間位置,t1,…,tj∈R代表時間位置,x(si,tj)=(x1(si,tj),…,xp(si,tj))T為p維解釋變量,y(si,tj)為響應變量。考慮時空變系數(shù)線性回歸

(1)

式中:βk(si,tj)是隨時空變化的回歸系數(shù);ε(si,tj)是均值為0、方差為σ2的獨立同分布隨機噪聲。當xp(si,tj)≡1時,βp(si,tj)代表隨時空變化的截距。

假定在同一時空位置僅有一次觀測,則可將式(1)寫作矩陣形式

y=Xβ+ε。

(2)

式中:β=(β1,…,βp)T;βk=(βk(s1,t1),…,βk(sn,t1),…,βk(s1,tm),…,βk(sn,tm));(n·m)×(n·m·p)維設計矩陣X=[diag(X1),…,diag(Xp)],Xk=(xk(s1,t1),…,xk(sn,t1),…,xk(s1,tm),…,xk(sn,tm))T。

對于p>1,n·m·p>n·m。因此,無法直接使用普通最小二乘擬合來估計式(2)中β的值,故需將其正則化。假定時空域中相鄰位置的β傾向于一致,則可通過最小化以下目標函數(shù)來估計β:

(3)

式中:P是一個懲罰矩陣,用于促使相鄰位置回歸系數(shù)之間的同質性;λ是懲罰參數(shù),決定懲罰的強度。接下來將討論P和λ具體的選擇方法。

1.1 選擇懲罰矩陣P

我們分別對空間上和時間上的相鄰點進行懲罰。在空間上,我們采用與SCC模型一致的空間懲罰矩陣,即通過最小生成樹將同一時刻的不同空間位置連接,被連接的兩個點被認為是鄰近的,對其回歸系數(shù)間的差異進行懲罰。在時間上,懲罰可采用不同的形式,如對時空位置(si,tj),j≠1,將其與(si′,tj-1)間的回歸系數(shù)差異進行懲罰,其中si′為si的空間位置相鄰點。此處我們采用最簡單的形式,即令i′=i。即

(4)

式中:H為(n-1)×n維的空間懲罰矩陣;O為(n-1)×n維的零矩陣;In為n×n維的單位矩陣。τ∈[0,∞)為時空懲罰比,用于控制時間懲罰與空間懲罰的相對大小。

當τ=0時,STCC模型相當于使用SCC模型對每個時刻的回歸系數(shù)進行估計;而當τ→∞時,STCC模型相當于使用Fused LASSO模型[7]對每個空間點上的回歸系數(shù)進行估計。對于一般的τ>0,懲罰矩陣P同時對空間和時間上的相鄰點的回歸系數(shù)差異進行懲罰。為同時促使回歸系數(shù)在空間和時間上的同質性,形成時空域上的團狀結構,我們取τ=1。

1.2 選擇懲罰參數(shù)λ

在STCC模型中,懲罰參數(shù)控制了回歸系數(shù)在時空域中的同質性程度。當λ→∞時,模型的回歸系數(shù)為常數(shù),即將整個時空域視為一個團;當λ→0時,STCC模型退化為逐點最小二乘擬合,即每個時空觀測點單獨成團。可采用廣義交叉驗證、AIC、BIC、EBIC等模型選擇準則來確定最優(yōu)的λ值。與Li和Sang[6]相同,我們這里采用BIC準則。

2 計算方法

為了便于表示,將針對p個解釋變量的懲罰矩陣P合寫為一個矩陣

(5)

式中:P=n·m·p,Q=(n-1)·m·p+n·(m-1)·p。此時式(3)可以改寫為

(6)

式(6)可看作是一種特殊的Fused LASSO形式,這意味著可以使用線性交替方向乘子法(Linear alternating direction method of multipliers, LADMM)算法來估計β。

(7)

式(7)的增廣拉格朗日函數(shù)是

(8)

第一,對迭代中β子問題進行求解。

(9)

(10)

(11)

第二,對迭代中γ子問題進行求解。

(12)

式中soft為軟閾值函數(shù),軟閾值函數(shù)可由下式定義:

soft(a,b)=sign(a)·max{0,|a|-b}。

(13)

式中sign為符號函數(shù)。

第三,對迭代中α子問題進行求解。

(14)

3 數(shù)值仿真研究

本節(jié)通過數(shù)值仿真實驗來探究STCC方法在不同情形下的穩(wěn)健性。為了便于展示,我們考慮二維空間(即d=2)中的三種情形。在前兩種情形中,回歸系數(shù)均在時空域中呈團狀結構,第三種情形模擬了隨時空連續(xù)變化的回歸系數(shù),用以探究STCC方法的靈活性。

在[0,1]×[0,1]的正方形區(qū)域中隨機生成n個空間位置,每個空間位置選取8個等間隔的時間節(jié)點。取n=250,n=500,n=1 000分別對應觀測樣本數(shù)的少、中、多三種不同情況。時空位置確定后,響應變量y與預測變量x1和x2間的關系由以下公式刻畫:

y(si,tj)=x1(si,tj)β1(si,tj)+
x2(si,tj)β2(si,tj)+β3(si,tj)+ε(si,tj)。

(15)

式中:β1(si,tj)和β2(si,tj)為隨時空變化的回歸系數(shù);ε(si,tj)～N(0,σ2)獨立同分布,令σ=0.1,0.2,0.4,分別對應信噪比強度的高、中、低的情況。

對于時空數(shù)據(jù),x1和x2通常具有明顯的時空結構。為了模擬這種特征,我們仿照Li和Sang[6],引入均值為0的輔助時空變量z1(si,tj)與z2(si,tj),其具有以下時空平穩(wěn)協(xié)方差矩陣:

(16)

在以下分析中,我們將對比GTWR、SCC和STCC三種不同方法。對于STCC模型,式(4)中的時空懲罰比τ設為1,即對時空懲罰是同等力度的。在LADMM算法中,設定式(9)中的μ=1,收斂精度ξ=5。對SCC與STCC模型,懲罰參數(shù)λ采用BIC收斂準則來選取。對GTWR方法,采用指數(shù)核函數(shù),其最佳帶寬由交叉驗證法選擇。為了量化每種估計方法的性能,考慮估計的均方誤差MSEβ,定義為

(17)

在每種情形下,我們分別比較了時空相關性、觀測樣本數(shù)、信噪比的變化對模型估計的影響,將100次獨立數(shù)值模擬實驗平均的MSEβ作為最終的MSEβ。比較時空相關性時,設定n=1 000,σ=0.1;比較觀測樣本數(shù)時,設定(φ1,φ2)=(0.3,4.2),σ=0.1;比較信噪比時,設定(φ1,φ2)=(0.3,4.2),n=500。

3.1 情形1

情形1模擬的是團狀結構的融合過程,即空間中的兩個相鄰團在某時間點后合二為一,新團的性質與之前的兩個團皆不相同(見圖1)。β1,β2,β3分別在時間節(jié)點t=2,4,6后發(fā)生融合。表1比較了情形1下三種方法的MSEβ。STCC和SCC的估計效果明顯優(yōu)于GTWR,特別是樣本觀測數(shù)多、信噪比高、時空相關性強時,STCC和SCC的MSEβ只有GTWR方法的1/10左右。此外,對任意程度的時空相關性、觀測樣本數(shù)、信噪比,STCC的MSEβ均小于SCC的MSEβ。綜上所述,STCC的估計效果最優(yōu)。

表1 情形1下SCC、STCC和GTWR方法的100次數(shù)值模擬實驗的平均MSEβ

(橫軸與縱軸為[0,1]×[0,1]的正方形的兩個維度,無量綱。The horizontal and vertical axes are two dimensions of the [0,1]×[0,1] square region, dimensionless.)

3.2 情形2

情形2模擬的是團狀結構的形變過程,即空間中的某個團在某時間點后形狀發(fā)生變化,但性質保持不變(見圖2)。β1,β2,β3分別在時間節(jié)點t=2,4,6后發(fā)生形變。表2比較了情形2下三種方法的MSEβ。與情形1類似,STCC的MSEβ最小,GTWR的MSEβ最大。

表2 情形2下SCC、STCC和GTWR方法的100次數(shù)值模擬實驗的平均MSEβ

(橫軸與縱軸為[0,1]×[0,1]的正方形的兩個維度,無量綱。The horizontal and vertical axes are two dimensions of the [0,1]×[0,1] square region, dimensionless.)

3.3 情形3

與情形1與情形2不同,本情形考慮的是回歸系數(shù)在時空上連續(xù)變化的情況,示意圖見圖3。表3比較了情形3下三種方法的MSEβ。在各種條件下,SCC的MSEβ均為最小。STCC和GTWR的相對優(yōu)劣依賴于時空相關性、樣本數(shù)量和信噪比。在高時空相關性、小樣本和低信噪比條件下,STCC的MSEβ小于GTWR的MSEβ。

表3 情形3下SCC、STCC和GTWR方法的100次數(shù)值模擬實驗的平均MSEβ

(橫軸與縱軸為[0,1]×[0,1]的正方形的兩個維度,無量綱。The horizontal and vertical axes are two dimensions of the [0,1]×[0,1] square region, dimensionless.)

3.4 仿真結果分析

數(shù)值仿真研究結果表明:當回歸系數(shù)在時空上呈現(xiàn)團狀結構時,STCC能夠較好的捕捉到回歸系數(shù)的時空結構,估計誤差比GTWR小得多,也優(yōu)于SCC。當回歸系數(shù)在時空上連續(xù)變化時,SCC的估計效果最優(yōu),STCC和GTWR估計效果整體相當。這說明了STCC具有較強的穩(wěn)健性和自適應能力。

4 海洋時空斷面數(shù)據(jù)分析

4.1 數(shù)據(jù)集介紹

本節(jié)使用STCC方法估計大西洋25°W斷面的溫度-鹽度關系,據(jù)此來檢測南極中層水。南極中層水在南半球的高中緯度地區(qū)形成,在大西洋經向翻轉環(huán)流影響下向北移動,對地球的氣候系統(tǒng)有重要影響。對南極中層水范圍的了解有助于推斷大西洋經向翻轉環(huán)流的強度。與其他水體不同,南極中層水的特點是溫度-鹽度關系為負,因此可以從溫度-鹽度關系中識別南極中層水。溫度和鹽度數(shù)據(jù)從美國國家海洋和大氣管理局發(fā)布的2018年世界海洋圖集(WOA2018)中取得。為了便于分析,本文選取了沿25°W,在60°S和60°N之間的一個溫度和鹽度的經線段進行研究,隨機選取了共10 000個地點的溫度和鹽度測量值。

圖4 25°W斷面每月溫度分布圖Fig.4 Temperature distribution map along 25°W per month

圖5 25°W斷面每月鹽度分布圖Fig.5 Salinity distribution map along 25°W per month

4.2 估計結果及分析

為了檢測溫度-鹽度關系的時空結構,建立如下的回歸模型:

S(si,tj)=T(si,tj)β1(si,tj)+
β0(si,tj)+ε(si,tj)。

(18)

式中:響應變量S(si,tj)表示時空位置(si,tj)的鹽度;T(si,tj)表示溫度;回歸系數(shù)β1(si,tj)衡量了溫度-鹽度關系;β0(si,tj)是截距項。將S與T進行無量綱化,即分別除以自身標準差σS,σT,得到S′與T′,建立回歸關系

S′(si,tj)=T′(si,tj)η1(si,tj)+
η0(si,tj)+ε′(si,tj)。

(19)

不難得出η與β之間的關系為

β1(si,tj)=η1(si,tj)σS/σT,

β0(si,tj)=η0(si,tj)σS。

(20)

使用STCC、SCC和GTWR模型估計出的β1的系數(shù)分布(見圖6、圖7和圖8)存在差異。STCC估計得到的β1在空間上的分布比較光滑,其識別出的負溫度-鹽度關系區(qū)域包括公認的南極中層水的生成地,以及被認為與南極中層水有關的低鹽水舌[9]。盡管SCC識別出的負溫度-鹽度關系區(qū)域與STCC類似,但是其估計得到的β1在某些區(qū)域表現(xiàn)出網(wǎng)格狀結構。這種網(wǎng)格狀結構與實際海洋中的溫度-鹽度關系結構不符,更可能是方法本身的不足所造成的虛假結構。GTWR模型的估計結果有大量零散分布的異常值,與真實海洋明顯不符。因此,相比于GTWR與SCC,STCC估計得到的β1更好地反映了實際海洋中的溫度-鹽度關系。

(橫軸和縱軸分別為無量綱的水平和垂直坐標,右側圖示中數(shù)字為回歸系數(shù)值。The horizontal and vertical axes are nondimensional coordinates of horizontal and vertical. The numbers in the right diagram are regression coefficient values.)

(橫軸和縱軸分別為無量綱的水平和垂直坐標,右側圖示中數(shù)字為回歸系數(shù)值。The horizontal and vertical axes are nondimensional coordinates of horizontal and vertical. The numbers in the right diagram are regression coefficient values.)

5 結語

本文提出了一種新的時空回歸方法,該方法稱為STCC方法。STCC方法可以捕獲回歸系數(shù)中的時空模式,特別是團狀模式?；诤Ｑ蟮臏囟?、鹽度觀測數(shù)據(jù),使用STCC方法檢測了南極中層水的時空位置,并給出了隨月份變化的大西洋25°W斷面海水的溫度-鹽度關系。本研究識別出的負溫度-鹽度關系區(qū)域有助于劃分南極中層水的影響范圍,為海洋的其他相關研究提供參考。