matlab經(jīng)驗(yàn)正交函數(shù)EOF(轉(zhuǎn)載)

1、A.7 EOF分析經(jīng)驗(yàn)正交函數(shù)分析方法(empirical orthogonal function.縮寫(xiě)為EOF),也稱特征 向量分析(eigenvector analysis),或者主成分分析(principal component analysis.縮 寫(xiě)PCA),是一種分析矩陣數(shù)據(jù)中的結(jié)構(gòu)特征，提取主要數(shù)據(jù)特征量的一種方 法。Lorenz在195()年代首次將其引入氣象和氣候研究，現(xiàn)在在地學(xué)及其他學(xué)科中 得到了非常廣泛的應(yīng)用。地學(xué)數(shù)據(jù)分析中通常特征向量對(duì)應(yīng)的是空間樣本，所以 也稱空間特征向量或者空間模態(tài)；主成分對(duì)應(yīng)的是時(shí)間變化，也稱時(shí)間系數(shù)。因 此地學(xué)中也將EOF分析稱為時(shí)空分解。原理與

2、算法選定要分析的數(shù)據(jù)，進(jìn)行數(shù)據(jù)預(yù)處理，通常處理成距平的形式。得到一個(gè)數(shù) 據(jù)矩陣乙心“計(jì)算X與其轉(zhuǎn)置矩陣XT的交叉積，得到方陣cmxm = -XxXTn如杲X是已經(jīng)處理成了距平的話，則C稱為協(xié)方差陣；如杲X已經(jīng)標(biāo)準(zhǔn) 化(即C中每行數(shù)據(jù)的平均值為0,標(biāo)準(zhǔn)差為1),則C稱為相關(guān)系數(shù)陣計(jì)算方陣C的特征根(兒凹)和特征向量二者滿足Xx JH Vr/iX/M X其中A是mxm維對(duì)角陣，即'Ai 0.00 A2 .0A =00 . Ani 一般將特征根入按從大到小順序排列，即入1 入2 . 丸“。因?yàn)閿?shù) 據(jù)X是真實(shí)的觀測(cè)值，所以入應(yīng)該大于或者等于()每個(gè)非()的特征根對(duì)應(yīng) 一列特征向量值，也稱EOF

3、。如入1對(duì)應(yīng)的特征向量值稱第一個(gè)EOF模態(tài)， 也就是U的第一列即EOFi = V(:.l);第九對(duì)應(yīng)的特征向量是U的第斤列， 即 EOFk = Vk).計(jì)算主成分。將EOF投影到原始資料矩陣X上，就得到所有空間特征向量對(duì) 應(yīng)的時(shí)間系數(shù)(即主成分)，即其中PC中每行數(shù)據(jù)就是對(duì)應(yīng)每個(gè)特征向量的時(shí)間系數(shù)。第一行PC(1,:)就是 第一個(gè)EOF的時(shí)間系數(shù)，其他類(lèi)推。上面是對(duì)數(shù)據(jù)矩陣X進(jìn)行計(jì)算得到的EOF和主成分(PC),因此利用EOF和PC也 可以完全恢復(fù)原來(lái)的數(shù)據(jù)矩陣X,即X = EOF x PC有時(shí)可以用前面最突出的幾個(gè)EOF模態(tài)就可以擬合出矩陣X的主要特征。此 外，EOF和PC都具有正交性的特點(diǎn)

4、，可以證明x PCT = A；即不同的PC之 間相關(guān)為0。E x ET = I. I為對(duì)角單位矩陣，即對(duì)用線上值為1,其他元素都 為0。這表明各個(gè)模態(tài)之間相關(guān)為0,是獨(dú)立的°由上面的計(jì)算過(guò)程可以看出，EOF分析的核心是計(jì)算矩陣C的特征根和特征向 量。計(jì)算矩陣特征根和特征向量的方法很多，下面具體給出Matlab中進(jìn)行EOF分 析的兩種不同的方法.具體步驟可參考下面兩個(gè)框圖中的實(shí)例。方法1:調(diào)用rEOF,El=eig(C),其中EOF為計(jì)算得到的空間特征向號(hào)，E為特 征根。然后計(jì)算主成分PC = EOFT x 需要指出的時(shí)，當(dāng)數(shù)據(jù)量很大時(shí)，例 如分析高分辨率的資料(如lkm分辨率的NDV

5、I資料)，空間范圍很大維數(shù)m很容易 超過(guò)數(shù)萬(wàn)個(gè)點(diǎn)，則矩陣C的維數(shù)是個(gè)巨大量，需要占用大量?jī)?nèi)存，也會(huì)導(dǎo)致計(jì)算 速度異常緩慢。而且很可能超出計(jì)算機(jī)的計(jì)算圾限而死機(jī)。方法2:直接對(duì)矩陣X進(jìn)行奇異值分解X =吃廠其中刀為奇異值對(duì)交陣(刀對(duì)角線上的元素為奇異值)，奇異值與特征根成倍數(shù)關(guān) 系.如果矩陣C=XXT, C的特征根為入，則有刀=応；如果矩陣c = xxT, c的特征根為入，則有e = yx；由于該方法是直接對(duì)矩陣x進(jìn)行分解，所以對(duì)內(nèi)存的要求遠(yuǎn)小于方法1.計(jì)算速度 很快。兩種方法對(duì)比練習(xí)。顯著性檢驗(yàn)可以證明mmmi=lfc=lt=l這說(shuō)明矩陣X的方差大小可以簡(jiǎn)單的用特征根的大小來(lái)表示。入越高說(shuō)明其

6、對(duì)應(yīng)的 模態(tài)越重要，對(duì)總方差的貢獻(xiàn)越大。第A個(gè)模態(tài)對(duì)總的方差解釋率為HI1=1x 100%44即使是隨機(jī)數(shù)或者虛假數(shù)據(jù)，放在一起進(jìn)行EOF分析，也可以將其分解成一 系列的空間特征向量和主成分。因此，實(shí)際資料分析中得到的空間模態(tài)是否是隨 機(jī)的，需要進(jìn)行統(tǒng)計(jì)檢驗(yàn).North等(1982)的研究指出，在95%置信度水平下的特 征根的誤差入是特征根，2是數(shù)據(jù)的有效自由度，這在前面相關(guān)系數(shù)分析中已經(jīng)有介紹(見(jiàn)4 頁(yè)相關(guān)內(nèi)容)。將入按順序依次檢查，標(biāo)上誤差范圍。如果前后兩個(gè)入之間誤差范圍 有重疊，那么他們之間沒(méi)有顯著差別。Q 圖A.16是對(duì)1949 - 2002年北半球1月平均海平面氣壓，做距平處理處理及

7、面積 加權(quán)后進(jìn)行EOF分析的結(jié)果。從特征根誤差范圍看，第一和第二模態(tài)存在顯著差 別，第二和第三模態(tài)之間也存在顯著差別。但是第三特征根和第四及以后的特征 根之間沒(méi)有顯著的差別。如果要分析主要的模態(tài)的話，最好只選擇前三個(gè)進(jìn)行分 析。練習(xí)：利用E,V =eig(C)計(jì)算矩陣X的特征向量和主成分X二2 6 1 5 2;9 4 0 5 4；X(l,:)=X(1,:)-nean(X(l,:); X(2,:)=X(2,:)-mean(X(2,:);得到X的距平值：X-1.202.80-2.201.80-1.204.60-0.40-4.400.60-0.40WJ co-variance matrixC-X&#

8、187;X>/5;協(xié)方差陣O3.760.920.928.24E0F,E-eig(C); X V: eigenvectors; E: eigenvaluesPC-EOFX;X% reverse the orderE=fliplr(flipud(E)lambd&"diag(E); X retain eigenvalues onlyEOF=fliplr(EOF)POflipud(PC)得到E0F=0.19-0.980.980.19得列轉(zhuǎn)征根8.42003.58得到主成分PO4.280.152.07282-4.741.310.941.65-0.621.1045#處checkEO

9、FEOF* % = I 栓查EOF的正交性得劉：1.00 001.00POPL/5 7. = lambda 檢查PC的正交ttflFfl:8.420.000.003.58E0FPC % X可以完全恢復(fù)X的距半值：1.80 -1.200.60-0.40-1.202.80-2.204.60-0.40-4.40練習(xí)：利用u,S, V =svd(X)計(jì)算矩陣X的特征向量和主成分X2 6 1 52；9 4 0 54；X(lf：)-X(l,:)-mean(X(l,:)>X(2,:)«X(2,:)-mean(X(2>:)9X的距平是：-1.202.80-2.201.801 204.60

10、0.40-4.400.60-0.40nJ>S>V-svd(X);得到u=0.190.980.980.19S=6.490 0 000 4.23 0 00V=0.66-0.490.560.09-0.060.020.670.63-0.320.22-0.73-0.310.530.25-0.160.140.390.030.910.06010-0.26-0.020.060.96EOF-U；PC=SV» 得到PO4.280.15-4.740.94062-2.072.82-1.311.65-1.10E-S.*2/5;lambdaE的軟值與上面得到的杵征根完全一樣即E乞8.42000003

11、.58000E0DPC %可以完全恢復(fù)X的距平值：1202.80-2.201.80-1.204.60-0.40-4.400.60-0.4046810Number1970198019902000YearE0F1 26.1%0.040.020-0.02-0.04-0.06-0.08圖A.1G:北半球1月海平面氣壓EOF分析的第一特征向量.為特征根及95%信度 誤差，（b）第一特征向量，（c）第一主成分，（1）第一主成分偏強(qiáng)+”時(shí)海平面氣壓的 變化量（hPa）. 1949 - 2002, NCEP/NCAR再分析資料結(jié)果展示通常情況下，主成分是有單位的，即反映的是矩陣X的單位，而空間特征向量是 無(wú)量

12、綱的。不過(guò)實(shí)際應(yīng)用中常常對(duì)EOF分析得到的主成分和特征向量進(jìn)行標(biāo)準(zhǔn)化 處理得到新的卩C*和EOF*PCk) =EOF® = EOF(k) y/X或者是簡(jiǎn)單地將PC標(biāo)準(zhǔn)化，使得其平均值為()，標(biāo)準(zhǔn)差為1。再將它與原始資 料矩陣X進(jìn)行回歸分析，這樣就得到PC變化一個(gè)單位時(shí)，變量X對(duì)應(yīng)的響應(yīng)的 空間特征及其強(qiáng)度.這樣得到的回歸系數(shù)的空間分布與空間特征向量的分布特 征空間分布特征是相似的，但是回歸系數(shù)可以看出相應(yīng)的變化的數(shù)量大小.如 圖 A.16(d)。空間模態(tài)應(yīng)該與主成分配合進(jìn)行分析。二者符號(hào)是相對(duì)應(yīng)的。分析中保留的模態(tài)的數(shù)目，沒(méi)有嚴(yán)格規(guī)定，還取決于分析目的。一般取滿 足North準(zhǔn)則；或

13、者有明確物理意義。數(shù)據(jù)性質(zhì)與預(yù)處理誤差(2) 資料的處理。原始場(chǎng)，距平場(chǎng)，與標(biāo)準(zhǔn)化場(chǎng)例子：我國(guó)160站夏季降水量的EOF分析(圖A.17)(3) 空間樣本點(diǎn)。大范圍的空間數(shù)據(jù)，特別需要注意資料空間牛表性。非均勻場(chǎng)與均勻分布場(chǎng)；空間抽樣；面積加權(quán)。北半球1月SLP例子時(shí)空轉(zhuǎn)換有時(shí)空間樣本Z72遠(yuǎn)大于時(shí)間序列長(zhǎng)度計(jì)算 2 X TH矩陣的特征根很困難，可以 考慮對(duì)其進(jìn)行時(shí)空轉(zhuǎn)換。矩陣A = XXf和3 = %用的特征根不同，但是特征向(a) EOF1 90.7%(c) EOF1 13.4%R0.160.140.120.10.08J0.060.040.020.150.10.050-0.05-0.1(b

14、) EOF1 15.8%0.150.10.050-0.05-0.1-0.15-0.2(d) EOF1 11.3%圖A.17:我國(guó)東部地區(qū)夏季降水量EOF分析第一特征向量。（町原始值、（b）距平 值，（c）距平百分率，（d）標(biāo)準(zhǔn)化值.1951 - 2002資料.量是一樣的。而可以證明c = x X，和c* = xr X有相同的特征根.但特征向量不 同.因此，通過(guò)時(shí)空轉(zhuǎn)換可以求X'X矩陣的特征根，進(jìn)而計(jì)算XX矩陣的特征向 量.即有c* X V* = V" X A"是C?的特征向量，A是特征根對(duì)角矩陣。根據(jù)廣是可以求出C的特征向量的, 首先計(jì)算匕=X X V*；對(duì)乞進(jìn)行處

15、理得到C的前77個(gè)特征向量以50#得到特征向量U后，就可以計(jì)算相應(yīng)的主成分PC = VT X X前面計(jì)算得到的EOF維數(shù)是m x m,而通過(guò)時(shí)空轉(zhuǎn)換得到的EOF維數(shù)只 有m x n.即只能得到前幾個(gè)特征向量。不過(guò)實(shí)際應(yīng)用中對(duì)結(jié)杲影響并不大，因?yàn)?通常我們只關(guān)心前幾個(gè)最重要的模態(tài)。下面是一個(gè)簡(jiǎn)單例子，有一個(gè)矩陣X,維數(shù)是5x2,先直接計(jì)算矩陣XXV勺5個(gè) 特征向量，然后再利用時(shí)空轉(zhuǎn)換方法計(jì)算其前2個(gè)特征向量。X= 一204.602.80 -0.40-2.20 -4.401.80 0.60-1.20 -0.40Vl>El=eig(X*X9； 7.7.Vl=fliplr(Vl);%El=fli

16、plr(flipud(El);% 得到特征向*V1=-0.660.49-0.45-0.02-0.67-0.140.730.31-0.56-0.15-0.15-0.15-0.32-0.72-0.17-0.14-0.39-0.42-0.580.570.100.260.53-0.77-0.19得到特征根El=42.110 000017.89 00000 00000 00000 000如果進(jìn)行時(shí)空轉(zhuǎn)換的話，計(jì)算結(jié)果是：V2,E2=eig(X,*X);%V2=fliplr(V2);%E2=fliplr(flipud(E2)得到特征向量V2=0.19-0.980.980.19得到特征根E2二42.1100

17、17.89可見(jiàn)Ei和E2是一樣的。再計(jì)算XX，矩陣的第一特征向量:Va=X*V2; %V_kl=Va(:,l)/sqrt(E2(l,l);得到：0.660.02-0.730.14-0.10計(jì)算XX，矩陣的第二特征向量是：V_k2=Va(:,2)/sqrt(E2(2,2);得到：0.4952-0.670.31-0.390.26可見(jiàn)，用時(shí)空轉(zhuǎn)換方法得到的2個(gè)特征向量，與前面直接計(jì)算矩陣XX和矩 陣£匕匕的得到的前兩個(gè)特征向量完全一致。當(dāng)數(shù)據(jù)量很大時(shí)，如對(duì)全球1月 份2.5。分辨率的再分析l()()0hPa高度場(chǎng)(G)進(jìn)行EOF分析時(shí)，空間點(diǎn)的數(shù)量是772 = 10512,時(shí)間長(zhǎng)度71 =

18、 54,則矩陣C =的維數(shù)是10512 x 10512,而如果用時(shí)空變換方法，則矩陣C* = 的維數(shù)是54 x 54,很快就可以計(jì)算出前54個(gè)特征向 量.高分辨率的遙感數(shù)據(jù)如NDVI,其空間維數(shù)遠(yuǎn)比氣象數(shù)據(jù)大，但其長(zhǎng)度通常只 有20年左右，因此進(jìn)行EOF分析時(shí)常需要借助吋空變換手段。A.7.1 REOF 分析算法：程序varimax.mo模態(tài)數(shù)目的選擇。A.8 SVD分析EOF分析中一次只分析了一個(gè)變量X,地學(xué)中常常涉及多個(gè)要素場(chǎng)之間的關(guān) 系。分析多個(gè)要素場(chǎng)關(guān)系的方法也有很多，包括混合EOF(coinbined empirical orthogonal function,縮寫(xiě)CEOF),奇異值

19、分解(singular value decomposition,縮 寫(xiě)S'D)分析，典型相關(guān)(canonical correlation analysis,縮寫(xiě)CCA)等。他們本質(zhì) 上是相同的的。這里主要介紹SVD分析。需要指出的是這里SVD分析只的是利 用SVD方法檢測(cè)兩個(gè)要素場(chǎng)相關(guān)模態(tài)和分析的過(guò)程，SVD本身只是對(duì)矩陣運(yùn)算求 其奇異值及廣義逆等，因此不要將二者混淆。算法兩個(gè)矩陣X和Y,維數(shù)分別是m x ri和p x "先計(jì)算他們的協(xié)方差陣C = £XYT, C的維數(shù)是mxp進(jìn)行SVD分解得到C = l/工廠,/是對(duì)應(yīng)X的空間模態(tài)，U是對(duì)應(yīng)Y的空間模態(tài)，刀對(duì)角線為

20、奇異值了Mat lab 中命令為U,S, V=svd(C)主成分。X的主成分是.4 = UTX, Y的主成分是B = VTY解釋率。丁與X和Y的協(xié)方差平方成正比，因此解釋率是牛2- X 100%結(jié)果解釋：實(shí)例分析1982-20()()年春李北半球NDVI和氣溫之間的耦合關(guān)系.首先將每一個(gè)格點(diǎn)上 的NDVI和溫度都處理成對(duì)1982-2000年的距平，再相乘得到協(xié)方差陣，對(duì)協(xié)方差 陣進(jìn)行SVD分析，可以得到奇異值，每一個(gè)奇異值對(duì)應(yīng)的NDVI和溫度的模態(tài)， 以及每一種模態(tài)的時(shí)間系數(shù).結(jié)果見(jiàn)圖A.18。春季植被NDVI對(duì)溫度的響應(yīng)信號(hào)非常強(qiáng)。二者之間的協(xié)方差高度集中在最前 面的幾對(duì)模態(tài)中。第一到第七對(duì)

21、模態(tài)，解釋率分別為42.6, 19.5, 10.3, 7.7, 5.0, 4.2和2.3%, 這7對(duì)模態(tài)的總解釋率高達(dá)91.6%,說(shuō)明整體上來(lái)看春季NDVI與溫度的關(guān)系是很 密切的。二者之間最主要也是最重要的耦合關(guān)系已經(jīng)包含在這前面幾個(gè)模態(tài)之中 了，這也表明我們只分析這幾個(gè)模態(tài)就已經(jīng)足夠了。其中最重要的第一對(duì)模態(tài)中心在西西伯利亞。第二對(duì)模態(tài)的主要特征是整個(gè) 北美大陸表現(xiàn)為相同符號(hào)的變化，中心在美國(guó)的東北部地區(qū)。這前兩對(duì)模態(tài)的空 間尺度都很大，屬于大陸尺度.第三及以后的各對(duì)模態(tài)尺度相對(duì)較小，都是區(qū) 域性的.而且這些模態(tài)表現(xiàn)出NDVI與溫度異常的高度一致性，正的溫度中心對(duì) 應(yīng)NDVI的正中心，負(fù)的溫度中心對(duì)應(yīng)NDV1的負(fù)中心。通常最強(qiáng)的NDVI