時(shí)空團(tuán)狀系數(shù)回歸模型及其在海洋時(shí)空斷面數(shù)據(jù)中的應(yīng)用*

付昊天, 李芙蓉

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

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

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

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

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

1 STCC模型介紹

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

(1)

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

假定在同一時(shí)空位置僅有一次觀測,則可將式(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)維設(shè)計(jì)矩陣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。因此,無法直接使用普通最小二乘擬合來估計(jì)式(2)中β的值,故需將其正則化。假定時(shí)空域中相鄰位置的β傾向于一致,則可通過最小化以下目標(biāo)函數(shù)來估計(jì)β:

(3)

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

1.1 選擇懲罰矩陣P

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

(4)

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

當(dāng)τ=0時(shí),STCC模型相當(dāng)于使用SCC模型對每個(gè)時(shí)刻的回歸系數(shù)進(jìn)行估計(jì);而當(dāng)τ→∞時(shí),STCC模型相當(dāng)于使用Fused LASSO模型[7]對每個(gè)空間點(diǎn)上的回歸系數(shù)進(jìn)行估計(jì)。對于一般的τ>0,懲罰矩陣P同時(shí)對空間和時(shí)間上的相鄰點(diǎn)的回歸系數(shù)差異進(jìn)行懲罰。為同時(shí)促使回歸系數(shù)在空間和時(shí)間上的同質(zhì)性,形成時(shí)空域上的團(tuán)狀結(jié)構(gòu),我們?nèi)ˇ?1。

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

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

2 計(jì)算方法

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

(5)

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

(6)

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

(7)

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

(8)

第一,對迭代中β子問題進(jìn)行求解。

(9)

(10)

(11)

第二,對迭代中γ子問題進(jìn)行求解。

(12)

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

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

(13)

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

第三,對迭代中α子問題進(jìn)行求解。

(14)

3 數(shù)值仿真研究

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

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

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í)空變化的回歸系數(shù);ε(si,tj)～N(0,σ2)獨(dú)立同分布,令σ=0.1,0.2,0.4,分別對應(yīng)信噪比強(qiáng)度的高、中、低的情況。

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

(16)

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

(17)

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

3.1 情形1

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

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

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

3.2 情形2

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

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

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

3.3 情形3

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

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

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

3.4 仿真結(jié)果分析

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

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

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

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

圖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 估計(jì)結(jié)果及分析

為了檢測溫度-鹽度關(guān)系的時(shí)空結(jié)構(gòu),建立如下的回歸模型:

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

(18)

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

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

(19)

不難得出η與β之間的關(guān)系為

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

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

(20)

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

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

(橫軸和縱軸分別為無量綱的水平和垂直坐標(biāo),右側(cè)圖示中數(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 結(jié)語

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