四、以Matlab來做「頻譜分析」

以「台灣地區河川流量站資料庫」中綠水觀測站1956至2004年的資料來看，對每月平均流量的記錄對頻率加以分析（扣除1956、1959、1990、2000這四年的資料，因為當中有些月份的資料有欠缺），希望可以得到某些訊息。

----------------------------------------------------------------------------------------------------
close all
clear all
f=load('1956to2004month.txt')
x=f(1:end, 1)
y=f(1:end, 2)

figure(1); plot(x,y)
xlabel('year'); % X軸註解，以1957年為0，每隔一個月，增加1/12，1958年為2，依此類推
ylabel('流量 cms');% Y軸註解
title('1956～2004年綠水觀測站月平均流量圖');
[Pxx,f]=periodogram(y,[],1024,1);
magnitude=abs(Pxx);
figure(2); plot(f,magnitude), grid


----------------------------------------------------------------------------------------------------

圖一為每月的流量分佈

可以在頻譜分析圖中，看到兩個峰值（peak），第一個峰值出現在頻率＝0處，可以忽略不管。第二個峰值出現在頻率約為0.08之處，換算成週期為12.5年。
意即每月的平均流量多寡，約每隔12.5年，會有一次循環。

----------------------------------------------------------------------------------------------------
曾於課堂上聽老師介紹過－自然界的諸多訊號可以用sin函數來檢視，
因此，底下嘗試用sin函數來fit上述的資料。

假設 y= a*sin(b*x) 為時間與流量的關係函數
則以matlab來算迴歸函數
----------------------------------------------------------------------------------------------------

clear all ; close all ;
y0 = load('1956to2004month.txt');
xdata = y0(1:end, 1);
ydata = y0(1:end, 2);
f = fittype('a*sin(b*x)');
fit1 = fit(xdata,ydata,f,'StartPoint',[1 1]);
fdata = feval(fit1,xdata);
I = abs(fdata - ydata) > 1.5*std(ydata);
outliers = excludedata(xdata,ydata,'indices',I);
figure(1); plot(fit1,'r-',xdata,ydata,'k.',outliers,'m*') %數值超過1.5個標準差以上，以粉紅色的*標示
text(4.5,-30,'Y = 5.244 × Sin ( 0.8625X )') %標記所算出的回歸方程式
----------------------------------------------------------------------------------------------------

所得之Sin函數為 Y = 5.244 × Sin ( 0.8625X )
算出週期約為7.6年( 2*pi/0.8625≒7.6)

意即每月的平均流量多寡，約每隔7.6年，會有一次循環。

推測與上一個頻譜分析所做出的結果不一樣的原因：（若有錯，歡迎指正）
1.兩個不一樣的方法，在資料處理上本來就有誤差。
2.用Sin函數去做迴歸時，信心程度只有95％，Coefficients (with 95% confidence bounds)
3.資料數仍不夠多，所找出的規則性不夠明確。
4.筆者本身的數學有待加強................




如果以不同的迴歸方式來計算迴歸線時，以下為例：
-------------------------------------------------------------------------------------------------
clear all ; close all ;
y0 = load('1956to2004month.txt');
xdata = y0(1:end, 1);
ydata = y0(1:end, 2);
f = fittype('a*sin(b*x)');
fit1 = fit(xdata,ydata,f,'StartPoint',[1 1]);
fdata = feval(fit1,xdata);
I = abs(fdata - ydata) > 1.5*std(ydata);
outliers = excludedata(xdata,ydata,'indices',I);
figure(1); plot(fit1,'r-',xdata,ydata,'k.',outliers,'m*') %數值超過1.5個標準差以上，以粉紅色的*標示
text(4.5,225,'紅色線 ： Y = 5.244 × Sin ( 0.8625 X )') %標記所算出的回歸方程式
hold on ;
fit2 = fit(xdata,ydata,f,'StartPoint',[1 1],'Exclude',outliers);
fit3 = fit(xdata,ydata,f,'StartPoint',[1 1],'Robust','on');
plot(fit2,'g-')
plot(fit3,'b:')
text(4.5,200,'綠色線 ： Y = 1.879 × Sin ( 1.059 X )') %標記所算出的回歸方程式
text(4.5,175,'藍色線 ： Y = 4.459 × Sin ( 0.7496 X )') %標記所算出的回歸方程式

--------------------------------------------------------------------------------------------------
所得到的圖如下：
紅色線為之前第一種所用的迴歸方式；
綠色線為Exclude Dada所算得的迴歸線；
藍色線為以Robust迴歸（強韌迴歸）方式來計算迴歸線，所算得的週期（約為8.4年）略大於第一種方式（7.6年）。

三、以「迴歸擬合線」來推算：立霧溪的輸砂量

1.標定曲線是描述 懸浮泥沙沉積量(L)和河水流量(Q)之間的關係，（Campbell和Bauder 1940; Crawford 1991)。這種關係通常以一個冪函數來表達：

沈積量＝k‧（流量）^b

2.立霧溪的輸砂量資料，取自於綠水觀測站。但輸砂量的紀錄卻不是一個連續的紀錄。
3.希望能透過matlab這個工具，希望已知的資料，來推算一個較長時間的總輸砂量和平均輸砂量。

---------------------------------------------------------------------------------------------------
以1986年的資料來計算（取自於水資源局1986年水文年報），其中某些觀測資料裡，只有水流量，而輸砂量的資料顯示為『0』，我將這些資料剔除，我的理由有二：
一、「0」無法取對數值。
二、這些資料會大幅干擾「迴歸擬合線」，R^2值也會大幅降低。

以下資料是已經剔除輸砂量為「0」的情形
lq=load('1986.txt') %取有輸砂量資料的欄位
Lq=log(lq(1:end ,4)) %流量取對數，相當於變數X
Ll=log(lq(1:end ,6)) %輸砂量取對數，相當於變數Y
figure(1); plot(Lq,Ll,'or') %畫出對數xy的點圖，以紅色圈圈表示
hold on %留住第一張圖，並把下面的也畫到第一張圖上
p=polyfit(Lq,Ll, 1) %以取完對數的資料點，算迴歸線的斜率p(1)與截距p(2)，乘冪為1次
plot(Lq, p(1)*Lq+p(2)) %以算出的斜率與截距，代入Y=mX+b ，畫迴歸線
title('1986年綠水站流量與輸砂量關係圖');%圖片標題
xlabel('Log Q cms'); % X軸註解
ylabel('Log L mt/d'); % Y軸註解
text(4.5,6,'R^2=0.8980'); %文字標註，前面兩位數為座標位置，另R2的值需要手動輸入
text(4.5,7,'Log(L)= 2.4843 × Log(Q) - 1.5303') %標記所算出的回歸方程式
text(4.5,5,'n=30'); %註記樣本數
sse = sum ((polyval(p, Lq) - Ll).^2) %殘差的平方和 sse=Σ（實際Y-模型預測Y）^2
ssr = sum ((polyval(p, Lq) - mean(Ll)).^2) % ssr=Σ（模型預測Y-平均值Y）^2
sst = sum ((Ll - mean(Ll)).^2) % 總變異量 sst=Σ（實際Y-平均值Y）^2
r_square = ssr/sst %決定係數R^2，也可以用1-sse/sst 來算


---------------------------------------------------------------------------------------------------
以下用同一年的所有觀測資料，不做剔除，但為了讓matlab能夠做對數運算，所以將原本輸砂量為0的地方，皆改成10。比較前後兩種方法所計算出的結果，有何差距？



lq=load('1986all.txt') %取有輸砂量資料的欄位
Lq=log(lq(1:end ,1)) %流量取對數，相當於變數X
Ll=log(lq(1:end ,2)) %輸砂量取對數，相當於變數Y
figure(1); plot(Lq,Ll,'or') %畫出對數xy的點圖，以紅色圈圈表示
hold on %留住第一張圖，並把下面的也畫到第一張圖上
p=polyfit(Lq,Ll, 1) %以取完對數的資料點，算迴歸線的斜率p(1)與截距p(2)，乘冪為1次
plot(Lq, p(1)*Lq+p(2)) %以算出的斜率與截距，代入Y=mX+b ，畫迴歸線
title('1986年綠水站流量與輸砂量關係圖'); %圖片標題
xlabel('Log Q cms'); % X軸註解
ylabel('Log L mt/d'); % Y軸註解
text(4.5,6,'R^2=0.7721'); %文字標註，前面兩位數為座標位置，另R2的值需要手動輸入
text(4.5,7,'Log(L)= 3.4699 × Log(Q) - 6.7825') %標記所算出的回歸方程式
text(4.5,5,'n=56'); %註記樣本數
sse = sum ((polyval(p, Lq) - Ll).^2) %殘差的平方和 sse=Σ（實際Y-模型預測Y）^2
ssr = sum ((polyval(p, Lq) - mean(Ll)).^2) %ssr=Σ（模型預測Y-平均值Y）^2
sst = sum ((Ll - mean(Ll)).^2) % 總變異量sst=Σ（實際Y-平均值Y）^2
r_square = ssr/sst %決定係數R^2，也可以用1-sse/sst 來算

可以看出R^2的值從原本的0.8980，下降到0.7721，迴歸線的斜率也從原本的2.4843，上升到3.4699。這樣的差異到底算大還是小？以這兩個不同的迴歸方程式，來計算年輸砂量。

算年輸砂量之前，要先有年總排水量，可到「台灣地區河川流量站資料庫」中取得，1986年綠水觀測站所測得的總排水量為14790.83 cms-day

以第一個迴歸方程式，反算年輸砂量，所得到的數值為4952.3百萬公噸。
MATLAB指令如下：（需延伸自上面的指令，才會有p(1)斜率和p(2)截距的值）

Lq=log(14790.83) %log內放入當年的總排水量，由水文年報中的日平均流量中，每月加總可得
Ll=p(1)*Lq+p(2) %將Lq帶入所算得的迴歸方程式，得出Ll
exp(Ll) %將Ll反算出實際的沈積量（算出結果為科學記號，可用excel看出原始的值）

以第二個迴歸方程式，反算年輸砂量，所得到的數值為334180百萬公噸。

前後兩者所算出的輸砂量相差約「67.48倍」，差距甚遠。
也就是說：有沒有剔除「那些輸砂量為0的資料」，在計算迴歸線時，影響頗大。

-------------------------------------------------------------------------------------------------

有了一年的資料，可以推算一年中的各月輸砂量及年總輸砂量，
亦可累積長時間的資料，來算得更佳的迴歸線，以推算出更有「根據的」輸砂量。

以下，僅以1986到1989年間，代表較長時間的綠水觀測站資料（已剔除那些輸砂量為0的資料），來計算迴歸線。如果有時間，應從1986年累積至今的觀測資料一併列入計算。
-------------------------------------------------------------------------------------------------
lq=load('86to89.txt') %取有輸砂量資料的欄位
Lq=log(lq(1:end ,4)) %流量取對數，相當於變數X
Ll=log(lq(1:end ,6)) %輸砂量取對數，相當於變數Y
figure(1); plot(Lq,Ll,'or') %畫出對數xy的點圖，以紅色圈圈表示
hold on %留住第一張圖，並把下面的也畫到第一張圖上
p=polyfit(Lq,Ll, 1) %以取完對數的資料點，算迴歸線的斜率p(1)與截距p(2)，乘冪為1次
plot(Lq, p(1)*Lq+p(2)) %以算出的斜率與截距，代入Y=mX+b ，畫迴歸線
title('1986年到1989年綠水站流量與輸砂量關係圖'); %圖片標題
xlabel('Log Q cms'); % X軸註解
ylabel('Log L mt/d'); % Y軸註解
text(4.5,6.5,'R^2=0.8416'); %文字標註，前面兩位數為座標位置，另R2的值需要手動輸入
text(4.5,7.5,'Log(L)= 2.1926 × Log(Q) - 0.2538') %標記所算出的回歸方程式
text(4.5,5.5,'n=69'); %註記樣本數
sse = sum ((polyval(p, Lq) - Ll).^2) %殘差的平方和 sse=Σ（實際Y-模型預測Y）^2
ssr = sum ((polyval(p, Lq) - mean(Ll)).^2) %ssr=Σ（模型預測Y-平均值Y）^2
sst = sum ((Ll - mean(Ll)).^2) %總變異量sst=Σ（實際Y-平均值Y）^2
r_square = ssr/sst %決定係數R^2，也可以用1-sse/sst 來算

-------------------------------------------------------------------------------------------------

二、圖示「綠水觀測站」相關的基本資料

1.以GMT畫出綠水附近的地圖。（在此特別感謝嘉鈴學姊提供立霧溪、和平溪與美崙溪的位置資料點）
------------------------------------------------------------------------------------------------
psbasemap -JM6i -X1.5 -R121/122/24/25 -B0.2NWes:."LU-SHUI": -V -K > LU-SHUI.ps

pscoast -JM6i -R121/122/24/25 -Df -W2 -Na -Ia -G162/205/90 -T123/34/1 -L123/32.5/27/100k -V -K -O >> LU-SHUI.ps

psxy hp.txt -JM6i -R121/122/24/25 -MX -W -V -H -K -O >> LU-SHUI.ps
psxy lw.txt -JM6i -R121/122/24/25 -MX -W -V -H -K -O >>LU-SHUI.ps
psxy mk.txt -JM6i -R121/122/24/25 -MX -W -V -H -K -O >> LU-SHUI.ps

psxy temp1.txt -JM6i -R121/122/24/25 -Sa0.25 -G255/0/0 -V -H -K -O >> LU-SHUI.ps
------------------------------------------------------------------------------------------------

補充說明：紅色星號處為綠水觀測站的位置，以下就不再贅述了。

2.畫出綠水觀測站附近的等高線圖
------------------------------------------------------------------------------------------------
makecpt -Crelief -T300/3000/100 -V -Z > high.cpt

psbasemap -JM6i -X1.5 -R121.4/121.6/24.1/24.3 -B0.2NWes:."LU-SHUI": -V -K > contour_lu-shui.ps

pscoast -JM6i -R121.4/121.6/24.1/24.3 -Df -W2 -Na -Ia -T123/34/1 -L123/32.5/27/100k -V -K -O >> contour_lu-shui.ps

grdcontour twdtm_small.grd -JM6i -R121.4/121.6/24.1/24.3 -C500 -A500 -V -K -O >> contour_lu-shui.ps

psxy lw.txt -JM6i -R121.4/121.6/24.1/24.3 -Sc0.05 -W0/200/255 -MX -W -V -H -K -O >>contour_lu-shui.ps

psxy temp1.txt -JM6i -R121.4/121.6/24.1/24.3 -Sa0.4 -G255/0/0 -V -H -K -O >> contour_lu-shui.ps
------------------------------------------------------------------------------------------------

補充說明：藍色線為立霧溪，紅色星號為綠水觀測站。

3.畫出綠水觀測站附近的立體圖
------------------------------------------------------------------------------------------------
makecpt -Crainbow -T0/400/40 -V -Z > depth.cpt

grd2cpt twdtm_small.grd -Crelief -R121.4/121.7/24.0/24.3 -V -Z > lushui_test.cpt
grdgradient twdtm_small.grd -Gtwdtm_small_temp_shade.grd -A90 -Ne0.9 -V


psbasemap -JM6i -X1.5 -R121.4/121.7/24.0/24.3 -B0.2NWes:."LU-SHUI": -V -K > shade_LU-SHUI.ps

grdimage twdtm_small.grd -Itwdtm_small_temp_shade.grd -JM6i -R121.4/121.7/24.0/24.3 -Clushui_test.cpt -V -K -O >> shade_LU-SHUI.ps

psxy lw.txt -JM6i -R121.4/121.7/24.0/24.3 -Sc0.05 -W0/200/255 -MX -W -V -H -K -O >>shade_LU-SHUI.ps

psxy temp1.txt -JM6i -R121.4/121.7/24.0/24.3 -Sa0.4 -G255/0/0 -V -H -K -O >> shade_LU-SHUI.ps

pscoast -JM6i -R121.4/121.7/24.0/24.3 -Dh -W30 -Na -Ia -W10/0/0/255 -T123/34/1 -L123/32.5/27/100k -V -K -O >> shade_LU-SHUI.ps

psscale -D17.5c/11.6c/10c/0.5c -Clushui_test.cpt -B500 -V -O >> shade_LU-SHUI.ps

------------------------------------------------------------------------------------------------

以地形高低和立霧溪走向來看，原始立霧溪的資料點位置似乎有點偏差。

經過修正後，（經度增加0.018，緯度增加0.045），得到上圖，看起來比較能符合地形。但最終應該還是得靠現地測量或是有其他可靠來源，才能驗證立霧溪資料正確性。




4.以3D立體圖來呈現綠水觀測站附近的地形高低起伏
------------------------------------------------------------------------------------------------
makecpt -Crainbow -T0/3500/100 -V -Z > high.cpt

grdgradient twdtm_small.grd -Gtwdtm_small_temp_shade.grd -A0 -Ne0.9 -V

psbasemap -P -JM6i -JZ2i -R121.3/121.55/24.0/24.3/0/3000 -B0.1/0.08/1000:."LU-SHUI": -E180/60 -T121.45/24.32/0.5 -V -K > 3D_LU-SHUI.ps

grdview twdtm_small.grd -JM6i -JZ2i -R121.3/121.6/24.0/24.3/0/3000 -Chigh.cpt -Itwdtm_small_temp_shade.grd -Qs -E180/60 -N0 -V -K -O >> 3D_LU-SHUI.ps

psxyz temp1.txt -JM6i -JZ6i -R121.3/121.6/24.0/24.3/0/3000 -Sa0.4 -G255/0/0 -V -H -K -O >> 3D_LU-SHUI.ps

psscale -D3/-0.2/5/0.15h -Chigh.cpt -B500 -V -K -O >> 3D_LU-SHUI.ps

------------------------------------------------------------------------------------------------

可以看出綠水觀測站附近的地形高低起伏落差相當大。

5.檢視綠水觀測站附近相關的地震，圓點的顏色代表深度（單位：公里），圓點大小代表地震規模。可以看出地震的發生的立體空間分布與規模。

------------------------------------------------------------------------------------------------
5.1 3D圖
從（http://tecdc.earth.sinica.edu.tw/TEC/db/）設定好選取位置的經緯度，選csv輸出
貼到txt檔，再用excel轉檔（排好欄位）後儲存為hijmsus.txt，qawk才能使用
經緯度範圍121/122/24/25，時間從1900到2009.02.28


gawk "{print $6, $5, -$7, -$7, 0.0003*$8^4}" hijmsus.txt > eq_lu-shui.txt

makecpt -Crainbow -T-150/0/3 -V -Z > eq.cpt
psbasemap -JM8i -Jz0.06c -R121/122/24/25/-150/0 -E175/15 -B0.2:Longitude:/0.2:Latitude:/30:Depth::."LU-SHUI EARTHQUAKE 3-D PLOT":10 -V -K > eq.ps
pscoast -JM8i -Jz0.06c -R121/122/24/25/-150/0 -G162/205/90 -E175/15 -Df -W -V -K -O >> eq.ps
psxyz eq_lu-shui.txt -JM8i -Jz0.06c -R121/122/24/25/-150/0 -E175/15 -Sc -Ceq.cpt -V -K -O >> eq.ps
psscale -D6c/-0.2/4/0.2h -Ceq.cpt -B50 -V -O >> eq.ps

上圖左下角有五個圓形，是為了比對地震規模之用，從上到下的地震規模分別為4、5、6、7、8。
------------------------------------------------------------------------------------------------
5.2 2D俯視圖
gawk "{print $6, $5, -$7, -$7, 0.0003*$8^4}" hijmsus2.txt > eq_lu-shui.txt

makecpt -Crainbow -T-150/0/3 -V -Z > eq.cpt
psbasemap -JM5i -Jz0.06c -X3 -R121/122/24/25/-150/0 -E180/90 -B0.2:Longitude:/0.2:Latitude:/30:Depth::."LU-SHUI EARTHQUAKE 2-D PLOT":6 -V -K > eq2.ps
pscoast -JM5i -Jz0.06c -R121/122/24/25/-150/0 -G162/205/90 -E180/90 -Df -W -V -K -O >> eq2.ps

psxyz eq_lu-shui.txt -JM5i -Jz0.06c -R121/122/24/25/-150/0 -E180/90 -Sc -Ceq.cpt -V -K -O >> eq2.ps

rem 以下三行標註和平溪、立霧溪、美崙溪的位置
psxy hp.txt -JM5i -R121/122/24/25 -MX -W -V -H -K -O >> eq2.ps
psxy lw.txt -JM5i -R121/122/24/25 -MX -W -V -H -K -O >> eq2.ps
psxy mk.txt -JM5i -R121/122/24/25 -MX -W -V -H -K -O >> eq2.ps

rem 以下標註綠水觀測站的位置
psxy temp1.txt -JM5i -R121/122/24/25 -Sa0.4 -G0/0/0 -V -H -K -O >> eq2.ps

psscale -D15c/2.5/5/0.2 -Ceq.cpt -B50 -V -O >> eq2.ps

從以上兩圖可以看出：該地的地震多發生在海陸交界處，且多為淺層地震。其餘地震發生在中央山脈地區，但相較而言，為數不多。

------------------------------------------------------------------------------------------------




6.希望由「震源機制解」圖形來觀察綠水觀測站附近地區的地震型態。
------------------------------------------------------------------------------------------------
EQ_CMT_Harvard_2006_Taiwan.dat由 裡設定選取區域的經緯度範圍
output格式選GMT psmeca input，貼在檔案裡，不用再轉檔，就直接可以用在psmeca裡
有關震源機制解可以參考裡的解說

psbasemap -JM6i -R121/122/24/25 -X3 -Y0.7 -B0.5:."LU-SHUI FOCAL MECHANISM SOLUTION PLOT": -V -K> focal_mech_LU-SHUI.ps
pscoast -JM6i -R121/122/24/25 -G162/205/90 -Df -W1 -V -K -O >> focal_mech_LU-SHUI.ps

rem -Sd0.7c/7 前面的0.7表示海灘球的大小尺寸，後面的/7表示海灘球註解文字的大小
psmeca eq_CMT_LU-SHUI.txt -JM6i -R121/122/24/25 -Sd1c/9 -H -V -K -O >> focal_mech_LU-SHUI.ps

rem 以下三行標註和平溪、立霧溪、美崙溪的位置
psxy hp.txt -JM6i -R121/122/24/25 -MX -W -V -H -K -O >> focal_mech_LU-SHUI.ps
psxy lw.txt -JM6i -R121/122/24/25 -MX -W -V -H -K -O >> focal_mech_LU-SHUI.ps
psxy mk.txt -JM6i -R121/122/24/25 -MX -W -V -H -K -O >> focal_mech_LU-SHUI.ps

rem 以下標註綠水觀測站的位置
psxy temp1.txt -JM6i -R121/122/24/25 -Sa0.4 -G255/0/0 -V -H -O >> focal_mech_LU-SHUI.ps

從上圖可以看出，此一地區的地震型態相當複雜，包含了正斷層、逆斷層、平移斷層、平移與逆斷層...等等的複合地震型態。

一、明年立霧溪的百年洪水發生率是多少？如果發生百年洪水，水量又是高到多少？

1. 到「台灣地區河川流量站資料庫」選取立霧溪的綠水觀測站，取出數據。以下選取範圍是從1960年至2004年。
2. Peak_va為當年最大日流量（單位為CMS），並以EXCEL對Peak_va做遞減排序，給予Rank值。
3. 計算Recurrence Interval; RI = (n+1)/m，其中n為資料筆數，m為magnitude。以本例來說：n+1＝46，m＝Rank值。
4. 計算洪水之「年發生機率」(Annual Exceedance Probability) Pe = (1/RI)*100
5. 以RI（取對數）對Peak作圖，可得一個有線性趨勢的圖，以EXCEL加入趨勢線，如下圖一
6. 以Peak對Pe作圖，如下圖二。

Year	Peak_va	Rank	RI	Peak	Pe
1982	2593.46	1	46	2593.46	2.173913
2000	2240.58	2	23	2240.58	4.347826
1997	1942.54	3	15.33333	1942.54	6.521739
1998	1521.52	4	11.5	1521.52	8.695652
1968	1450	5	9.2	1450	10.86957
1973	1450	6	7.666667	1450	13.04348
1994	1439.04	7	6.571429	1439.04	15.21739
1975	1325	8	5.75	1325	17.3913
1962	1310	9	5.111111	1310	19.56522
1961	1250.5	10	4.6	1250.5	21.73913
1971	1160	11	4.181818	1160	23.91304
1981	1078.46	12	3.833333	1078.46	26.08696
1967	1000	13	3.538462	1000	28.26087
1965	950	14	3.285714	950	30.43478
1992	941.01	15	3.066667	941.01	32.6087
1990	905.71	16	2.875	905.71	34.78261
2004	804.34	17	2.705882	804.34	36.95652
1987	751.26	18	2.555556	751.26	39.13043
1979	725.25	19	2.421053	725.25	41.30435
1969	723	20	2.3	723	43.47826
2003	670.98	21	2.190476	670.98	45.65217
1986	617.65	22	2.090909	617.65	47.82609
1996	557.26	23	2	557.26	50
1989	543.12	24	1.916667	543.12	52.17391
2001	531.63	25	1.84	531.63	54.34783
1999	462.94	26	1.769231	462.94	56.52174
1980	444.77	27	1.703704	444.77	58.69565
1966	432	28	1.642857	432	60.86957
1974	415	29	1.586207	415	63.04348
1977	369.13	30	1.533333	369.13	65.21739
1985	364.81	31	1.483871	364.81	67.3913
1960	355	32	1.4375	355	69.56522
1991	324.94	33	1.393939	324.94	71.73913
1964	280	34	1.352941	280	73.91304
1963	223	35	1.314286	223	76.08696
1983	196.33	36	1.277778	196.33	78.26087
1988	193.26	37	1.243243	193.26	80.43478
1972	187	38	1.210526	187	82.6087
1984	186.97	39	1.179487	186.97	84.78261
1970	172	40	1.15	172	86.95652
1995	168.17	41	1.121951	168.17	89.13043
1976	168	42	1.095238	168	91.30435
1978	149.04	43	1.069767	149.04	93.47826
1993	79.68	44	1.045455	79.68	95.65217
2002	53.75	45	1.022222	53.75	97.82609

圖一：百年洪水日流量約為3200cms

圖二：明年發生百年洪水的可能機率為2.17%

使用EXCEL、GMT與Matlab來處理「立霧溪綠水觀測站的水文資料」

自1956年起，台電公司在太魯閣-綠水附近，設置一長期觀測「立霧溪」的觀測站。
此一觀測站，提供了許多立霧溪的水文資料。
如何讓這些「水文資料」提供我們更多有關地球科學方面的訊息，便是我們所要做的努力：

一、明年立霧溪的百年洪水發生率是多少？如果發生百年洪水，水量又是高到多少？
二、圖示「綠水觀測站」相關的基本資料
三、以「迴歸擬合線」來推算：立霧溪的輸砂量
四、以Matlab來做「頻譜分析」