針對傳統(tǒng)空間Moran’s I只適用于分析單一變量的局限性，提出基于Moran’s I的多變量空間自相關性分析理論.首先，借助傳統(tǒng)空間Moran’s I的向量定義，推導出適用于分析多變量空間自相關性的Moran’s I矩陣，并通過蒙特卡洛法模擬研究Moran’s I矩陣中元素的分布情況，結果顯示：在樣本量較小時，只有非主對角線上的元素服從正態(tài)分布，在樣本量較大時，任一元素都服從正態(tài)分布.基于上述分布結論，可對Moran’s I矩陣中元素進行顯著性檢驗.其次，當空間權重矩陣為正定矩陣時，證明Moran’s I矩陣服從Wishart分布.然后，根據(jù)Moran’s I矩陣的代數(shù)意義提出適用于多變量空間自相關理論的若干綜合評價指標.最后，結合多維空氣污染數(shù)據(jù)進行空間自相關分析.

空間Moran’s I; Moran’s I矩陣; 蒙特卡洛模擬; 正態(tài)分布; Wishart分布; 綜合評價

O211.9

A

0818-12

06.012

0 引言

空間自相關分析理論的發(fā)展來源之一是生物計量學的研究[1-7].Moran[1-2]基于生物現(xiàn)象的空間分析將一維空間概念的相關系數(shù)推廣到二維空間，從而定義了Moran’s I.后來，Geary[3]利用回歸分析的方法并結合Durbin-Watson統(tǒng)計量提出了Geary系數(shù)的概念.在Moran’s I和Geary系數(shù)的基礎上，廣義G統(tǒng)計量[8]又被提出.在3個統(tǒng)計量中，廣義G統(tǒng)計量一般用來檢驗小范圍內(nèi)的局部空間依賴性，存在一定局限性.Geary系數(shù)對鄰近單元之間數(shù)據(jù)變化的絕對差較為敏感[9].而Moran’s I與Geary系數(shù)相比，不易受偏離正態(tài)分布的影響，具有更好的分布特征[10]，因此Moran’s I應用更廣.

空間Moran’s I是一種用于衡量空間自相關性的統(tǒng)計指標，它可以檢測空間數(shù)據(jù)中的空間聚集、空間離散和空間隨機性等特征.在過去的幾十年中，空間Moran’s I已經(jīng)成為空間數(shù)據(jù)分析中的一個重要工具，被廣泛應用于地理信息系統(tǒng)、環(huán)境科學、社會學、經(jīng)濟學等領域的研究中.在地理信息系統(tǒng)領域中，眾多學者利用空間Moran’s I來分析地理數(shù)據(jù)的空間分布、聚集程度和空間自相關性等問題.如：Anselin等[11]介紹了名為GeoDa的開源軟件，它可以用于空間數(shù)據(jù)的可視化和分析，并且可以計算空間Moran’s I. Jiao等[12]利用空間Moran’s I從區(qū)域角度關注城市土地價格的空間自相關和空間分布，并采用基于一種復合距離的空間聚類方法來探究土地價格的聚類特征.Li等[13]利用空間Moran’s I分析了中國城市綠地的空間分布和聚集特征，并探討其影響因素.在環(huán)境科學領域中，眾多學者利用空間Moran’s I來分析環(huán)境數(shù)據(jù)的空間分布、污染源的空間分布和環(huán)境質量的空間變化等問題.如：Guo等[14]利用空間Moran’s I分析了四川省宜賓市重金屬的濃度和來源，并評估土壤環(huán)境質量.Liang等[15]以湖南省漣源市典型煤礦城市表土為研究對象，研究了重金屬的污染程度和空間分布，并結合單因素方差分析，確定了不同土地利用模式的濃度差異有統(tǒng)計學意義.Zhang[16]采用Moran’s I和LISA散點圖分析霧霾污染的空間相關性.在社會學和經(jīng)濟學領域中，利用空間Moran’s I來分析社會經(jīng)濟數(shù)據(jù)的空間分布、貧困率的空間分布和城市化進程的空間分布等問題.如：Anselin[17]概述了空間計量經(jīng)濟學模型規(guī)范的分類法，并為空間回歸的結構形式推導出一系列熟悉的規(guī)范.Rey等[18]介紹了一種基于空間Moran’s I的空間-時間自相關分析方法，用于研究區(qū)域經(jīng)濟系統(tǒng)的時空演化.文獻[19]作為空間計量經(jīng)濟學的經(jīng)典著作，介紹了空間Moran’s I的理論分析和計算方法.此外，王雪青等[20]利用CV、全局Moran’s I、Moran’s" I散點圖、LISA聚類圖等方法探索空間分布特征，提出房地產(chǎn)開發(fā)的優(yōu)化方法.除以上研究外，也有學者對空間Moran’s I的定義、性質和假設檢驗等理論進行研究.如：文獻[21]介紹了基于距離統(tǒng)計量的空間自相關分析方法，包括空間Moran’s I的計算方法和應用.甘茂林等[22]為了減少信息損失，提出了一種適用于空間單元含有多個觀測數(shù)據(jù)的改進的空間Moran’s I，并推導出該指數(shù)的期望和方差，Carrijo等[23]提出了一種適用于小樣本數(shù)據(jù)的修正的空間Moran’s I，并給出了其服從的分布及范圍.

雖然空間相關分析已經(jīng)廣泛應用于各種領域，但大多數(shù)研究者通過研究某一種變量的Moran’s I值的分布[23-24]，并結合假設檢驗來判斷該變量是否在空間上存在聚集現(xiàn)象，進而為空間自回歸模型[25]提供理論前提.然而，空間自回歸模型通常是建立在2種或者多種空間型數(shù)據(jù)的基礎之上，因此討論多變量的空間自相關性是有必要的.不僅如此，研究多變量的空間自相關性能為空間相關性[26]、空間回歸模型[27]、空間插值模型[28]提供理論前提.如對于空間回歸模型而言，如果只考慮單一變量的空間自相關性，并由此選擇空間自回歸模型進行預測[29]，可能結果不盡如人意.事實上，大多數(shù)情況應該根據(jù)2種[30-31]或多種變量是否具有空間自相關性和空間聚集程度的大小來選擇合理策略或模型.針對上述問題，本文提出了用于分析多變量空間自相關性的Moran’s I矩陣，并通過蒙特卡洛法模擬對Moran’s I矩陣中元素的分布進行研究，以便對矩陣中的元素進行顯著性檢驗.然后推導并證明了在一定假設下Moran’s I矩陣的分布情況.然后通過類比傳統(tǒng)空間Moran’s I的理論范圍，給出了Moran’s I矩陣中非主對角線上元素的取值范圍，并基于此范圍給出聚集類型的分類方法.

本文第1節(jié)給出空間Moran’s I的基本理論.第2節(jié)提出適用于多變量空間數(shù)據(jù)的多維Moran’s I，通過將數(shù)據(jù)向量拓展為數(shù)據(jù)矩陣，推導出用于衡量多變量空間數(shù)據(jù)聚集性的Moran’s I矩陣.第3節(jié)在不同樣本量下，利用蒙特卡洛法模擬并檢驗Moran’s I矩陣中元素的分布，并利用相關統(tǒng)計知識，在一定假設下推導出Moran’s I矩陣服從Wishart分布.第4節(jié)提出若干Moran’s I矩陣的評價指標，且具有很直觀的解釋性.第5節(jié)基于2022年11月四川省的空氣污染數(shù)據(jù)做實證分析.

1 傳統(tǒng)空間Moran’s I

空間全局Moran’s I是由Patrick Moran提出的，它衡量的是全局空間的自相關性，檢驗某種屬性或現(xiàn)象在整個空間的聚集效應.但是之后有學者發(fā)現(xiàn)全局Moran’s I只能判斷整個區(qū)域的空間特征，如果存在局部區(qū)域與全局空間特征相反的情況，全局Moran’s I則不再適用.1995年，Anselin[32]提出了空間局部Moran’s I，用于測度每個空間單元與其周圍空間區(qū)域的自相關性.隨后大量學者使用局部Moran’s I進行了應用研究[33-34].衡量空間數(shù)據(jù)之間的相關性依賴于空間權重矩陣，其常見的定義方式有：二進制的鄰接空間權重矩陣、地理空間權重矩陣、經(jīng)濟地理空間權重矩陣[35].假設指定區(qū)域有n個空間單元，其中最常用的是采用二進制的鄰接空間權重矩陣W*=（w*ij）n×n，即w*ij（i，j=1，2，…，n）表示空間單元i與空間單元j之間的依賴程度大小，則w*ij=1，如果空間單元i與空間單元j相鄰;w*ij=0，如果空間單元i與空間單元j不相鄰;wii=0.然而為了統(tǒng)一空間單元的外在影響，常常對空間權重矩陣進行矩陣歸一化，即

wij=wij∑ni=1∑nj=1wij.

通常，將歸一化之后的空間權重矩陣W=（wij）n×n稱為嚴格意義上的空間權重矩陣.

假設指定區(qū)域中n個空間單元的觀測為x=（x1，x2，…，xn）T，結合上述空間權重矩陣，進而空間全局Moran’s I的向量定義如下：

I=zTWzw0，

其中

z=x－s， =1n∑ni=1xi，

s=1n∑ni=1（xi－）2， w0=∑ni=1∑nj=1wij.

空間局部Moran’s I的定義如下：

Ii=zi∑nj=1wijzjw0， i=1，2，…，n，

其中

zi=xi－s， zj=xj－s.

理論上空間全局Moran’s I的取值范圍為[-1，1].當觀測數(shù)據(jù)樣本較大時，空間全局Moran’s I服從正態(tài)分布[36].因此可以通過構造Z（I）統(tǒng)計量[37-38]對其進行顯著性檢驗，空間局部Moran’s I的顯著性檢驗同理.

若空間Moran’s I在一定的顯著性水平下通過了檢驗，則可依據(jù)空間Moran’s I的取值范圍來對指定區(qū)域的聚集情況進行分類，具體如下：Igt;0，說明整個空間區(qū)域存在空間正相關性，即該區(qū)域呈現(xiàn)的是擁有高（低）觀測值的單元被擁有高（低）觀測值的單元所包圍的聚集效應.Ilt;0，說明整個空間區(qū)域存在空間負相關性，即該空間區(qū)域呈現(xiàn)的是擁有高（低）觀測值的單元被擁有低（高）觀測值的單元所包圍的聚集效應.I=0，表示整個空間區(qū)域不存在空間相關性，即各空間單元的觀測值在整個空間的分布呈現(xiàn)隨機性.

2 Moran’s I矩陣

經(jīng)研究發(fā)現(xiàn)：空間單元中某一變量的聚集性變化會受到臨近單元其他變量的影響[39-40].而且，當對某變量進行空間自回歸建模時，如果不考慮多種變量的聚集情況，那么預測的結果將不夠理想.因此，單變量的Moran’s I對于多變量的數(shù)據(jù)已不再適用，而多變量的空間自相關性的研究應運而生.

假設指定區(qū)域中有n個空間單元，每個空間單元包含p個變量，那么觀測數(shù)據(jù)矩陣為

X=（x（1），x（2），…，x（p））=

x11x12…x1p

x21x22…x2p

xn1xn2…xnp

，

其中，x（i）=（x1i，x2i，…，xni）T表示第i個變量在n個空間單元中的觀測值，i=1，2，…，p.進而，列標準化后的觀測數(shù)據(jù)矩陣為

Z=（Z（1），Z（2），…，Z（p）），

其中

Z（i）=x（i）－（i）S（i）， （i）=1n∑nj=1xji，

S（i）=1n∑nj=1（xji－（i））2， i=1，2，…，p，

（i）表示第i個變量的平均值，S（i）表示第i個變量的標準差.

根據(jù)單變量的空間全局Moran’s I的向量式定義，多變量的空間全局Moran’s I的向量式定義為

I*p×p=ZTWZw0=

（Z（1），Z（2），…，Z（p））TW（Z（1），Z（2），…，Z（p））w0=

Z（1）TWZ（1）w0Z（1）TWZ（2）w0…Z（1）TWZ（p）w0

Z（2）TWZ（1）w0Z（2）TWZ（2）w0…Z（2）TWZ（p）w0

Z（p）TWZ（1）w0Z（p）TWZ（2）w0…Z（p）TWZ（p）w0

p×p=

I*11I*12…I*1pI*21I*22…I*2p

I*p1I*p2…I*pp

p×p，

其中

I*ij=Z（i）TMZ（j）w0， 1≤i，j≤p.

在矩陣I*p×p中，當i≠j時，I*ij表示不同變量之間空間全局Moran’s I;當i=j時，I*ii表示同一變量空間全局Moran’s I.I*ij表示空間單元中變量i與周圍空間單元中變量j之間的聚集情況.顯然當p=1時，I*p×p就退化為單變量情況的空間全局Moran’s I.因此，傳統(tǒng)Moran’s I只是Moran’s I矩陣在一維情況下的特殊情況，它們的本質區(qū)別主要在于考慮的空間維度和數(shù)據(jù)結構.

在應用方面，傳統(tǒng)Moran’s I是一種用于測量單一變量在空間上的自相關性的方法.它通常用于分析單一屬性或變量在地理空間上的分布模式，例如人口密度、降雨量等.Moran’s I矩陣則是針對多維數(shù)據(jù)（即多個相關變量）在空間上的自相關性進行建模和測量的方法.在現(xiàn)實世界中，很多情況下需要同時考慮多個屬性之間的空間分布模式，例如城市發(fā)展中同時考慮人口、經(jīng)濟、環(huán)境等多個方面的數(shù)據(jù).Moran’s I矩陣考慮了多個屬性之間的關聯(lián)性，能夠更全面地揭示多個相關變量之間的空間關系.

由于I*ij（i≠j）的取值范圍的不確定性，導致無法根據(jù)該指數(shù)的取值來確定指定區(qū)域所屬的聚集類型，下面借助方差非負性的知識，給出非主對角線上元素取值范圍的定理.

定理 1 理論上，I*ij（i≠j）的取值范圍為

I*ij∈[－w20+n22w20，w20+n22w20].

當對空間權重矩陣進行歸一化時，即w0=1，則I*ij∈[－n22，n22].

當對空間權重矩陣進行列歸一化時，即w0=n，則I*ij∈[－1，1].

證明

假設指定區(qū)域中第i個變量在第k個空間單元觀測數(shù)據(jù)為xki，k=1，2，…，n，i=1，2，…，p，并且同一個變量的觀測都是獨立同分布的隨機變量，即

xki～N（ui，σ2i）， i=1，2，…，p.

令

Yi=xki－（i）， Yj=xkj－（j），

那么

cov（Yi，Yj）=1n∑nk=1（xki－（i））（xkj－（i））

表示隨機變量Yi與隨機變量Yj之間的協(xié)方差，而

var（Yi）=1n∑mk=1（xki－）2=σ2i

表示隨機變量Yi的方差，

進而，I*ij（i≠j）等價于

I*ij=Z（i）TWZ（j）w0=n∑ni=1（xki－（i））∑nj=1wij（xkj－（j））w0[∑ni=1（xki－（i））2]12[∑ni=1（xkj－（j））2]12=ncov（Yi，Yj）w0σiσj.

又由于

var（Yiσi+nYjw0σj）≥0，

即

1σ2ivar（Yi）+n2w20σ2jvar（Yj）+2ncov（Yi，Yj）w0σiσj≥0，

1+n2w20+2I*ij≥0，

I*ij≥－w20+n22w20.

同理

var（Yiσi－nYjw0σj）≥0，

則

I*ij≤w20+n22w20，

即

I*ij∈[－w20+n22w20，w20+n22w20].

若空間權重矩陣進行歸一化，即使得空間權重矩陣的所有元素之和等于1，則有w0=1，則

I*ij∈[－n22，n22].

若對空間權重矩陣進行列歸一化，即使得空間權重矩陣的每列之和等于1，則有w0=n，則

I*ij∈[－1，1].

由于定理1給出了I*ij（i≠j）的取值范圍.若非主對角線上的空間Moran’s I在一定的顯著性水平下通過了檢驗，那么類比傳統(tǒng)空間Moran’s I的分類方法，則基于I*ij的分類方法如下：I*ijgt;0，說明整個空間區(qū)域存在空間正相關性，即該區(qū)域呈現(xiàn)的是擁有高（低）觀測值的單元被擁有高（低）觀測值的單元所包圍的聚集效應.I*ijlt;0，說明整個空間區(qū)域存在空間負相關性，即該空間區(qū)域呈現(xiàn)的是擁有高（低）觀測值的單元被擁有低（高）觀測值的單元所包圍的聚集效應.I*ij=0，表示整個空間區(qū)域不存在空間相關性，即各空間單元的觀測值在整個空間的分布呈現(xiàn)隨機性，且I*ij越接近范圍的邊界值，表明空間單元的聚集效應越強，反之則越弱.

然而，對于空間全局Moran’s I而言，顯著性檢驗離不開該指數(shù)的分布，因為當觀測樣本量較大時，并且觀測數(shù)據(jù)獨立且服從同一正態(tài)分布，空間全局Moran’s I服從的是正態(tài)分布[36]，即Moran’s I矩陣主對角線上的元素服從正態(tài)分布，所以才能夠對其進行顯著性檢驗，因此研究Moran’s I矩陣及矩陣中元素的分布是尤為重要的.

3 Moran’s I矩陣元素及矩陣分布

傳統(tǒng)空間Moran’s I服從正態(tài)分布已有理論上的證明，但由于Moran’s I矩陣中非主對角線上元素分布的證明相對復雜，因此本文采用計算機模擬的方式來驗證Moran’s I矩陣中非主對角線上元素的分布.區(qū)別在于理論證明過程的核心是利用大數(shù)定律、中心極限定理等統(tǒng)計定理進行證明，而計算機模擬利用計算機實現(xiàn)不同樣本量下的頻數(shù)統(tǒng)計及可視化的能力.而聯(lián)系在于理論證明作為計算機模擬驗證的先驗知識，才有了進一步的研究可能.

下面將利用蒙特卡洛法模擬[41]研究Moran’s I矩陣中元素I*ij（1≤i，j≤p）服從的分布，進而實現(xiàn)對Moran’s I矩陣中元素的假設檢驗.下面假設有一個包含（k×k）個空間單元的規(guī)則系統(tǒng)，即指定系統(tǒng)中每個空間單元是矩形區(qū)域，整個系統(tǒng)由k行k列的空間單元組成.若用p代表每個空間單元觀測數(shù)據(jù)的維數(shù)，則觀測數(shù)據(jù)矩陣的維數(shù)為（k×k）×p.由此可以根據(jù)如下算法來研究Moran’s I矩陣中元素的分布情況，詳細的算法步驟如下：

算法 Moran’s I矩陣元素的模特卡洛模擬：

1） （k×k）×p維的觀測數(shù)據(jù)矩陣分別由任意p個正態(tài)分布隨機產(chǎn)生;

2） 觀測數(shù)據(jù)被隨機分配至（k×k）個空間單元中;

3） 計算Moran’s I矩陣;

4） 步驟1～3重復10 000次.

下面，令p=4，k=6，20，100，并分別用X1、X2、X3、X4代表4維的觀測變量.采用基于Queen鄰接[42-44]的空間權重矩陣，則可利用R語言畫出Moran’s I矩陣元素的經(jīng)驗分布圖，并且每一個圖按照Moran’s I矩陣的行、列進行排列，其中，橫坐標表示改進的Moran’s I的值，縱坐標表示出現(xiàn)的頻數(shù)，詳細結果見圖1～3（主對角線上表示單變量Moran’s I的經(jīng)驗分布圖，其他表示雙變量Moran’s I的經(jīng)驗分布圖）.另外，為了使驗證結果更加可靠，本文基于SW統(tǒng)計量對模擬結果進行檢驗，詳細結果見表1～3.

根據(jù)表1～3顯示，當樣本量較小時，Moran’s I矩陣中主對角線上的元素不服從正態(tài)分布，而非主對角線上的元素服從正態(tài)分布.當樣本量較大時，Moran’s I矩陣中所有元素都服從正態(tài)分布.因此，對于主對角線上的改進的空間全局Moran’s I，在樣本量較大時才能用于衡量指定區(qū)域中是否存在聚集.而非主對角線上改進的空間全局Moran’s I適用于任何情況.

需要注意的是：Queen鄰接的空間權重矩陣具有二元連接關系、對稱性、不考慮距離等特點，并且本節(jié)主要目的在于驗證Moran’s I矩陣中的元素僅在大樣本觀測下服從正態(tài)分布，采用Queen鄰接的空間權重矩陣是因為其具有一般性的特點，同時參考了文獻[33]中的模擬方法.因此，如果采用基于Rook鄰接的空間權重矩陣、地理空間權重矩陣或經(jīng)濟地理空間權重矩陣，Moran’s I矩陣中的元素分布不變，具體實現(xiàn)見第5節(jié)中增加實驗數(shù)據(jù)以及詳細代碼的獲取方式.

另外，根據(jù)多元統(tǒng)計分析理論知，

下面，用x（·）表示資料矩陣X的行向量，并假設x（·）獨立且服從同一p元正態(tài)分布，即

x（·）～Np（μ，Σ），

并且每一次觀測向量相互獨立.設Z（·）為列標準化后的觀測數(shù)據(jù)矩陣Z中的行向量，并根據(jù)Z與觀測數(shù)據(jù)矩陣X之間的關系，進而向量Z（·）與x（·）有關系

ZT（·）=AxT（·）-B，

其中

A=

S（1）0…0

0S（2）…0

00…S（p）

－1p×p，

B=

（1）S（i）（2）S（2）…（p）S（p）

Tp×1.

進而，根據(jù)多分正態(tài)分布的性質可知，隨機向量的任何線性變換仍然服從正態(tài)分布，即

Z（·）～Np（0，Σ*）， Σ*=AΣAT，

且彼此相互獨立.

由于Moran’s I矩陣作為衡量多個變量之間空間自相關性的指標，為了進一步研究該指標的實際意義，通常需要研究該指標的分布情況，下面假設空間權重矩陣為正定矩陣，那么Moran’s I矩陣服從Wishart分布[45-46]，詳細證明見定理2.

定理 2 若假設空間權重矩陣為正定矩陣，那么空間權重矩陣W為對稱且正定的矩陣，則Moran’s I矩陣服從Wishart分布.

證明

在上述假設條件下，根據(jù)隨機矩陣的定義，去中心化后資料矩陣的分布為

Z～Nn×p（0，Σ*In），

其中，表示Kronecker乘積.

又由W為正定且對稱的矩陣，則存在正交矩陣U，使得

UTWU=

λ10…0

0λ2…0

00…λn

n×n，

進而，存在可逆矩陣P，使得

PTUTWUP=In，

其中，In為n階單位矩陣，且

P=

1λ10…0

01λ2…0

00…1λn

n×n，

即

QTWQ=In，

其中，Q=UP.

進而

W=（QT）－1InQ－1=（Q－1）TQ－1，

即

I*p×p=1w0ZT（Q－1）TQ－1Z=1w0RTR，

其中，R=Q－1Z，且

R～Nn×p（0，（Q－1）TΣ*Q－1In）.

故根據(jù)Wishart分布的定義知

I*p×p=1w0RTR～Wp（n，（Q－1）TΣ*Q－1）.

需要注意的是：Moran’s I矩陣的分布取決于空間權重矩陣是不是對稱正定.若空間權重矩陣是對稱正定，則存在正交矩陣U，使得空間權重矩陣W可對角化，即

UTWU=

λ10…0

0λ2…0

00…λn

n×n，

進而決定著Moran’s I矩陣是否可以表示為2個矩陣的乘積（Wishart分布的定義）.

以本文提到的空間權重矩陣為例，因為二進制的鄰接空間權重矩陣、地理空間權重矩陣本身是實對稱的空間權重矩陣，進而若滿足正定的條件，那么Moran’s I矩陣便符合定理2的結論.另外，因為經(jīng)濟地理空間權重矩陣不一定為對稱矩陣，則對基于經(jīng)濟地理空間權重矩陣的Moran’s I矩陣而言，定理2便不一定滿足.

4 Moran’s I矩陣評價

在此之前，已經(jīng)通過蒙特卡洛法證明了Moran’s I矩陣中元素服從正態(tài)分布.那么便可以通過經(jīng)典的假設檢驗方法，即構造Z（I）統(tǒng)計量對Moran’s I矩陣元素進行顯著性檢驗.本文設檢驗后的Moran’s I矩陣為

Ip×p=（Iij）p×p，

該矩陣僅保留通過顯著性檢驗的全局空間Moran’s I，即

Iij=

I*ij， I*ij通過顯著性檢驗;

0， I*ij未通過顯著性檢驗.

然后根據(jù)檢驗后的Moran’s I矩陣Ip×p的代數(shù)意義及應用場景，提出如下適用于空間自相關理論的若干綜合評價指標.

4.1 Ip×p的M范數(shù)

矩陣絕對值最大的元素取絕對值，即

‖Ip×p‖M=max{|Iij|，1≤i，j≤p}.

該范數(shù)可用于判斷首要聚集特征.當Iijgt;0，表示首要聚集特征呈現(xiàn)高高、低低聚集現(xiàn)象;當Iijlt;0，表示首要聚集特征呈現(xiàn)高低聚集現(xiàn)象;當Iij=0，表示無首要聚集特征，且所有特征沒有聚集現(xiàn)象.

4.2 Ip×p的1-范數(shù)

矩陣列向量中絕對值之和的最大值，即

‖Ip×p‖1=max{∑pi=1|Iij|，j=1，2，…，p}.

該范數(shù)用于判斷“綜合首要”聚集特征，“綜合首要”即某特征和其他所有特征存在最多的聚集現(xiàn)象.通常“綜合首要”聚集特征在多變量特征中具有主成分的含義，即對綜合聚集現(xiàn)象的貢獻最大.

4.3 Ip×p的跡

主對角線元素求和，即

tr（Ip×p）=∑pIii.

該指標用于判斷所有特征是否存在綜合聚集現(xiàn)象.當tr（Ip×p）gt;0，則表明綜合考慮，特征之間存在高高、低低聚集現(xiàn)象（此處不能根據(jù)tr（Ip×p）判斷2種特征之間是否存在聚集現(xiàn)象，更不能判斷某地區(qū)的2種特征之間是否存在聚集現(xiàn)象）.當tr（Ip×p）lt;0，則表明綜合考慮，特征之間存在高低聚集現(xiàn)象.當tr（Ip×p）=0，則表明綜合考慮，特征之間不存在聚集現(xiàn)象.

4.4 Ip×p的L0范數(shù)

非零元素的個數(shù).可根據(jù)r（非零元素的個數(shù)占矩陣維數(shù)的百分比）判斷聚集等級.具體參考標準正態(tài)分布的一倍標準差占比68.3%，則聚集等級為三級（一般聚集）;二倍標準差占比95.5%，則聚集等級為二級（較聚集）;三倍標準差占比97.7%，則聚集等級為一級（顯著聚集），具體見表4.

5 案例研究

本節(jié)以2022年11月四川省為研究區(qū)域，該區(qū)域包含21個市（州），將每個空間單元的空氣污染數(shù)據(jù)用于實證分析.涉及圖件文本數(shù)據(jù)包括：1） 四川省2015年縣界限Shp格式矢量圖與地形高程圖.2） 各市（州）空氣污染物數(shù)據(jù)來源于四川省生態(tài)環(huán)境廳數(shù)據(jù)庫（http：//sthjt.sc.gov.cn/sthjt/sjxz/list_w.shtml），包括O3、PM10、PM2.5和SO2的觀測濃度值.

（，具體濃度分布見圖4，和圖同時刪）

軟件使用基于R語言的Rstudio，詳細代碼見個人開源的github代碼倉庫（https：//github.com/Spatsta/TradiSpa.git）.

表5顯示了2022年11月四川省關于多種空氣污染物的空間全局Moran’s I以及對應p值.由

于本次案例的實驗數(shù)據(jù)僅僅包含21個空間單元，因此對于空間全局Moran’s I的顯著性檢驗不能通過構造Z（I）統(tǒng)計量實現(xiàn)，而是通過隨機置換檢驗[47]這種更一般的檢驗方式開展的.

表6顯示了2022年11月四川省21個市（州）關于各變量檢驗后的Moran’s I矩陣Ip×p.

由表6可以看出，單變量空間全局Moran’s I都為正且通過顯著性檢驗，表明O3、PM10、PM2.5和SO2這4種變量在四川省21個市（州）具有正的空間自相關，即該區(qū)域呈現(xiàn)的是擁有高（低）觀測值的單元被擁有高（低）觀測值的單元所包圍的聚集效應.雙變量空間全局Moran’s I只有SO2分別與O3、PM10和PM2.5未通過顯著性檢驗，其余雙變量空間全局Moran’s I都為正且通過了顯著性檢驗.

圖5非常直觀地展示了變量之間的聚集現(xiàn)象.

另外，對于Moran’s I矩陣整體而言，顯然，大部分都通過了顯著性檢驗.

接下來，根據(jù)上一節(jié)中空間自相關理論的綜合評價指標，對多變量空間自相關分析進行進一步探索，結果見表7.

為了更好地了解不同變量之間的聚集情況以及各污染物對空氣質量指數(shù)AQI[48]的影響，圖5直觀地展示了相關信息.

由表7知，Moran’s I矩陣的M范數(shù)為0.510gt;0，即為O3的空間全局Moran’s I的絕對值，因此首要聚集特征為O3.首要聚集特征表明某地區(qū)聚集現(xiàn)象最顯著的變量，在進行污染防控時，該污染物應給予高度的重視.根據(jù)表5可知O3具有正的空間自相關性，即表現(xiàn)出高-高型、低-低型聚集現(xiàn)象.Moran’s I矩陣的1范數(shù)為0.893gt;0，即為變量PM2.5所在列，因此“綜合首要”聚集特征為PM2.5.“綜合首要”聚集特征表明了某地區(qū)中與其他污染物聚集“最多”的變量，也即最容易與其他污染物形成聚集的變量，該特征一般在綜合指標中占主要地位.

該特征在空氣質量指標AQI中占主導地位，這一點可以很明顯的從圖5中觀察到.

Moran’s I矩陣的跡為1.386gt;0，表明案例中所有污染物總的來看具有正的空間自相關，即上述4種污染物之間一定存在聚集性.若要研究他們之間的關聯(lián)或影響，綜合考慮每項污染物的影響是合理的.Moran’s I矩陣的L0范數(shù)為10，且占比為r=10/16=0.625，查表4可知，該地區(qū)的聚集等級為“一般聚集”.

6 結論

本文的性質屬于空間自相關理論及方法的拓展研究，主要目的是探究多變量空間自相關性的研究方法.主要拓展了以下內(nèi)容：

1） 提出并推導了適用于多變量空間自相關分析的Moran’s I矩陣，該矩陣在空間上包含了更多聚集信息，彌補了單變量空間Moran’s I不能解釋多變量在空間中存在的關聯(lián)聚集，同時也能作為選擇空間回歸模型、空間插值模型等的理論前提.并且，在理論上Moran’s I矩陣中任一元素的取值范圍為[-1，1].

2） 模擬并檢驗了在不同樣本量下，Moran’s I矩陣中元素的分布情況，并在一定條件下論證了Moran’s I矩陣服從Wishart分布.雖然本文推導出了Moran’s I矩陣的分布情況，但Moran’s I矩陣整體的檢驗統(tǒng)計量仍需進一步研究.

3） 提出了基于Moran’s I矩陣若干綜合聚集指標，而且在實際應用中具有很直觀的解釋性，豐富了空間自相關理論的實際意義.

參考文獻

[1] MORAN P A P. The interpretation of statistical maps[J]. Journal of the Royal Statistical Society： Series B（Methodological），1948，10（2）：243-251.

[2] MORAN P A P. Notes on continuous stochastic phenomena[J]. Biometrika，1950，37（1/2）：17-23.

[3] GEARY R C. The contiguity ratio and statistical mapping[J]. The Incorporated Statistician，1954，5（3）：115-146.

[4] SOKAL R R， ODEN N L. Spatial autocorrelation in biology： 1 methodology[J]. Biological Journal of the Linnean Society，1978，10（2）：199-228.

[5] CLIFF A D， ORD J K. The problem of spatial autocorrelation[M]. London： Pion，1969.

[6] CLIFF A D， ORDJ K. Spatial autocorrelation[J]. Journal of Economic Literature，1976，14（3）：924-925.

[7] CLIFF A D， ORD J K. Spatial processes： models and applications[M]. London： Pion，1981.

[8] SERFLING R J. Generalized l-， m-， and r-statistics[J]. The Annals of Statistics，1984，12（1）：76-86.

[9]" 戴芹，馬建文，陳雪. 北京環(huán)線建設驅動的土地利用變化遙感檢測與分析[J]. 遙感學報，2005，9（3）：314-322.

[10]" 王遠飛，何洪林. 空間數(shù)據(jù)分析方法[M]. 北京：科學出版社，2007.

[11] ANSELIN L， SYABRI I， KHO Y. GeoDa： an introduction to spatial data analysis[J]. Geographical Analysis，2009，38（1）：73-89.

[12] JIAO L M， LIU Y L. Analyzing the spatial autocorrelation of regional urban datum land price[J]. Geo-Spatial Information Science，2012，15（4）：263-269.

[13] LI F， SUN W W， YANG G， et al. Investigating spatiotemporal patterns of surface urban heat islands in the Hangzhou metropolitan area， China， 2000—2015[J]. Remote Sensing，2019，11（13）：1553.

[14] GUO G H， WU F C， XIE F Z， et al. Spatial distribution and pollution assessment of heavy metals in urban soils from southwest China[J]. Journal of Environmental Sciences，2012，24（3）：410-418.

[15] LIANG J， FENG C T， ZENG G M， et al. Spatial distribution and source identification of heavy metals in surface soils in a typical coal mine city， Lianyuan， China[J]. Environmental Pollution，2017，225（17）：681-690.

[16] ZHANG Y Y. Analysis of the spatial effects of inter-provincial air pollution in China[J]. E3S Web of Conferences，2020，194（1）：1-4.

[17] ANSELIN L. Spatial externalities， spatial multipliers， and spatial econometrics[J]. International Regional Science Review，2003，26（2）：153-166.

[18] REY S J， JANIKAS M V. STARS： space-time analysis of regional systems [J]. Geographical Analysis，2006，38（1）：67-86.

[19] LESAGE J， PACE R K. Introduction to spatial econometrics[M]. New York： Chapman and Hall，2009.

[20] 王雪青，陳媛，劉炳勝. 中國區(qū)域房地產(chǎn)經(jīng)濟發(fā)展水平空間統(tǒng)計分析：全局Moran’s I、Moran散點圖與LISA集聚圖的組合研究[J]. 數(shù)理統(tǒng)計與管理，2014，33（1）：59-71.

[21] GETIS A， ORD J K. The analysis of spatial association by use of distance statistics[J]. Geographical Analysis，1992，24（3）：189-206.

[22] 甘茂林，呂王勇，符璐. 基于改進Moran’s I指數(shù)的成都市PM2.5的空間統(tǒng)計分析[J]. 環(huán)境科學與技術，2016，39（9）：187-193.

[23] CARRIJO T B， DA SILVA A R. Modified Moran’s I for small samples[J]. Geographical Analysis，2017，49（4）：451-467.

[24] DIAWARA N， WALLER L， KING R， et al. Simulations of local Moran’s Index in a spatio-temporal setting[J]. Communications in Statistics-Simulation and Computation，2019，48（6）：1849-1859.

[25] 呂燁. 空間自回歸模型的穩(wěn)健估計[D]. 昆明：云南財經(jīng)大學，2018.

[26] 姜磊. 空間回歸模型選擇的反思[J]. 統(tǒng)計與信息論壇，2016，31（10）：10-16.

[27] 朱會義，劉述林，賈紹鳳. 自然地理要素空間插值的幾個問題[J]. 地理研究，2004，23（4）：425-432.

[28] 何紅艷，郭志華，肖文發(fā). 降水空間插值技術的研究進展[J]. 生態(tài)學雜志，2005，24（10）：1187-1191.

[29] 趙永，劉旭華，孫騰達. 基于空間自回歸模型的中國耕地面積變化預測[J]. 干旱區(qū)資源與環(huán)境，2013，27（8）：1-5.

[30] LEVINE N. Crime mapping and the crimestat program[J]. Geographical Analysis，2006，38（1）：41-56.

[31] XUE B， XIAO X， LI J Z. Identification method and empirical study of urban industrial spatial relationship based on POI big data： a case of Shenyang City， China[J]. Geography and Sustainability，2020，1（2）：152-162.

[32] ANSELIN L. Local indicators of spatial association： LISA[J]. Geographical Analysis，1995，27（2）：93-115.

[33] 丁偉，高大帥，駱華松，等. “一帶一路”沿線國家或地區(qū)人類發(fā)展水平空間分異研究[J]. 資源開發(fā)與市場，2020，36（11）：1219-1226.

[34] BARBOSA C C， BONFIM C V D， DE BRITO C M G， et al. Spatial analysis of reported new cases and local risk of leprosy in hyper-endemic situation in northeastern Brazil [J]. Tropical Medicine amp; International Health，2018，23（7）：748-757.

[35] 王薇，劉韞寧，殷鵬，等. 不同空間權重矩陣對我國心血管疾病死亡空間自相關分析的影響[J]. 中華流行病學雜志，2021，42（8）：1437-1444.

[36] CLIFF A D， ORD J K. Classics in human geography revisited[J]. Progress in Human Geography，1995，19（2）：245-249.

[37] 陶吉興，傅偉軍，姜培坤，等. 基于Moran’s I 和地統(tǒng)計學的浙江森林土壤有機碳空間分布研究[J]. 南京林業(yè)大學學報（自然科學版），2014，38（5）：97-101.

[38] JIN F， LEE L F. On the bootstrap for Moran’s I test for spatial dependence[J]. Journal of Econometrics，2015，184（2）：295-314.

[39] 申婷，李飛雪，陳振杰. 基于多源數(shù)據(jù)的城市活力評價與空間關聯(lián)性分析：以常州市主城區(qū)為例[J]. 長江流域資源與環(huán)境，2022，31（5）：1006-1015.

[40] 馮軍，陳太好，戴麗芳，等. 貴州省2020—2021年流感季6月齡～5歲兒童流感疫苗接種率與社會經(jīng)濟指標的雙變量空間自相關分析[J]. 中國疫苗和免疫，2022，28（2）：199-203.

[41] 邵偉. 蒙特卡洛方法及在一些統(tǒng)計模型中的應用[D]. 濟南：山東大學，2012.

[42] 羅正鋒. 基于空間面板模型的我國主要城市房價影響因素分析[D]. 成都：西南財經(jīng)大學，2020.

[43] 楊書，范博凱，顧蕓. 投資型環(huán)境規(guī)制對綠色全要素生產(chǎn)率的非線性影響[J]. 中國人口·資源與環(huán)境，2022，32（5）：120-131.

[44] WANG P R， ZENG C， SONG Y， et al.The spatial effect of administrative division on land-use intensity[J]. Land，2021，10（5）：543.

[45] OUIMET F. A symmetric matrix-variate normal local approximation for the Wishart distribution and some applications[J].Journal of Multivariate Analysis，2022，189（1）：1-17.

[46] CARO-LOPERA F J， GONZLEZ FARAS G， BALAKRISHNAN N. Matrix variate distribution theory under elliptical models-V： the non-central Wishart and inverted Wishart distributions[J]. Mathematical Methods of Statistics，2022，31（1）：18-42.

[47] JAYA I G N M， ANDRIYANA Y， TANTULAR B， et al. Spatiotemporal dengue disease clustering by means local spatiotemporal Moran’s Index[J]. IOP Conference Series： Materials Science and Engineering，2019，621（1）：1-17.

[48] FANG， CL， LIU， HM; LI， GD， SUN， DQ， MIAO， Z. Estimating the Impact of Urbanization on Air Quality in China Using Spatial Regression Models.[J].Sustainability，2015，7（11）：15570-15592.

Research and Application of Multivariate Spatial Autocorrelation Based on Moran’s I

ZHANG Ce， L Wangyong， ZHANG Ping， SONG Jiacheng

（School of Mathematical Sciences， Sichuan Normal University， Chengdu 610066， Sichuan）

In response to the limitations of traditional spatial Moran’s I which only applies to the analysis of a single variable， this paper proposes a multivariate spatial autocorrelation analysis theory based on Moran’s I. Firstly， utilizing the vector definition of traditional spatial Moran’s I， a Moran’s I matrix suitable for analyzing multivariate spatial aggregation is derived. And through Monte Carlo simulation， the distribution of elements in the Moran’s I matrix is studied. The results indicate that only off-diagonal elements follow a normal distribution when the sample size is small. However， all elements follow a normal distribution when the sample size is large. Therefore the elements in Moran’s I matrix can be tested for significance. Secondly， when the spatial weight matrix is a positive definite matrix， the Moran’s I matrix follows a Wishart distribution. Thirdly， several comprehensive evaluation indicators applicable to the multivariate spatial autocorrelation theory are proposed based on the algebraic significance of the Moran’s I matrix. Finally， spatial autocorrelation study based on multidimensional air pollution data was carried out.

spatial Moran’s I; Moran’s I matrix; Monte Carlo simulation; normal distribution; Wishart distribution; comprehensive evaluation

2020 MSC：62P12

（編輯 周 ?。?/p>