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

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 來算

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