四、水文資料分析(二)

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含沙量。