到目前為止僅單就水流量及含沙量兩筆資料做比較，
未來可以在加入更多的資料做分析，
想像是甚麼因子會影響水文資料的記錄?
流速、颱風、地震都可以再加入做討論。
詳細記錄的水文資料得來不易，
此筆台電綠水觀測站的記錄資料可以再做更多的應用 :
1. 將流速與水流量互相對比，尋找是否有相互關係。
2. 把流速和含沙量拿來做相同的分析，流速對於含沙量的影響
可能更甚於水流量。
3. 颱風會帶來強烈的降雨，將颱風的時間加入資料當中，
分析在颱風過後，水流量是否大幅上升，含沙量的變化又是如何。
4. 強地震可能造成山區崩塌，而崩塌會增加水中的沉積物通量，
將地震的時間加入資料當中，分析在地震過後是否立即反應這個過程，
或是在地震後的一段長時間中，觀察含沙量是否呈現較高的趨勢。
2011年6月27日 星期一
五、水文資料分析(三)
6.1986~2008水流量與含沙量分佈
水流量和含沙量之間的關係無法以線性表示，
那麼兩者之間的分佈是否有相對應關係?
因此我將水流量和含沙量的分佈放在一起比較。
axes(ax(1));ylabel('含沙量ppm');
title('1986~2008水流量與含沙量比較圖');
axes(ax(2));ylabel('水流量cms');
![]()
為了能同時比對水流量和含沙量，
使用plotyy可以同時將兩筆資料畫在同一個x軸上，
水流量為綠色線，含沙量為藍色線。
水流量與含沙量雖無法以線性方程式表現，
但在我們的想像當中，兩者應該還是有關係存在，
水流量大的時後，含沙量的記錄應該也會比較高，
由圖觀察卻發現，兩者不一定呈現正相關，有幾個現象:
1.水流量大時不一定含沙量會跟著增加
水流量(綠)往上升，但含沙量(藍)依然很低。
![]()
2.有時後水流量並無明顯提升，含沙量卻大幅增加
約在1999年前後，水流量(綠)的分佈屬正常起伏，
但含沙量(藍)卻呈現不規則波動，有幾次大幅高於水流量，
剛好1999年發生台灣近年來最大的集集地震，
是不是因為此次地震造成這個現象，也值得深入討論。
![]()
![]()
![]()
水流量和含沙量之間的關係無法以線性表示，
那麼兩者之間的分佈是否有相對應關係?
因此我將水流量和含沙量的分佈放在一起比較。
程式碼
clear all;close all;
s=load ('taipower.txt');
t=julian(s(: ,2),s(: ,4),s(: ,5),s(: ,6))-julian(1986,1,1,0);%julian day
vol=s(:,16);%流量cms
sand=s(:,18);%含砂量ppm
[ax,h1,h2]=plotyy(t/365.25+1986,sand,t/365.25+1986,vol);
%將水流量和含沙量畫在同一個時間橫軸上s=load ('taipower.txt');
t=julian(s(: ,2),s(: ,4),s(: ,5),s(: ,6))-julian(1986,1,1,0);%julian day
vol=s(:,16);%流量cms
sand=s(:,18);%含砂量ppm
[ax,h1,h2]=plotyy(t/365.25+1986,sand,t/365.25+1986,vol);
axes(ax(1));ylabel('含沙量ppm');
title('1986~2008水流量與含沙量比較圖');
axes(ax(2));ylabel('水流量cms');
為了能同時比對水流量和含沙量，
使用plotyy可以同時將兩筆資料畫在同一個x軸上，
水流量為綠色線，含沙量為藍色線。
水流量與含沙量雖無法以線性方程式表現，
但在我們的想像當中，兩者應該還是有關係存在，
水流量大的時後，含沙量的記錄應該也會比較高，
由圖觀察卻發現，兩者不一定呈現正相關，有幾個現象:
1.水流量大時不一定含沙量會跟著增加
水流量(綠)往上升，但含沙量(藍)依然很低。
2.有時後水流量並無明顯提升，含沙量卻大幅增加
約在1999年前後，水流量(綠)的分佈屬正常起伏，
但含沙量(藍)卻呈現不規則波動，有幾次大幅高於水流量，
剛好1999年發生台灣近年來最大的集集地震，
是不是因為此次地震造成這個現象，也值得深入討論。
3.當然還是有兩者的記錄值同時提高的時後
1995年前後及1991年的資料，
水流量(綠)及含沙量(藍)兩者有同步提高的現象。
小結 - 水流量與含沙量的相對應關係比想像中的複雜，
或許水流量的大小不是影響含沙量的主要因子，
與前一段得到相同的小結，我們需要更多資料的加入，
才有辦法解釋含沙量數值高低的意義。
四、水文資料分析(二)
4.1986~2008含沙量分佈圖
程式碼
clear all; close all
f=load('taipower.txt')
x=julian(f(: ,2),f(: ,4),f(: ,5),f(: ,6))-julian(1986,1,1,0);%julian day換算
y=f(1:end, 18)
figure(1); plot(x/365.25+1986,y)%繪製xy分佈圖
xlabel('year'); % X軸註解
ylabel('含沙量 ppm');% Y軸註解
title('1986～2008年立霧溪含沙量分佈圖');%給予標題
a=load('sand2.txt')
b=julian(a(: ,1),a(: ,2),a(: ,3),a(: ,4))-julian(1986,1,1,0);%julian day換算
c=log(a(:, 5))%y軸取log
figure(2); plot(b/365.25+1986,c)
xlabel('year'); % X軸註解
ylabel('Log 含沙量 ppm');% Y軸註解
title('Log 1986～2008年立霧溪含沙量分佈圖');給予標題
![]()
![]()
第一張圖為正常分佈圖，第二張圖將y軸取Log。
先前提到含沙量的資料並非每天一筆，資料的時間間隔也不固定，
因此使用ulian day換算，將每筆資料標記在正確的時間序列上，
2000年後的資料密度遠低於先前，所以圖上顯示的資料較不準確，
2000~2004的資料極少，且多為0。
如分佈圖所示，大部分的資料都落在8000ppm以下，
含沙量最高的峰值出現在1997年，1997年後的含沙量明顯比先前高，
1999年921地震附近的資料也較背景值高。
5.水流量與含砂量的關係
含沙量的資料從1986年才開始記錄，每筆資料除了含沙量之外，
也同時記錄了水流量、水位、水溫、流速，
因此我把水流量和含沙量的資料擷取出來相互比對，給予趨勢線，
由於許多筆含沙量的資料為0，先把0的資料拿掉繪第一條線(紅)，
同時使用robustfit繪製第二條線(綠)(由電腦去除較不可信的資料)。
程式碼
clear all; close all
fs=load('flow&sand-0.txt') %讀取資料,含沙量去除0的資料
Lf=log(fs(1:end ,1)) %流量取log，x
Ls=log(fs(1:end ,2)) %含沙量取log，y
figure(1); plot(Lf,Ls,'o') %畫xy分佈圖
hold on %保留第一張圖
p=polyfit(Lf,Ls,1) %計算迴歸線的斜率p(1)與截距p(2),乘冪次1
plot(Lf, p(1)*Lf+p(2),'r') %代入y=ax+b,畫迴歸線
b=robustfit(Lf,Ls) %用robustfit畫另一條回歸線
figure(1);hold on;plot(Lf,b(2)*Lf+b(1),'g') %套到相同圖
title('1986-2008年立霧溪流量與輸砂量關係圖');%給予標題
xlabel('Log 流量 cms'); % X軸註解
ylabel('Log 含沙量 ppm'); % Y軸註解
sse = sum ((polyval(p, Lf) - Ls).^2)
% Residual sum of squares剩餘平方何 SSE=Σ(y-y預測值)^2
ssr = sum ((polyval(p, Lf) - mean(Ls)).^2)
% Regression sum of squares回歸平方和 SSR=Σ(y預測值-y平均值)^2
sst = sum ((Ls - mean(Ls)).^2)
% Total sum of squares總平方和 SST=Σ(y-y平均值)^2
r_square = ssr/sst %R^2回歸線解釋資料的程度SSR/SST
sse1 = sum ((polyval(b, Lf) - Ls).^2)
% Residual sum of squares剩餘平方何 SSE=Σ(y-y預測值)^2
ssr1 = sum ((polyval(b, Lf) - mean(Ls)).^2)
% Regression sum of squares回歸平方和 SSR=Σ(y預測值-y平均值)^2
sst1 = sum ((Ls - mean(Ls)).^2)
% Total sum of squares總平方和 SST=Σ(y-y平均值)^2
r_square1 = ssr/sst %R^2回歸線解釋資料的程度SSR/SST
![]()
紅色為正常的趨勢線，綠色為robustfit所繪的趨勢線。
polyfit 斜率p1=0.7702 截距p2=-0.7939
sse=282.5970 ssr=176.0673 sst=458.6642
R^2=0.3839
robustfit 斜率b1=0.7836 截距b2=-0.9117
sse1=8.1651 e+004 ssr1=8.0775 e+0.04 sst1=458.6642
R^2=0.3839
polyfit和robustfit所畫的趨勢線相差不多，
整體資料分散，趨勢線解釋資料的程度並不好，
代表水流量與含沙量之間的關係無法用簡單的線性去表示。
若我們不將0的資料刪除，而用1取代，畫出來的趨勢會如何?
程式碼
clear all; close all
fs=load('flow&sand.txt') %讀取資料,含沙量去除0的資料
Lf=log(fs(1:end ,1)) %流量取log，x
Ls=log(fs(1:end ,2)) %含沙量取log，y
figure(1); plot(Lf,Ls,'o') %畫xy分佈圖
hold on %保留第一張圖
p=polyfit(Lf,Ls,1) %計算迴歸線的斜率p(1)與截距p(2),乘冪次1
plot(Lf, p(1)*Lf+p(2),'r') %代入y=ax+b,畫迴歸線
b=robustfit(Lf,Ls) %用robustfit話另一條回歸線
figure(1);hold on;plot(Lf,b(2)*Lf+b(1),'g') %套到相同圖
title('1986-2008年立霧溪流量與輸砂量關係圖');%給予標題
xlabel('Log 流量 cms'); % X軸註解
ylabel('Log 含沙量 ppm'); % Y軸註解
sse = sum ((polyval(p, Lf) - Ls).^2)
% Residual sum of squares剩餘平方何 SSE=Σ(y-y預測值)^2
ssr = sum ((polyval(p, Lf) - mean(Ls)).^2)
% Regression sum of squares回歸平方和 SSR=Σ(y預測值-y平均值)^2
sst = sum ((Ls - mean(Ls)).^2)
% Total sum of squares總平方和 SST=Σ(y-y平均值)^2
r_square = ssr/sst %R^2回歸線解釋資料的程度SSR/SST
sse1 = sum ((polyval(b, Lf) - Ls).^2)
% Residual sum of squares剩餘平方何 SSE=Σ(y-y預測值)^2
ssr1 = sum ((polyval(b, Lf) - mean(Ls)).^2)
% Regression sum of squares回歸平方和 SSR=Σ(y預測值-y平均值)^2
sst1 = sum ((Ls - mean(Ls)).^2)
% Total sum of squares總平方和 SST=Σ(y-y平均值)^2
r_square1 = ssr/sst %R^2回歸線解釋資料的程度SSR/SST
![]()
polyfit 斜率p1=2.7541 截距p2=-25.5587
sse=5.6348 e+003 ssr=6.2041 e+003 sst=1.1839 e+004
R^2=0.5240
robustfit 斜率b1=2.9301 截距b2=-27.5598
sse1=6.9533 e+007 ssr1=6.9397 e+0.07 sst1=1.1839 e+004
R^2=0.5240
polyfit和robustfit所畫的趨勢線依然相差不多，
R^2值雖然提高，但由圖看來，資料分佈更為分散。
小結- 水流量和含沙量的關係頗為複雜，無法用簡單的線性公式代表，
應該還需要考慮其他的因子，例如: 流速? 水溫? 岩性? 坡度?
在公式中加上其他參數，
或許能夠使用公式回推出一個可供參考的水流量or含沙量。
程式碼
clear all; close all
f=load('taipower.txt')
x=julian(f(: ,2),f(: ,4),f(: ,5),f(: ,6))-julian(1986,1,1,0);%julian day換算
y=f(1:end, 18)
figure(1); plot(x/365.25+1986,y)%繪製xy分佈圖
xlabel('year'); % X軸註解
ylabel('含沙量 ppm');% Y軸註解
title('1986～2008年立霧溪含沙量分佈圖');%給予標題
a=load('sand2.txt')
b=julian(a(: ,1),a(: ,2),a(: ,3),a(: ,4))-julian(1986,1,1,0);%julian day換算
c=log(a(:, 5))%y軸取log
figure(2); plot(b/365.25+1986,c)
xlabel('year'); % X軸註解
ylabel('Log 含沙量 ppm');% Y軸註解
title('Log 1986～2008年立霧溪含沙量分佈圖');給予標題
先前提到含沙量的資料並非每天一筆，資料的時間間隔也不固定，
因此使用ulian day換算，將每筆資料標記在正確的時間序列上，
2000年後的資料密度遠低於先前，所以圖上顯示的資料較不準確，
2000~2004的資料極少，且多為0。
如分佈圖所示，大部分的資料都落在8000ppm以下，
含沙量最高的峰值出現在1997年，1997年後的含沙量明顯比先前高，
1999年921地震附近的資料也較背景值高。
5.水流量與含砂量的關係
含沙量的資料從1986年才開始記錄，每筆資料除了含沙量之外，
也同時記錄了水流量、水位、水溫、流速，
因此我把水流量和含沙量的資料擷取出來相互比對，給予趨勢線，
由於許多筆含沙量的資料為0，先把0的資料拿掉繪第一條線(紅)，
同時使用robustfit繪製第二條線(綠)(由電腦去除較不可信的資料)。
程式碼
clear all; close all
fs=load('flow&sand-0.txt') %讀取資料,含沙量去除0的資料
Lf=log(fs(1:end ,1)) %流量取log，x
Ls=log(fs(1:end ,2)) %含沙量取log，y
figure(1); plot(Lf,Ls,'o') %畫xy分佈圖
hold on %保留第一張圖
p=polyfit(Lf,Ls,1) %計算迴歸線的斜率p(1)與截距p(2),乘冪次1
plot(Lf, p(1)*Lf+p(2),'r') %代入y=ax+b,畫迴歸線
b=robustfit(Lf,Ls) %用robustfit畫另一條回歸線
figure(1);hold on;plot(Lf,b(2)*Lf+b(1),'g') %套到相同圖
title('1986-2008年立霧溪流量與輸砂量關係圖');%給予標題
xlabel('Log 流量 cms'); % X軸註解
ylabel('Log 含沙量 ppm'); % Y軸註解
sse = sum ((polyval(p, Lf) - Ls).^2)
% Residual sum of squares剩餘平方何 SSE=Σ(y-y預測值)^2
ssr = sum ((polyval(p, Lf) - mean(Ls)).^2)
% Regression sum of squares回歸平方和 SSR=Σ(y預測值-y平均值)^2
sst = sum ((Ls - mean(Ls)).^2)
% Total sum of squares總平方和 SST=Σ(y-y平均值)^2
r_square = ssr/sst %R^2回歸線解釋資料的程度SSR/SST
sse1 = sum ((polyval(b, Lf) - Ls).^2)
% Residual sum of squares剩餘平方何 SSE=Σ(y-y預測值)^2
ssr1 = sum ((polyval(b, Lf) - mean(Ls)).^2)
% Regression sum of squares回歸平方和 SSR=Σ(y預測值-y平均值)^2
sst1 = sum ((Ls - mean(Ls)).^2)
% Total sum of squares總平方和 SST=Σ(y-y平均值)^2
r_square1 = ssr/sst %R^2回歸線解釋資料的程度SSR/SST
紅色為正常的趨勢線，綠色為robustfit所繪的趨勢線。
polyfit 斜率p1=0.7702 截距p2=-0.7939
sse=282.5970 ssr=176.0673 sst=458.6642
R^2=0.3839
robustfit 斜率b1=0.7836 截距b2=-0.9117
sse1=8.1651 e+004 ssr1=8.0775 e+0.04 sst1=458.6642
R^2=0.3839
polyfit和robustfit所畫的趨勢線相差不多，
整體資料分散，趨勢線解釋資料的程度並不好，
代表水流量與含沙量之間的關係無法用簡單的線性去表示。
若我們不將0的資料刪除，而用1取代，畫出來的趨勢會如何?
程式碼
clear all; close all
fs=load('flow&sand.txt') %讀取資料,含沙量去除0的資料
Lf=log(fs(1:end ,1)) %流量取log，x
Ls=log(fs(1:end ,2)) %含沙量取log，y
figure(1); plot(Lf,Ls,'o') %畫xy分佈圖
hold on %保留第一張圖
p=polyfit(Lf,Ls,1) %計算迴歸線的斜率p(1)與截距p(2),乘冪次1
plot(Lf, p(1)*Lf+p(2),'r') %代入y=ax+b,畫迴歸線
b=robustfit(Lf,Ls) %用robustfit話另一條回歸線
figure(1);hold on;plot(Lf,b(2)*Lf+b(1),'g') %套到相同圖
title('1986-2008年立霧溪流量與輸砂量關係圖');%給予標題
xlabel('Log 流量 cms'); % X軸註解
ylabel('Log 含沙量 ppm'); % Y軸註解
sse = sum ((polyval(p, Lf) - Ls).^2)
% Residual sum of squares剩餘平方何 SSE=Σ(y-y預測值)^2
ssr = sum ((polyval(p, Lf) - mean(Ls)).^2)
% Regression sum of squares回歸平方和 SSR=Σ(y預測值-y平均值)^2
sst = sum ((Ls - mean(Ls)).^2)
% Total sum of squares總平方和 SST=Σ(y-y平均值)^2
r_square = ssr/sst %R^2回歸線解釋資料的程度SSR/SST
sse1 = sum ((polyval(b, Lf) - Ls).^2)
% Residual sum of squares剩餘平方何 SSE=Σ(y-y預測值)^2
ssr1 = sum ((polyval(b, Lf) - mean(Ls)).^2)
% Regression sum of squares回歸平方和 SSR=Σ(y預測值-y平均值)^2
sst1 = sum ((Ls - mean(Ls)).^2)
% Total sum of squares總平方和 SST=Σ(y-y平均值)^2
r_square1 = ssr/sst %R^2回歸線解釋資料的程度SSR/SST
polyfit 斜率p1=2.7541 截距p2=-25.5587
sse=5.6348 e+003 ssr=6.2041 e+003 sst=1.1839 e+004
R^2=0.5240
robustfit 斜率b1=2.9301 截距b2=-27.5598
sse1=6.9533 e+007 ssr1=6.9397 e+0.07 sst1=1.1839 e+004
R^2=0.5240
polyfit和robustfit所畫的趨勢線依然相差不多，
R^2值雖然提高，但由圖看來，資料分佈更為分散。
小結- 水流量和含沙量的關係頗為複雜，無法用簡單的線性公式代表，
應該還需要考慮其他的因子，例如: 流速? 水溫? 岩性? 坡度?
在公式中加上其他參數，
或許能夠使用公式回推出一個可供參考的水流量or含沙量。
三、水文資料分析(一)
1.月均流量Histogarm
程式碼
clear all; close all
bio=dlmread('flow.txt');%讀取資料
flow=bio(1:end,2);figure(1);%選擇資料區間
[n1,x1]=hist(flow);h1=bar(x1,n1);%繪製長條圖
title('Histogram of 1969~2008 立霧溪月均流量');%加標題
xlabel('水流量 cms')%x軸標示
ylabel('數量 month')%y軸標示
mean(flow)%平均
std(flow)%標準差
![]()
資料區間 1959~2008，每月一筆。
mean = 3.2319e+004
std = 3.0804e+004
月均流量資料的初步呈現，
約70%的月均流量小於16570cms。
2.含沙量Histogram
程式碼
clear all; close all
bio=dlmread('sand.txt');%讀取資料
sand=bio(1:end,2);figure(1);%選擇資料區間
[n1,x1]=hist(sand);h1=bar(x1,n1);%繪製長條圖
title('Histogram of 1986~2008 立霧溪含沙量');%顯示標題
xlabel('含沙量 ppm')%x軸標示
ylabel('次數')%y軸標示
mean(sand)%平均
std(sand)%標準差
![]()
mean = 1.4952e+003
std = 4.0156e+003
資料區間1986~2008，記錄無固定間隔，
超過90%的含沙量小於5000ppm，
含沙量資料初步呈現。
3.1959~2008立霧溪月均流量分佈圖
程式碼
clear all; close all
f=load('flow.txt')
x=f(1:end, 1)
y=f(1:end, 2)
z=log(f(:,2))
figure(1); plot(x,y);
xlabel('時間序列1959~2008'); % X軸註解,一月一筆資料,以1959年1月為0、2月為1/12，
1960年1月為1，以此類推
ylabel('流量 cms');% Y軸註解
title('1959～2008年立霧溪月平均流量圖');
figure(2);plot(x,z);
xlabel('時間序列1959~2008');
ylabel('Log 流量 cms');% Y軸註解
title('Log 1959～2008年立霧溪月平均流量圖');
![]()
![]()
第一張圖為月均流量分佈圖，第二張圖為y軸取log。
一個月一筆資料以1959年1月為0、2月為1/12，1960年則為1，以此類推，
水流量的最大值出現在1999年，50年來並沒有特別乾旱的現象，
僅2003年的數值較往年稍小，表示立霧溪水量較為穩定，
水流量的分佈似遵循某個規律的循環，因為台灣乾濕季分明，
降雨明顯集中於7~10月的颱風季節，可能與此種季節降雨相關。
程式碼
clear all; close all
bio=dlmread('flow.txt');%讀取資料
flow=bio(1:end,2);figure(1);%選擇資料區間
[n1,x1]=hist(flow);h1=bar(x1,n1);%繪製長條圖
title('Histogram of 1969~2008 立霧溪月均流量');%加標題
xlabel('水流量 cms')%x軸標示
ylabel('數量 month')%y軸標示
mean(flow)%平均
std(flow)%標準差
資料區間 1959~2008，每月一筆。
mean = 3.2319e+004
std = 3.0804e+004
月均流量資料的初步呈現，
約70%的月均流量小於16570cms。
2.含沙量Histogram
程式碼
clear all; close all
bio=dlmread('sand.txt');%讀取資料
sand=bio(1:end,2);figure(1);%選擇資料區間
[n1,x1]=hist(sand);h1=bar(x1,n1);%繪製長條圖
title('Histogram of 1986~2008 立霧溪含沙量');%顯示標題
xlabel('含沙量 ppm')%x軸標示
ylabel('次數')%y軸標示
mean(sand)%平均
std(sand)%標準差
mean = 1.4952e+003
std = 4.0156e+003
資料區間1986~2008，記錄無固定間隔，
超過90%的含沙量小於5000ppm，
含沙量資料初步呈現。
3.1959~2008立霧溪月均流量分佈圖
程式碼
clear all; close all
f=load('flow.txt')
x=f(1:end, 1)
y=f(1:end, 2)
z=log(f(:,2))
figure(1); plot(x,y);
xlabel('時間序列1959~2008'); % X軸註解,一月一筆資料,以1959年1月為0、2月為1/12，
1960年1月為1，以此類推
ylabel('流量 cms');% Y軸註解
title('1959～2008年立霧溪月平均流量圖');
figure(2);plot(x,z);
xlabel('時間序列1959~2008');
ylabel('Log 流量 cms');% Y軸註解
title('Log 1959～2008年立霧溪月平均流量圖');
第一張圖為月均流量分佈圖，第二張圖為y軸取log。
一個月一筆資料以1959年1月為0、2月為1/12，1960年則為1，以此類推，
水流量的最大值出現在1999年，50年來並沒有特別乾旱的現象，
僅2003年的數值較往年稍小，表示立霧溪水量較為穩定，
水流量的分佈似遵循某個規律的循環，因為台灣乾濕季分明，
降雨明顯集中於7~10月的颱風季節，可能與此種季節降雨相關。
二、研究區域簡圖(GMT)
台灣位於板塊交界，因造山運動形成中央山脈、雪山山脈等高山地區，
中央山脈東翼太魯閣地區因水樣豐沛的立霧溪及其支流持續不斷切割岩層，
加上板塊擠壓形成的抬升作用，於此處形成獨特的峽谷景觀。
立霧溪為太魯閣國家公園主要河川，長55公里，流域面積616平方公里，
利用台電綠水觀測站所記錄的立霧溪水文資料，
配合matlab和gmt工具，進行資料的分析工作。
首先以GMT繪製研究區域相關簡圖。
1.台灣島附近海底地形、斷層位置及地震分佈圖
![]()
使用地形為ETOPO1，繪製範圍117/125/19/27，
將海岸山脈斷層、琉球海溝、馬尼拉海溝、呂宋島弧、
梨山斷層、西部麓山戴變形前緣以斷層線標記，
黃點為1986年~2008年在範圍內，規模超過6級的地震震央，
(因為水文資料含沙量的資料區間為1986~2008，故地震資料也取相同區間)
地震的分佈多靠近台灣東部或外海，海岸山脈至和平之間是最頻繁的區域。
由海底地形可得知，台灣西部台灣海峽是一個平坦的大陸棚，
東部海岸的高度落差就大很多，屬於斷層海岸。
程式碼
gmtset PLOT_DEGREE_FORMAT ddd:mm:ssF
makecpt -Crainbow -T-7000/4000/200 -Z > 1.cpt
grdgradient ETOPO1_taiwan.grd -A310 -GETOPO1.int -Ne1.0 -M -V
grdimage ETOPO1_taiwan.grd -R117/125/19/27 -JM6i -C1.cpt -IETOPO1.int
-Bf2a2/f2a2WSne:."1986-2008 earthquakes 6<Mw": -V -P -K -Y3i > 1.ps
pscoast -R117/125/19/27 -J -Df -A50 -O -P -V -W0.7p -K >> 1.ps
grdcontour ETOPO1_taiwan.grd -C200 -J -L-15000/0 -R -W0.5p -O -K -V >> 1.ps
psxy 1.txt -R -JM -Gyellow -Sc0.25c -W0.1p -O -K >> 1.ps
pstext taroko.txt -R117/125/19/27 -J -Gred -V -O -K >> 1.ps
psxy ryukyu_trench.asc -R -J -Sf0.5i/0.1ilt -Gblack -Wthick -O -V -K>> 1.ps
psxy south_eastern_fault.asc -R -J -Sf0.5i/0.1ilt -Gblack -Wthick -K -O -V >> 1.ps
psxy deformation_front.asc -R -J -Sf0.5i/0.1irt -Gblack -Wthick -K -O -V >> 1.ps
psxy malina_trench.asc -R -J -Sf0.5i/0.1ilt -Gblack -Wthick -K -O -V >> 1.ps
psxy Lishan_laonung_fault.asc -R -J -Sf0.5i/0.1irt:0.5i -Gblack -Wthick -K -O -V >> 1.ps
psxy coastal_range_fault.asc -R -J -Sf0.5i/0.1irt -Gblack -Wthick -K -O -V >> 1.ps
psscale -D3i/2i/4.5i/0.25ih -O -C1.cpt -I -B2000:BATHYMETRY:/:m: -Y-3i >> 1.ps
2.立霧溪外海海底地形
立霧溪為水量豐沛之溪流，水流亦攜帶大量沉積物，
河流流出山區後形成一個標準形狀的沖積扇，但面積不大，
代表河流所攜帶的沉積物被帶往其他地方沉積，
因此繪製立霧溪出海口外海的地形剖面，
觀察此處是否有海底沖積扇的地形產生。
![]()
上圖為海地地形剖面，下圖為剖面位置，
由剖面可見至10公里外海的海岸地形皆為連續的陡峭下坡，
10公里處一個隆起後繼續更為陡峭的下坡，
發現此處並無明顯的海底沖積扇存在。
那麼立霧溪所攜帶的大量沉積物被攜帶至何處?
是因為海底地形原本就非常陡峭，所以無法形成堆積地形，
或者被沿岸流、洋流帶至其他地方堆積，值得繼續探討。
程式碼
makecpt -Csealand -T-10000/4000/500 -V -Z > 0620.cpt
project -C121.65/24.14 -E123.48/23.52 -G0.5 -N -Q > temp1.d
grdtrack temp1.d -GTaiDBv624.grd -V > temp3.d
grdgradient TaiDBv624.grd -TaiDBv624_shade.grd -A0/30 -Ne0.5 -V
grdimage -JM6i -R119/125/21/27 TaiDBv624.grd -TaiDBv624_shade.grd
-C0620.cpt -B2nwSE -P -V -K > 0620.ps
cat <<END> profile.dat
121.65 24.14
123.48 23.52
END
psxy -JM6i -R119/125/21/27 profile.dat -W10/250/0/0 -V -K -O >> 0620.ps
pstext taroko.txt -R119/125/21/27 -J -G -V -O -K >> 0620.ps
pscoast -JM6i -R119/125/21/27 -Df -W5 -V -K -O >> 0620.ps
grdcontour TaiDBv624.grd -C500 -J -A500 -L-15000/0 -R -O -K -V >> 0620.ps
awk '{print $3,$4}' temp3.d | psxy -JX6i/3i -R0/250/-5000/-600 -Sc0.1c
-B20/500/SW -Y7i -V -O >> 0620.ps
3.太魯閣附近地區地震三維分佈圖
由一些研究文獻可得知，河流含沙量的供應大多來自於山區的崩塌事件，
崩塌為地表侵蝕的重要機制，而崩塌的主要起因為暴雨或地震，
台灣地區位於板塊交界帶，每年有頻繁的地震，
同時也位於颱風路徑所必經的位置，每年因颱風帶來大量的降雨。
將1986年~2008年3~8級的發生於太魯閣地區附近，
範圍121/123/23/25的地震資料，以3D的方式呈現。
![]()
由圖可見，太魯閣附近地區的地震同樣非常頻繁，
規模並沒有一個特定的分佈，地震深度大多不超過50公里，
較深層的地震多靠近較北方的位置。
程式碼
awk '{if (NR<1599) print $7, $6, -$8, -$8, 0.0035*$9^3}' test.txt > eq.txt
makecpt -Crainbow -T-350/0/3 -V -Z > eq.cpt
psbasemap -JM4i -Jz0.02c -R121/123/23/25/-350/0 -E150/40 -B1:Longitude:/1:Latitude:/30:Depth::."Taroko 1986~2008 4<Mw<8 3-DPLOT":10 -X2i -V -K > eq.ps
psxyz eq.txt -JM4i -Jz0.02c -R121/123/23/25/-350/0 -E150/40 -Sc -Ceq.cpt -V -K -O >> eq.ps
pscoast -JM4i -Jz0.02c -R121/123/23/25/-350/0 -E150/40 -Df -W10 -V -K -O >> eq.ps
pstext taroko.txt -JM4i -Jz0.02c -R121/123/23/25/-350/0 -E150/40 -Gred -V -O -K >> eq.ps
psscale -D20c/8/5/0.25 -Ceq.cpt -B50 -V -X-1.8i -O >> eq.ps
4.太魯閣附近地區地震的震源機制解
由震測海灘球看地震類型，為了避免資料太過於混亂，
僅選取1986~2008規模大於6的地震，範圍121/123/23/25。
![]()
由圖可見，此地區的地震型態並不固定，正、逆、平移、逆衝斷層皆有。
程式碼
ps=taroko_cmt.ps
pscoast -R121/123/23/25 -JM6i -B0.5:."Taroko 1986~2008 Mw>6 CMT":
a1g0.5f0.5/a1g0.5f0.5WSen \
-X2.0 -Y2.0 -Df -P -K -V -S100/100/200 -G100/205/100 -W1/0/255/0 > $ps
psmeca tarokocmt.txt -JM6i -R121/123/23/25 -Sd1.2c/12 -G200/0/50 -H -V -K -O >> $ps
pstext taroko.txt -R121/123/23/25 -J -Gred -V -O >> $ps
中央山脈東翼太魯閣地區因水樣豐沛的立霧溪及其支流持續不斷切割岩層，
加上板塊擠壓形成的抬升作用，於此處形成獨特的峽谷景觀。
立霧溪為太魯閣國家公園主要河川，長55公里，流域面積616平方公里，
利用台電綠水觀測站所記錄的立霧溪水文資料，
配合matlab和gmt工具，進行資料的分析工作。
首先以GMT繪製研究區域相關簡圖。
1.台灣島附近海底地形、斷層位置及地震分佈圖
使用地形為ETOPO1，繪製範圍117/125/19/27，
將海岸山脈斷層、琉球海溝、馬尼拉海溝、呂宋島弧、
梨山斷層、西部麓山戴變形前緣以斷層線標記，
黃點為1986年~2008年在範圍內，規模超過6級的地震震央，
(因為水文資料含沙量的資料區間為1986~2008，故地震資料也取相同區間)
地震的分佈多靠近台灣東部或外海，海岸山脈至和平之間是最頻繁的區域。
由海底地形可得知，台灣西部台灣海峽是一個平坦的大陸棚，
東部海岸的高度落差就大很多，屬於斷層海岸。
程式碼
gmtset PLOT_DEGREE_FORMAT ddd:mm:ssF
makecpt -Crainbow -T-7000/4000/200 -Z > 1.cpt
grdgradient ETOPO1_taiwan.grd -A310 -GETOPO1.int -Ne1.0 -M -V
grdimage ETOPO1_taiwan.grd -R117/125/19/27 -JM6i -C1.cpt -IETOPO1.int
-Bf2a2/f2a2WSne:."1986-2008 earthquakes 6<Mw": -V -P -K -Y3i > 1.ps
pscoast -R117/125/19/27 -J -Df -A50 -O -P -V -W0.7p -K >> 1.ps
grdcontour ETOPO1_taiwan.grd -C200 -J -L-15000/0 -R -W0.5p -O -K -V >> 1.ps
psxy 1.txt -R -JM -Gyellow -Sc0.25c -W0.1p -O -K >> 1.ps
pstext taroko.txt -R117/125/19/27 -J -Gred -V -O -K >> 1.ps
psxy ryukyu_trench.asc -R -J -Sf0.5i/0.1ilt -Gblack -Wthick -O -V -K>> 1.ps
psxy south_eastern_fault.asc -R -J -Sf0.5i/0.1ilt -Gblack -Wthick -K -O -V >> 1.ps
psxy deformation_front.asc -R -J -Sf0.5i/0.1irt -Gblack -Wthick -K -O -V >> 1.ps
psxy malina_trench.asc -R -J -Sf0.5i/0.1ilt -Gblack -Wthick -K -O -V >> 1.ps
psxy Lishan_laonung_fault.asc -R -J -Sf0.5i/0.1irt:0.5i -Gblack -Wthick -K -O -V >> 1.ps
psxy coastal_range_fault.asc -R -J -Sf0.5i/0.1irt -Gblack -Wthick -K -O -V >> 1.ps
psscale -D3i/2i/4.5i/0.25ih -O -C1.cpt -I -B2000:BATHYMETRY:/:m: -Y-3i >> 1.ps
2.立霧溪外海海底地形
立霧溪為水量豐沛之溪流，水流亦攜帶大量沉積物，
河流流出山區後形成一個標準形狀的沖積扇，但面積不大，
代表河流所攜帶的沉積物被帶往其他地方沉積，
因此繪製立霧溪出海口外海的地形剖面，
觀察此處是否有海底沖積扇的地形產生。
上圖為海地地形剖面，下圖為剖面位置，
由剖面可見至10公里外海的海岸地形皆為連續的陡峭下坡，
10公里處一個隆起後繼續更為陡峭的下坡，
發現此處並無明顯的海底沖積扇存在。
那麼立霧溪所攜帶的大量沉積物被攜帶至何處?
是因為海底地形原本就非常陡峭，所以無法形成堆積地形，
或者被沿岸流、洋流帶至其他地方堆積，值得繼續探討。
程式碼
makecpt -Csealand -T-10000/4000/500 -V -Z > 0620.cpt
project -C121.65/24.14 -E123.48/23.52 -G0.5 -N -Q > temp1.d
grdtrack temp1.d -GTaiDBv624.grd -V > temp3.d
grdgradient TaiDBv624.grd -TaiDBv624_shade.grd -A0/30 -Ne0.5 -V
grdimage -JM6i -R119/125/21/27 TaiDBv624.grd -TaiDBv624_shade.grd
-C0620.cpt -B2nwSE -P -V -K > 0620.ps
cat <<END> profile.dat
121.65 24.14
123.48 23.52
END
psxy -JM6i -R119/125/21/27 profile.dat -W10/250/0/0 -V -K -O >> 0620.ps
pstext taroko.txt -R119/125/21/27 -J -G -V -O -K >> 0620.ps
pscoast -JM6i -R119/125/21/27 -Df -W5 -V -K -O >> 0620.ps
grdcontour TaiDBv624.grd -C500 -J -A500 -L-15000/0 -R -O -K -V >> 0620.ps
awk '{print $3,$4}' temp3.d | psxy -JX6i/3i -R0/250/-5000/-600 -Sc0.1c
-B20/500/SW -Y7i -V -O >> 0620.ps
3.太魯閣附近地區地震三維分佈圖
由一些研究文獻可得知，河流含沙量的供應大多來自於山區的崩塌事件，
崩塌為地表侵蝕的重要機制，而崩塌的主要起因為暴雨或地震，
台灣地區位於板塊交界帶，每年有頻繁的地震，
同時也位於颱風路徑所必經的位置，每年因颱風帶來大量的降雨。
將1986年~2008年3~8級的發生於太魯閣地區附近，
範圍121/123/23/25的地震資料，以3D的方式呈現。
由圖可見，太魯閣附近地區的地震同樣非常頻繁，
規模並沒有一個特定的分佈，地震深度大多不超過50公里，
較深層的地震多靠近較北方的位置。
程式碼
awk '{if (NR<1599) print $7, $6, -$8, -$8, 0.0035*$9^3}' test.txt > eq.txt
makecpt -Crainbow -T-350/0/3 -V -Z > eq.cpt
psbasemap -JM4i -Jz0.02c -R121/123/23/25/-350/0 -E150/40 -B1:Longitude:/1:Latitude:/30:Depth::."Taroko 1986~2008 4<Mw<8 3-DPLOT":10 -X2i -V -K > eq.ps
psxyz eq.txt -JM4i -Jz0.02c -R121/123/23/25/-350/0 -E150/40 -Sc -Ceq.cpt -V -K -O >> eq.ps
pscoast -JM4i -Jz0.02c -R121/123/23/25/-350/0 -E150/40 -Df -W10 -V -K -O >> eq.ps
pstext taroko.txt -JM4i -Jz0.02c -R121/123/23/25/-350/0 -E150/40 -Gred -V -O -K >> eq.ps
psscale -D20c/8/5/0.25 -Ceq.cpt -B50 -V -X-1.8i -O >> eq.ps
4.太魯閣附近地區地震的震源機制解
由震測海灘球看地震類型，為了避免資料太過於混亂，
僅選取1986~2008規模大於6的地震，範圍121/123/23/25。
由圖可見，此地區的地震型態並不固定，正、逆、平移、逆衝斷層皆有。
程式碼
ps=taroko_cmt.ps
pscoast -R121/123/23/25 -JM6i -B0.5:."Taroko 1986~2008 Mw>6 CMT":
a1g0.5f0.5/a1g0.5f0.5WSen \
-X2.0 -Y2.0 -Df -P -K -V -S100/100/200 -G100/205/100 -W1/0/255/0 > $ps
psmeca tarokocmt.txt -JM6i -R121/123/23/25 -Sd1.2c/12 -G200/0/50 -H -V -K -O >> $ps
pstext taroko.txt -R121/123/23/25 -J -Gred -V -O >> $ps
一、處理資料簡介
研究區域設定為花蓮縣北方立霧溪，嘗試由水文資料看立霧溪的逕流特性，
例如-流量變化情形、水流量與含沙量之相對應關係...等。
資料的來源為台電公司綠水觀測站，此站水文資料記錄開始於1956年，
資料內容包含流量、含沙量、水位、流速、水溫，
自1956年1月1日起有每日流量之記錄，1981年起增加時水位及時流量記錄，
1986年開始不定期記錄含沙量、水溫、流速，記錄沒有固定的記錄間距，
2000年前的記錄較頻繁，2000~2008年的記錄則較缺乏。
分析的資料主要擷取1959年至2008的每日日流量，
將每日資料平均計算為月平均，
因儀器損毀而缺記錄的月份，以50年的當月平均取代，
50年總計600筆月平均資料。
另外取1986~2008的含沙量資料做分析，
測定含沙量時會同時記錄水量、水溫、流速、水位，
選取水量及含沙量的部分，23年總計884筆。
例如-流量變化情形、水流量與含沙量之相對應關係...等。
資料的來源為台電公司綠水觀測站，此站水文資料記錄開始於1956年，
資料內容包含流量、含沙量、水位、流速、水溫，
自1956年1月1日起有每日流量之記錄，1981年起增加時水位及時流量記錄，
1986年開始不定期記錄含沙量、水溫、流速，記錄沒有固定的記錄間距，
2000年前的記錄較頻繁，2000~2008年的記錄則較缺乏。
分析的資料主要擷取1959年至2008的每日日流量，
將每日資料平均計算為月平均，
因儀器損毀而缺記錄的月份，以50年的當月平均取代，
50年總計600筆月平均資料。
另外取1986~2008的含沙量資料做分析，
測定含沙量時會同時記錄水量、水溫、流速、水位，
選取水量及含沙量的部分，23年總計884筆。
訂閱：
意見 (Atom)