拓端tecdat|R語(yǔ)言非線性回歸nls探索分析河流階段性流量數(shù)據(jù)和評(píng)級(jí)曲線、流量預(yù)測(cè)可視
原文鏈接:http://tecdat.cn/?p=24761?
原文出處:拓端數(shù)據(jù)部落公眾號(hào)
本文檔通過(guò)一些探索性數(shù)據(jù)分析來(lái)制定河流的評(píng)級(jí)曲線和流量預(yù)測(cè)。目的是利用 (1) 在底部安裝單元的定期部署期間測(cè)量的瞬時(shí)流量和 (2) 來(lái)自長(zhǎng)期部署在河流中的水位數(shù)據(jù)記錄器的瞬時(shí)深度測(cè)量,以創(chuàng)建和更新評(píng)級(jí)曲線。額定曲線將用于計(jì)算 HOBO 壓力傳感器部署期間(大約 1 年)的流量。所得數(shù)據(jù)將用于創(chuàng)建和驗(yàn)證河流 10-15 年期間的回歸和 DAR 流量估計(jì)。
介紹
評(píng)級(jí)曲線
由于與河流流量連續(xù)測(cè)量相關(guān)的高成本,最好使用河流高度測(cè)量來(lái)估計(jì)流量。使用壓力傳感器可以連續(xù)測(cè)量水流高度。當(dāng)輔以周期性的流量測(cè)量時(shí),冪函數(shù)可以關(guān)聯(lián)河流高度和流量(Venetis):
其中:Q代表穩(wěn)態(tài)排放,H代表流高(階段),H0是零排放階段;K 和 z 是評(píng)級(jí)曲線常數(shù)。按照慣例,Q 和 H 通常在參數(shù)估計(jì)之前進(jìn)行對(duì)數(shù)變換。
當(dāng)河流水位過(guò)程線的上升和下降階段導(dǎo)致相同河流高度的不同流量時(shí),就會(huì)發(fā)生不穩(wěn)定流。由此產(chǎn)生的受滯后影響的評(píng)級(jí)曲線將呈現(xiàn)為一個(gè)循環(huán)而不是一條線。Petersen-?verleir 描述的修正瓊斯公式和 Zakwan?可能用過(guò)了:
其中 Q?是流量,h 是河流高度。偏一階導(dǎo)數(shù)?
?使用有限差分近似為 J:
?
其中 ht 是時(shí)間 t 的水流高度,Δt 是時(shí)間區(qū)間。這可以被認(rèn)為是河流高度和時(shí)間之間函數(shù)的斜率或瞬時(shí)變化率,它是使用測(cè)量的河流高度值估計(jì)的。K、a、n?和 x?是評(píng)級(jí)曲線常數(shù)。
許多不同的方法可用于求解額定曲線參數(shù)。我們使用非線性最小二乘回歸來(lái)最小化評(píng)級(jí)曲線參數(shù)的殘差平方和 (SSE)。殘差 SSE 計(jì)算如下:
?
其中:X?是測(cè)量值,Y?是預(yù)測(cè)值。非線性優(yōu)化方法搜索參數(shù)組合以最小化目標(biāo)函數(shù)(在這種情況下為殘差 SSE)。彼得森 應(yīng)用 Nelder-Mead 算法求解瓊斯公式。扎關(guān)?使用廣義減少梯度和遺傳算法提出非線性優(yōu)化方法。大多數(shù)方法需要仔細(xì)規(guī)劃有點(diǎn)接近全局最小值的參數(shù)起始值,或者存在識(shí)別替代局部最小值的風(fēng)險(xiǎn)。為了減少局部最小值收斂的可能性,?R
?提供了在許多不同的起始值上迭代非線性最小二乘優(yōu)化的功能(Padfield 和 Matheson).
未控制的流量估計(jì)
評(píng)級(jí)曲線允許在部署水流深度數(shù)據(jù)記錄器的時(shí)間段內(nèi)開發(fā)每日水流記錄。然而,當(dāng)站點(diǎn)未啟用時(shí),對(duì)每日流量的估計(jì)需要額外的信息。統(tǒng)計(jì)信息傳遞和經(jīng)驗(yàn)回歸是兩種相對(duì)簡(jiǎn)單的方法,可用于估算測(cè)量不當(dāng)?shù)牧饔蛑械牧髁俊=y(tǒng)計(jì)傳輸程序使用面積和徑流之間的假設(shè)關(guān)系,簡(jiǎn)單地將流量持續(xù)時(shí)間曲線或每日流量值從有測(cè)量的流域傳輸?shù)轿礈y(cè)量的流域。最常用的統(tǒng)計(jì)傳遞方法是流域面積比。使用流域面積比,通過(guò)將面積比與日流量相乘,日流量從一個(gè)流域轉(zhuǎn)移到另一個(gè)流域:
其中,?
?是未測(cè)量盆地 y 和時(shí)間 t 的流量,?
?是測(cè)量盆地 x 和時(shí)間 tt 處的流量,和?
?是盆地的面積比。參數(shù) ? 通常等于 1(Asquith、Roussel 和 Vrabel?)。然而,阿斯奎斯、魯塞爾和弗拉貝爾?提供了在德克薩斯州應(yīng)用時(shí)用于流域面積比的 ? 的經(jīng)驗(yàn)估計(jì)值。有了可用的短期流量記錄,可以使用排水面積比方法評(píng)估各種流量?jī)x表的性能。此外,可以使用非線性最小二乘法開發(fā) ? 的局部值。如果主要輸出是流量持續(xù)時(shí)間曲線,則主要關(guān)注的是候選量具有相似的徑流因變量并且在未治理流域的合理距離內(nèi)。但是,如果主要輸出包括每日流量估計(jì),則具有具有相同流量超出概率時(shí)間的候選量具更為重要。
基于經(jīng)驗(yàn)回歸的方法需要一段時(shí)間的測(cè)量流量和一些預(yù)測(cè)變量來(lái)估計(jì)徑流因變量。通常,使用日降雨量數(shù)據(jù)將回歸模型擬合到測(cè)量的流量數(shù)據(jù):
其中 Qi是第 i?天的預(yù)測(cè)排放量,β?是第 j?個(gè)變量的系數(shù),x 是第 i 天的預(yù)測(cè)變量值。假設(shè)誤差項(xiàng) ?i 正態(tài)分布在均值零附近。使用簡(jiǎn)單線性或多元線性回歸 Q通常在估計(jì)回歸系數(shù)之前進(jìn)行對(duì)數(shù)變換。如果預(yù)測(cè)變量和因變量之間的關(guān)系預(yù)期為非線性多項(xiàng)式,則可以包括項(xiàng)。然而,稱為廣義加性模型的線性回歸的擴(kuò)展允許將這些非線性項(xiàng)相對(duì)容易地?cái)M合到數(shù)據(jù)中。對(duì)于廣義加性模型,因變量取決于應(yīng)用于每個(gè)預(yù)測(cè)變量的平滑函數(shù)的總和。此外,廣義加性模型可以擬合具有非正態(tài)分布的誤差分布的因變量。然而,與線性或多元線性回歸相比,廣義加性模型由于缺乏單一模型系數(shù)而更難以解釋。因此,每個(gè)單獨(dú)的平滑函數(shù)對(duì)因變量均值的影響通常以圖形方式傳達(dá)。
方法
數(shù)據(jù)采集
數(shù)據(jù)來(lái)源于水位數(shù)據(jù)記錄器。部署了一個(gè)額外的數(shù)據(jù)記錄器,為部署在水下的數(shù)據(jù)記錄器提供環(huán)境大氣壓力校正。從 2020-03-02 到 2021-..-...,水位數(shù)據(jù)記錄器幾乎連續(xù)部署,并設(shè)置為每隔 15 分鐘記錄一次水位。
## 制作要導(dǎo)入的文件列表
list.files(path = here("Data
##創(chuàng)建一個(gè)空白的tibble來(lái)填充
tibble()
## 遍歷文件路徑以讀取每個(gè)文件
for (i in fihs) {
x <- read_csv(
copes
bind_rows(hf, x)
rm(x)
表 1:在每個(gè)站點(diǎn)測(cè)量的 15 分鐘流級(jí)別的匯總統(tǒng)計(jì)數(shù)據(jù)
ggplot(hdf) +
geom_line+
facet_wrap +
scale_x_datetim+
scale_y +
guid +
th +
theme
圖 1:每個(gè)站點(diǎn)的 15 分鐘區(qū)間流深度記錄。
使用底部安裝的流量計(jì)在每個(gè) SWQM 站點(diǎn)進(jìn)行定期流量測(cè)量。測(cè)量橫截面積、水流高度和速度。使用這些測(cè)量值,該設(shè)備利用指數(shù)速度方法來(lái)報(bào)告瞬時(shí)流量。流量測(cè)量設(shè)備一次部署幾天,在每個(gè)站的不同流量條件下捕獲完整的水文過(guò)程線。只有兩個(gè)流量計(jì)可用,因此在站點(diǎn)之間輪流部署。此外,一臺(tái)設(shè)備停止工作并進(jìn)行了幾個(gè)月的維修。以 15 分鐘的間隔記錄流量。
在數(shù)據(jù)探索過(guò)程中,每個(gè)站點(diǎn)的低流量數(shù)據(jù)中明顯存在過(guò)多噪聲。在停滯或接近停滯條件期間,多普勒流量計(jì)記錄高度可變的流速并報(bào)告不切實(shí)際的流量。由于過(guò)多的數(shù)據(jù)噪聲,從數(shù)據(jù)記錄中清除了極低或停滯的流量時(shí)期。未來(lái)的部署將需要考慮在什么條件下長(zhǎng)期部署是合適的。對(duì)于像這樣的小流,定期的風(fēng)暴流部署可能是最合適的部署。
## 制作要導(dǎo)入的文件列表
file_paths <- paste0(he ".csv"))
##創(chuàng)建一個(gè)空白的tibble來(lái)填充
iq <- tibble()
## 遍歷文件路徑以讀取每個(gè)文件
for (i in file_paths) {
x <- read_csv
file <- i
idf <- bind_rows
rm
}
iqdf <- iqdf %>%
## 使用 `dplyr::` 指定要使用的重命名函數(shù),以防萬(wàn)一
dplyr::rename(Sam)
ggplot(iqdf)+
geom_point(aes(Dme, Flow), alpha = 0.2) +
facet_wrap
圖 2:每個(gè)站點(diǎn)記錄數(shù)據(jù)的周期。
## 為了將測(cè)量深度與IQ的流速測(cè)量結(jié)合起來(lái)
## ##我們需要插值測(cè)量深度到每分鐘,因?yàn)樯疃仁瞧?。然后我們就可以連接這些數(shù)據(jù)。我們將使用線性插值。
##使用purrr::map在每個(gè)站點(diǎn)上運(yùn)行插值運(yùn)算
hdf %>%
split%>%
map %>%
bind_row %>%
as_tibble
##這就是我們要開發(fā)評(píng)級(jí)曲線的數(shù)據(jù)框架。
idf %>%
left_join %>%
select %>%
mutate -> iqdf
評(píng)級(jí)曲線發(fā)展
由于植物生長(zhǎng)、植物枯死、沉積、侵蝕和其他過(guò)程,河內(nèi)和河岸的條件全年都會(huì)發(fā)生變化。這些不斷變化的條件可能需要開發(fā)多個(gè)評(píng)級(jí)曲線。探索性數(shù)據(jù)分析用于確定變化時(shí)期和滯后影響評(píng)級(jí)曲線的可能性。一旦確定了評(píng)級(jí)曲線周期和適當(dāng)?shù)墓?,公式中的評(píng)級(jí)曲線參數(shù)?(1)?和?(2)?通過(guò)非線性最小二乘估計(jì)回歸使用??R
?(Padfield?)。該方法利用 Levenberg-Marquardt 算法和多個(gè)起始值來(lái)尋找全局最小 SSE 值。
單獨(dú)的評(píng)級(jí)曲線用于使用測(cè)量的河流高度估計(jì)河流流量。Nash-Sutcliffe 效率 (NSE) 和歸一化均方根誤差用于評(píng)估測(cè)量和估計(jì)流量之間的擬合優(yōu)度。NSE 是歸一化統(tǒng)計(jì)量,用于評(píng)估相對(duì)于測(cè)量數(shù)據(jù)方差的相對(duì)殘差方差,計(jì)算公式如下:
?
其中?
?是觀察到的排放量的平均值,?
?是 t時(shí)刻的估計(jì)流量,Qt 是 t時(shí)刻觀察到的流量。NSE 的值范圍從 ?∞ 到 1,其中 1 表示完美的預(yù)測(cè)性能。NSE 為零表示模型具有與數(shù)據(jù)集均值相同的預(yù)測(cè)性能。
nRMSE 是一個(gè)基于百分比的指標(biāo),用于描述預(yù)測(cè)和測(cè)量的排放值之間的差異:
其中
?其中 Qt 是在時(shí)間 t 觀察到的流量,?
?是 t 時(shí)刻的估計(jì)排放量,n是樣本數(shù),?
?和?
?是觀察到的最大和最小排放量。產(chǎn)生的 nRMSE 計(jì)算是一個(gè)百分比值。
結(jié)果
站點(diǎn)?
基于探索性分析,為站點(diǎn)制定了兩條評(píng)級(jí)曲線。評(píng)級(jí)曲線周期為2020-03-03至2020-11-30和2020-12-01至2021-01-31。由于觀察到的水層中存在明顯的不穩(wěn)定流動(dòng),我們應(yīng)用了瓊斯公式(公式(2))。兩個(gè)時(shí)間段都產(chǎn)生了 NSE 大于 0.97 和 nRMSE 小于 3% 的評(píng)級(jí)曲線,表明非常適合(表?2; 數(shù)字?3)。數(shù)字?3?確實(shí)表明在極低流量測(cè)量中存在一些有偏差的流量估計(jì)。這歸因于多普勒流量計(jì)在低流量時(shí)記錄的流量變化。
## 為站點(diǎn) 制作數(shù)據(jù)框
if %>%
group_split %>%
## 刪除最大流量未超過(guò) 10 cfs 的事件
imap %>%
bind_rows
## 為站點(diǎn) 2020 年 12 月 14 日至 2021 年 1 月制作數(shù)據(jù)框
idf %>%
diff_time %>%
group_split %>%
imap(~mutate) %>%
bind_rows %>%
##使用nls估計(jì)瓊斯公式中的參數(shù)
## 一些起始參數(shù)。 nls_multstart 將使用多個(gè)
##起始參數(shù)和模型選擇查找
##全局最小值
stlower
stupper
##適合nls
rc<- nls(jorm,
suors = "Y")
##設(shè)置參數(shù)起始限制
##nls
nls(jorm
suors = "Y")
## 使用參數(shù)結(jié)果設(shè)置數(shù)據(jù)框。 將使用它來(lái)報(bào)告參數(shù)和 GOF 指標(biāo)
dfresu <- tibble(Site ?)
## 開發(fā)評(píng)級(jí)曲線預(yù)測(cè)和
## 使用 GOF 指標(biāo)創(chuàng)建表
df %>%
dfres %>%
NSE(predicted, as.numeric(Flow))-> dfres
##顯示表
kable(dfres)
表 2:站點(diǎn) 的評(píng)級(jí)曲線參數(shù)估計(jì)和擬合優(yōu)度指標(biāo)
##繪制評(píng)級(jí)曲線結(jié)果
p1 <- ggplot +
geom_pointlpha = 0.25) +
geom_abline +
scale_x_continuous
theme()
圖 3:站點(diǎn) 處每個(gè)額定曲線周期的額定曲線估計(jì)流量與實(shí)測(cè)流量(A、C)和階段流量預(yù)測(cè)(B、D)的散點(diǎn)圖。
站點(diǎn) 16397
探索性分析表明,在該站點(diǎn)使用冪函數(shù)(公式?(1)) 因?yàn)樵谒倪^(guò)程圖中沒(méi)有觀察到不穩(wěn)定的流動(dòng)條件。評(píng)級(jí)曲線預(yù)測(cè)導(dǎo)致 NSE 大于 0.95,表明非常適合(表?2)。nRMSE 小于 5%,這對(duì)于在該站獲得的較小樣本量來(lái)說(shuō)可能是一個(gè)很好的結(jié)果,并且可能受到觀察到的低流量方差的影響(表?2; 圖?3).
## 設(shè)置數(shù)據(jù)框以將評(píng)級(jí)曲線擬合到 1697
##冪函數(shù)
##參數(shù)起始限制
rc_7 <- nls(prm,
stalower
supors = "Y")
## 使用參數(shù)結(jié)果設(shè)置數(shù)據(jù)框。 稍后將使用它來(lái)報(bào)告參數(shù)和 GOF 指標(biāo)
dfres <- tibble(Sit,
K,
H0 ,
Z["Z"]])
df167 %>%
mutate(pr = exp(predict(rc))) -> df
表 3:站點(diǎn) 167 的評(píng)級(jí)曲線參數(shù)估計(jì)和擬合優(yōu)度指標(biāo)
##繪制評(píng)級(jí)曲線結(jié)果
p1 <- df7 %>%
ggplot() +
geom_point +
geom_abline+
scale_x_continuous +
scale_y_continuous
p1 + p2 + plot_annotation
圖 4:167 站的額定流量估計(jì)流量與實(shí)測(cè)流量 (A) 和階段流量預(yù)測(cè) (B) 的散點(diǎn)圖。
站點(diǎn) 162
基于探索性分析,為站點(diǎn)166制定了3條評(píng)級(jí)曲線。評(píng)級(jí)曲線周期為2020-03-03至2020-05-30;2020-06-01 至 2020-10-31;和 2020-11-01 至 2021-01-31。由于觀察到的水層中存在明顯的不穩(wěn)定流動(dòng),我們應(yīng)用了瓊斯公式(公式(2))。3 月至 9 月的結(jié)果表明評(píng)級(jí)曲線具有非常好的擬合(NSE > 0.96,nRMSE <6%;表?4)。低流量下觀測(cè)值和預(yù)測(cè)值之間的巨大差異可歸因于具有極快的水流高度變化(> 1.5 英尺/小時(shí))的事件,參數(shù)估計(jì)難以擬合(圖?5?)。其余評(píng)級(jí)曲線的擬合優(yōu)度指標(biāo)有所下降,但仍表明性能良好(表?4)。測(cè)得的中低流量值的高方差影響評(píng)級(jí)曲線性能(圖?5).
## 制作 3 個(gè)不同的數(shù)據(jù)框以適應(yīng)瓊斯公式。
iqpdf %>%
filter(DaTime > as.POSIXct)
#filter %>%
group_split %>%
imap %>%
bind_rows
## 估計(jì)每個(gè)數(shù)據(jù)集的參數(shù)
##設(shè)置參數(shù)起始限制
stlower
stupper
rc<- nls(jonrm
surs = "Y")
##設(shè)置參數(shù)起始限制
stlower
stupper
rc <- nls(jorm,
stawer,
suprs = "Y")
##設(shè)置參數(shù)起始限制
stalower
## 使用參數(shù)結(jié)果設(shè)置數(shù)據(jù)框。稍后將使用它來(lái)報(bào)告參數(shù)和 GOF 指標(biāo)
dfres <- tibble(Site
K ,
a ,
n ,
x )
## 開發(fā)評(píng)級(jí)曲線預(yù)測(cè)和
## 使用 GOF 指標(biāo)創(chuàng)建表
##顯示表
kable(df)
表 4:站點(diǎn) 16882 的評(píng)級(jí)曲線參數(shù)估計(jì)和擬合優(yōu)度指標(biāo)
##繪制評(píng)級(jí)曲線結(jié)果
p1 <- ggplot(df_03) +
geom_point +
geom_abline+
scale_x_continuous +
scale_y_continuous +
theme_ms()
圖 5:16882 站每個(gè)額定曲線周期的額定曲線估計(jì)流量與實(shí)測(cè)流量(A、C、E)和階段流量預(yù)測(cè)(B、D、F)的散點(diǎn)圖。
每日流量估算
# 使用原始數(shù)據(jù)集
# 按日期使用評(píng)級(jí)曲線估計(jì)流量
# 聚合表示每日流量,報(bào)告匯總統(tǒng)計(jì)數(shù)據(jù)。
hodf %>%
dplyr::select%>%
group_split(站點(diǎn)) %>%
bind_rows()
## 制作模型的數(shù)據(jù)框,預(yù)測(cè)數(shù)據(jù),然后映射預(yù)測(cè)函數(shù),并取消嵌套數(shù)據(jù)框。
tibble) %>%
~exp( newdata = .y))
)) %>%
tidyr::unnest%>%
as_tsibble
##繪制數(shù)據(jù)
ggplot() +
geom_line+
geom_point +
facet_wrap +
scale_y_continuous +
theme_ms()
圖 6:15 分鐘流量估算與疊加測(cè)量流量值
## 計(jì)算平均每日流量并報(bào)告統(tǒng)計(jì)數(shù)據(jù)
dplyr::select(Site, DTime) %>%
group_by_key() %>%
index_by(date = ~as_date(.)) %>%
## 報(bào)告摘要統(tǒng)計(jì)
meflow %>%
as_tibble() %>%
dplyr::select %>%
tbl_summary %>%
as_kable()
表 5:每個(gè)站點(diǎn)平均日流量估計(jì)的匯總統(tǒng)計(jì)數(shù)據(jù)
ggplot(mean_daily_flow %>% fill_gaps()) +
geom_line(aes(date, mean_daily)) +
facet_wrap(~Site, ncol = 1, scales = "free_y") +
scale_x_date("Date", date_breaks = "1 month",
date_labels = "%F") +
labs(y = "Mean daily streamflow [cfs]") +
theme_ms() +
guides(x = guide_axis(angle = 90)) +
theme(axis.text.x = element_text(size = 8))
圖 7:平均每日流量估計(jì)。
##導(dǎo)出數(shù)據(jù)
melow %>%
writcsv("modf.csv")
參考
Asquith, William H., Meghan C. Roussel, and Joseph Vrabel. 2006. “Statewide Analysis of the Drainage-Area Ratio Method for 34 Streamflow Percentile Ranges in Texas.” 2006-5286. U.S. Geological Survey.
最受歡迎的見解
1.R語(yǔ)言基于ARMA-GARCH-VaR模型擬合和預(yù)測(cè)實(shí)證研究
2.R語(yǔ)言時(shí)變參數(shù)VAR隨機(jī)模型
3.R語(yǔ)言估計(jì)時(shí)變VAR模型時(shí)間序列的實(shí)證研究
4.R語(yǔ)言基于ARMA-GARCH過(guò)程的VAR擬合和預(yù)測(cè)
5.GARCH(1,1),MA以及歷史模擬法的VaR比較
6.R語(yǔ)言用向量自回歸(VAR)進(jìn)行經(jīng)濟(jì)數(shù)據(jù)脈沖響應(yīng)
7.R語(yǔ)言實(shí)現(xiàn)向量自動(dòng)回歸VAR模型
8.R語(yǔ)言隨機(jī)搜索變量選擇SSVS估計(jì)貝葉斯向量自回歸(BVAR)模型
9.R語(yǔ)言VAR模型的不同類型的脈沖響應(yīng)分析