matlab經(jīng)驗正交函數(shù)EOF(轉載)

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

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

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

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

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

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

7、面積 加權后進行EOF分析的結果。從特征根誤差范圍看，第一和第二模態(tài)存在顯著差 別，第二和第三模態(tài)之間也存在顯著差別。但是第三特征根和第四及以后的特征 根之間沒有顯著的差別。如果要分析主要的模態(tài)的話，最好只選擇前三個進行分 析。練習：利用E,V =eig(C)計算矩陣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得列轉征根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可以完全恢復X的距半值：1.80 -1.200.60-0.40-1.202.80-2.204.60-0.40-4.40練習：利用u,S, V =svd(X)計算矩陣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 %可以完全恢復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）第一主成分偏強+”時海平面氣壓的 變化量（hPa）. 1949 - 2002, NCEP/NCAR再分析資料結果展示通常情況下，主成分是有單位的，即反映的是矩陣X的單位，而空間特征向量是 無量

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

13、者有明確物理意義。數(shù)據(jù)性質與預處理誤差(2) 資料的處理。原始場，距平場，與標準化場例子：我國160站夏季降水量的EOF分析(圖A.17)(3) 空間樣本點。大范圍的空間數(shù)據(jù)，特別需要注意資料空間牛表性。非均勻場與均勻分布場；空間抽樣；面積加權。北半球1月SLP例子時空轉換有時空間樣本Z72遠大于時間序列長度計算 2 X TH矩陣的特征根很困難，可以 考慮對其進行時空轉換。矩陣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:我國東部地區(qū)夏季降水量EOF分析第一特征向量。（町原始值、（b）距平 值，（c）距平百分率，（d）標準化值.1951 - 2002資料.量是一樣的。而可以證明c = x X，和c* = xr X有相同的特征根.但特征向量不 同.因此，通過時空轉換可以求X'X矩陣的特征根，進而計算XX矩陣的特征向 量.即有c* X V* = V" X A"是C?的特征向量，A是特征根對角矩陣。根據(jù)廣是可以求出C的特征向量的, 首先計算匕=X X V*；對乞進行處

15、理得到C的前77個特征向量以50#得到特征向量U后，就可以計算相應的主成分PC = VT X X前面計算得到的EOF維數(shù)是m x m,而通過時空轉換得到的EOF維數(shù)只 有m x n.即只能得到前幾個特征向量。不過實際應用中對結杲影響并不大，因為 通常我們只關心前幾個最重要的模態(tài)。下面是一個簡單例子，有一個矩陣X,維數(shù)是5x2,先直接計算矩陣XXV勺5個 特征向量，然后再利用時空轉換方法計算其前2個特征向量。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如果進行時空轉換的話，計算結果是：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可見Ei和E2是一樣的。再計算XX，矩陣的第一特征向量:Va=X*V2; %V_kl=Va(:,l)/sqrt(E2(l,l);得到：0.660.02-0.730.14-0.10計算XX，矩陣的第二特征向量是：V_k2=Va(:,2)/sqrt(E2(2,2);得到：0.4952-0.670.31-0.390.26可見，用時空轉換方法得到的2個特征向量，與前面直接計算矩陣XX和矩 陣£匕匕的得到的前兩個特征向量完全一致。當數(shù)據(jù)量很大時，如對全球1月 份2.5。分辨率的再分析l()()0hPa高度場(G)進行EOF分析時，空間點的數(shù)量是772 = 10512,時間長度71 =

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

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

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

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