🥦Python流体数据统计模型和浅水渗流平流模型模拟

关键词

Python | 流体 | 数据 | 统计 | 数学 | 拟合 | 模型 | 算法 | 曲线 | 回归 | 预测 | 水位 | 绘图 | 湿度 | 湿度 | 三维 | 二维 | 非线性 | 平稳性 | 趋势 | 季节性 | 周期性 | 序列 | 时间 | 分量 | 误差 | 显着性测试 | 方差 | 估计 | 曼宁方程 | 明渠 | 流体动力 | 矩形 | 三角形 | 圆形 | 偏微分 | 方程 | 粗糙 | 坡度 | 有限体积 | 浅水 | 渗流 | 平流 | 畜水层 | 压力 | 速度 | 网格

Python流体数据统计模型和浅水渗流平流模型模拟 | Notion亚图跨际 on Notion

🏈指点迷津 | Brief

🎯要点

1. 流体数据统计模型：

   1. 🎯曲线拟合和回归模型预测水流：🖊建立流量水位回归模型，模型拟合，绘制结果，🖊以为水位相对湿度的函数，建立多线性回归模型，绘制三维结果。🖊以水位，水位和相对湿度乘积的函数，建立非线性回归模型，绘制三维结果。

   2. 🎯平稳性、趋势、季节性和周期性数据序列分析：🖊时间戳的水流数据分解为趋势和季节性分量，使用增强迪基富勒测试水流数据平稳性，绘制二维趋势图，🖊建立水流自回归、自回归移动平均线、自回归综合移动平均线和简单指数平滑模型，使用增强迪基富勒测试测试平稳性，均方根误差量化预测误差，绘制二维趋势图评估有效性。

   3. 🎯显着性测试两处水文流量差异方法：🖊单因素和双因素方差分析，🖊t-测试，F-测试，柯尔莫哥洛夫-斯米尔诺夫检验，曼惠特尼测试。

   4. 🎯估计不确定性：🖊非特定分布区间下，估计可能水文下降的范围，🖊自举置信区间估计，生成可替换多重采样水文，🖊对称置信区间估计一维水流，🖊分位数非参数置信区间估计和预测水文范围，🖊分位数回归推理和绘图检测水文超出指定分位数的异常值，🖊蒙特卡罗不确定性传播模拟曼宁方程，估计流体动力明渠中的流速，绘制二维图示。

2. 流体​模型​：

   1. 🎯模拟矩形、三角形、圆形通道运动学波动偏微分方程：🖊设置库朗弗里德里希-路易条件，定义粗糙系数，坡度，降雨量或侧向流入量，流入时间，绘制排放量趋势图。

   2. 🎯二维浅水方程偏微分方程：🖊有限体积法模拟溃坝，绘制三维流图。

   3. 🎯多孔介质水力传导率及其梯度关联模型：🖊有限体积法模拟，计算压力条件下流体网格，设置狄利克雷边界条件，创建线性求解器，绘制压力和流速趋势图。

   4. 🎯畜水层水流及分布预测模型：🖊有限体积法模拟多孔介质流的数学模型，绘制三维流域图。

   5. 🎯平流传输模型：🖊纳维-斯托克斯方程二维模拟，软件生成网格元素，代码设置边界条件，指定流体运动粘度和密度，绘制二维动画结果。

🍇Python洪泛区潜在危险分类建模和图形显示

溃坝洪水建模的第一步是对要建模的河流长度进行经验计算。 水库中储存的最大水量是确定下游到感受到破裂波影响的距离的最相关因素。 通过分析大坝溃坝的真实案例，建立了一条回归曲线，将蓄水量与潜在受影响区域下游的最大距离联系起来。

为了使该曲线连续，从而实现程序的自动化，施加了 6.7 公里的最小限制，因为在这个距离上至少有 50% 的死亡发生，最多下游 100 公里，被认为足以进行分类的值。 最小和最大限制如图中的虚线所示。

回归由以下方程定义。该方程对于0 $hm ^3$ 和 1000 $hm ^3$ 之间的体积有效，

$$L=8.87^* 10^8 V^3-2.6^* 10^{-4} V^2+0.265 V+6.74$$

其中，

L：发生破裂时受影响的沿河长度（公里）

V：大坝体积（hm）

如果 L<5 则 L=5；如果 L >100 则 L =100。

大坝下游的破裂流和阻尼流计算都是使用文献中提供的经验方程进行的。 对于大坝的最大破裂流，根据大坝高度和体积之间的关系，采用美国陆军工程兵团的弗罗因德利希方程和 MMC ( 绘图、建模和结果生产中心) 方程。

对于与高度相关的高体积，使用 MMC 方程方法：

$$Q_{\max }=0.0039\left(V_{\max }^{0.8122}\right)$$

否则，使用弗罗因德利希方程方法：

$$Q_{\max }=0.607\left(V_{\max }^{0.295} H_{\max }^{1.24}\right)$$

其中，

$Q _{\max }$: 大坝溃坝流量，$m ^3 / s$

$V_{\max }$：大坝体积，$m ^3$

$H _{\max }$：大坝高度，$m$

对于沿山谷的阻尼流，使用上述两个方程。 对于体积大于 6.2 hm3 的大坝，遵循美国垦务局提出的建议，得出以下公式，其中阻尼仅取决于截面距大坝的距离“x” 。

$$Q_x=Q_{\max } 10^{-0.01243 x}$$

其中，

$Q _{\max }$: 大坝溃坝流量，$m ^3 / s$

$x$：距大坝的距离，$m$

$Q _{ x }$ ：距大坝 x 米的流量，单位为 $m ^3 / s$

对于低于 6.2 $hm^3$ 的体积，采用世界银行提出的算盘，其中阻尼取决于距离 x 和水库体积 $V_{\max }$。 根据算盘，定义了以下方程。

$$\frac{Q_x}{Q_{\max }}=a e^{b x}$$

其中，

$\begin{aligned} & a=0.002 \ln \left(V_{\max }\right)+0.9626 \\ & b=-0.20047\left(V_{\max }+25000\right)^{-0.5979}\end{aligned}$

参数a和b是由下图所示的5条曲线通过多参数回归得到的。

Python和PyQt5代码实现

图形界面

 from PyQt5 import QtCore, QtGui, QtWidgets
 import matplotlib.pyplot as plt
 from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas
 from matplotlib.backends.backend_qt5agg import NavigationToolbar2QT as NavigationToolbar
 import pyvista as pv
 import pandas as pd
 import random
 ​
 class Ui_Dialog(object):
     def setupUi(self, Dialog):
         
         self.fig, self.ax = plt.subplots(figsize=(10,10))
         self.canvas = FigureCanvas(self.fig)
         
         self.plotter =QtInteractor()
     
         Dialog.setObjectName("Dialog")
         Dialog.resize(900, 900)
         self.gridLayout = QtWidgets.QGridLayout(Dialog)
         self.gridLayout.setObjectName("gridLayout")
         self.tabWidget = QtWidgets.QTabWidget(Dialog)
         self.tabWidget.setEnabled(True)
         self.tabWidget.setObjectName("tabWidget")
         self.tab = QtWidgets.QWidget()
         self.tab.setObjectName("tab")
         self.gridLayout_3 = QtWidgets.QGridLayout(self.tab)
         self.gridLayout_3.setObjectName("gridLayout_3")
         self.label_3 = QtWidgets.QLabel(self.tab)
         self.label_3.setObjectName("label_3")
         self.gridLayout_3.addWidget(self.label_3, 2, 0, 1, 1)
         self.hval = QtWidgets.QDoubleSpinBox(self.tab)
         self.hval.setMaximum(100000000000.0)
         self.hval.setObjectName("hval")
         self.gridLayout_3.addWidget(self.hval, 2, 2, 1, 1)
         self.label = QtWidgets.QLabel(self.tab)
         self.label.setObjectName("label")
         self.gridLayout_3.addWidget(self.label, 0, 0, 1, 1)
         self.longval = QtWidgets.QDoubleSpinBox(self.tab)
         self.longval.setMaximum(1000000.0)
         self.longval.setMinimum(-1000000.0)
         self.longval.setDecimals(4)
         self.longval.setObjectName("longval")
         self.gridLayout_3.addWidget(self.longval, 1, 2, 1, 1)
         self.label_4 = QtWidgets.QLabel(self.tab)
         self.label_4.setObjectName("label_4")
         self.gridLayout_3.addWidget(self.label_4, 3, 0, 1, 1)
         self.latval = QtWidgets.QDoubleSpinBox(self.tab)
         self.latval.setMaximum(1000000.0)
         self.latval.setMinimum(-1000000.0)
         self.latval.setDecimals(4)
         self.latval.setObjectName("latval")
         self.gridLayout_3.addWidget(self.latval, 0, 2, 1, 1)
         self.tab1calc = QtWidgets.QPushButton(self.tab)
         self.tab1calc.setObjectName("tab1calc")
         self.gridLayout_3.addWidget(self.tab1calc, 4, 2, 1, 1)
         self.label_2 = QtWidgets.QLabel(self.tab)
         self.label_2.setObjectName("label_2")
         self.gridLayout_3.addWidget(self.label_2, 1, 0, 1, 1)
         self.vval = QtWidgets.QDoubleSpinBox(self.tab)
         self.vval.setMaximum(100000000000.0)
         self.vval.setObjectName("vval")
         self.gridLayout_3.addWidget(self.vval, 3, 2, 1, 1)
         self.crioval = QtWidgets.QLabel(self.tab)
         self.crioval.setObjectName("crioval")
         self.gridLayout_3.addWidget(self.crioval, 5, 2, 1, 1)
         self.label_8 = QtWidgets.QLabel(self.tab)
         self.label_8.setObjectName("label_8")
         self.gridLayout_3.addWidget(self.label_8, 5, 0, 1, 1)
         spacerItem = QtWidgets.QSpacerItem(20, 40, QtWidgets.QSizePolicy.Minimum, QtWidgets.QSizePolicy.Expanding)
         self.gridLayout_3.addItem(spacerItem, 6, 0, 1, 1)
         self.tabWidget.addTab(self.tab, "")
         self.tab_2 = QtWidgets.QWidget()
         self.tab_2.setObjectName("tab_2")
         self.gridLayout_4 = QtWidgets.QGridLayout(self.tab_2)
         self.gridLayout_4.setObjectName("gridLayout_4")
         self.nretasval = QtWidgets.QSpinBox(self.tab_2)
         self.nretasval.setEnabled(True)
         self.nretasval.setProperty("value", 8)
         self.nretasval.setObjectName("nretasval")
         self.gridLayout_4.addWidget(self.nretasval, 2, 1, 1, 1)
         self.matplotlibwidget = self.canvas
         self.matplotlibwidget.setObjectName("matplotlibwidget")
         self.gridLayout_4.addWidget(self.matplotlibwidget, 8, 0, 1, 2)
         self.mtoolbar = NavigationToolbar(self.canvas, self.matplotlibwidget)
         self.tab2calc = QtWidgets.QPushButton(self.tab_2)
         self.tab2calc.setEnabled(True)
         self.tab2calc.setObjectName("tab2calc")
         self.gridLayout_4.addWidget(self.tab2calc, 7, 1, 1, 1)
         self.carregartracado = QtWidgets.QPushButton(self.tab_2)
         self.carregartracado.setEnabled(True)
         self.carregartracado.setObjectName("carre")
         self.gridLayout_4.addWidget(self.carregartracado, 1, 1, 1, 1)
         self.carregarsrtm = QtWidgets.QPushButton(self.tab_2)
         self.carregarsrtm.setEnabled(True)
         self.carregarsrtm.setObjectName("carregarsrtm")
         self.gridLayout_4.addWidget(self.carregarsrtm, 0, 1, 1, 1)
 ​
 ​
     def retranslateUi(self, Dialog):
         _translate = QtCore.QCoreApplication.translate
         Dialog.setWindowTitle(_translate("Dialog"))
         self.label_3.setText(_translate("Dialog"))
         self.label.setText(_translate("Dialog", "Latitude"))
         self.label_4.setText(_translate("Dialog", "Volume (m³)"))
         self.tab1calc.setText(_translate("Dialog", "Calcular"))
         self.label_2.setText(_translate("Dialog", "Longitude"))
         self.crioval.setText(_translate("Dialog", "0"))
         self.label_8.setText(_translate("Dialog", "Comprim (km)"))
         self.tabWidget.setTabText(self.tabWidget.indexOf(self.tab), _translate("Dialog", "Pas 1"))
         self.tab2calc.setText(_translate("Dialog", "Calcular"))
         self.carregartracado.setText(_translate("Dialog", "Carregar traçado"))
         self.carregarsrtm.setText(_translate("Dialog", "Carregar SRTM"))
         self.label_5.setText(_translate("Dialog", "Nmero"))
         self.label_6.setText(_translate("Dialog", "Comprime(m)"))
         self.exportsec.setText(_translate("Dialog", "Export"))
         self.importsec.setText(_translate("Dialog", "Import"))
         self.tabWidget.setTabText(self.tabWidget.indexOf(self.tab_2), _translate("Dialog", "Pas 2"))
         self.saveshape.setText(_translate("Dialog", "shape file"))
         self.tab3calc.setText(_translate("Dialog", "Calcular"))
         self.savereport.setText(_translate("Dialog", "Salva"))
         self.tabWidget.setTabText(self.tabWidget.indexOf(self.tab_3), _translate("Dialog", "Pas 3"))
    
     def calcrio(self):
         print('Calculan...')
         longitude = self.longval.value()
         latitude =  self.latval.value()
   
         self.ponto_informado = Point((longitude, latitude))
 ​
         self.h = self.hval.value()
         self.v = self.vval.value()
 ​
         self.criov = crio(self.v)
         self.criov = round(self.criov, 2)
 ​
         self.crioval.setText(str(self.criov))
         
     def carregar_srtm(self):
         file_path = QtWidgets.QFileDialog.getOpenFileName(None, "Selec", "", "tif files (*.tif)")
         srtm_arquivo = file_path[0] 
         self.srtm = rasterio.open(srtm_arquivo)
         print('SRTM!')
         
     def carregar_tracado(self):
         file_path = QtWidgets.QFileDialog.getOpenFileName(None, "Selec", "", "kml files (*.kml)")
         tracado_arquivo = file_path[0] 
         self.tracado = geopandas.read_file(tracado_arquivo, driver='KML')
         print('Tra!')
 ​
     def calcular_perpendiculares(self):
         print('Calcu...')
         self.ax.clear()
         n = self.nretasval.value()
         comp = self.comsecval.value()
 ​
         self.tracado_simplificado = simplificar_tracado(self.tracado, n)
 ​
         self.s, self.ds = secoes_perpendiculares(self.tracado_simplificado, n=21, comprimento=comp)
         self.s.crs = f'EPSG:{datum}'
         st = self.s.to_crs(epsg=4326)
 ​
         
     def exportr_seces(self):
         self.shape_flname = QtWidgets.QFileDialog.getSaveFileName(None, "Selec", "", "shp files (*.shp)")
         self.shape_flname = self.shape_flname[0]
         exportar_geopandas(self.s, nome_do_arquivo=self.shape_flname)
         
     def importar_secoes(self):
         self.s = geopandas.read_file(self.shape_flname)
         
         self.ax.clear()
         show(self.srtm, ax=self.ax)
         self.s.plot(ax=self.ax, color='red')
         self.ax.scatter(self.ponto_informado.x, self.ponto_informado.y, color='red', label='Bar')
         
         self.s = self.s.to_crs(epsg=datum)
         
     def save_shp(self):
         self.shape_flname_mancha = QtWidgets.QFileDialog.getSaveFileName(None, "Selec", "", "shp files (*.shp)")
         flname = self.shape_flname_mancha[0]
         surfaces_to_kml(self.surf_surface, self.surf_water, flname)
         print('Shape file!')
 ​
     def save_report(self):
         nomes = ['seção {}'.format(i) for i in range(21)]
         alturas_a = [self.alturas[idx]-self.ct[idx] for idx in range(21)]
         data_array = np.array(
             [nomes,
             self.ds,
             self.qs,
             alturas_a]
         ).T
       
   
 ​
     def calcular(self):
         print('Calcul...')
         fc = 1
 ​
         qmax_barr = qmax_barragem(self.h, self.v)
         cotas(self.ponto_informado, self.srtm, self.h)
 ​
         c, dp, xs, ys = cotas_secoes(self.s, self.srtm)
         self.ct = [i[40] for i in c]
 ​
         self.alturas, self.qs = altura_de_agua_secoes(self.ds, dp, c, qmax_barr, self.v, self.h, fc)
 ​
         x_all = []
         y_all = []
         h_all = []
         for idx, vv in enumerate(self.alturas):
             for idx1 in range(len(xs[idx])):
                 h_all.append(vv)
                 x_all.append(xs[idx][idx1])
                 y_all.append(ys[idx][idx1])
 ​
         int1 = random.randrange(0,9,1)
         int2 = random.randrange(0,9,1)
         flname= 'srtm_cortado_'+str(int1)+str(int2)
 ​
         clip_raster(self.s, self.srtm, flname)
         clipado = rasterio.open(flname)
 ​
         xcoords, ycoords, z = get_coordinates(clipado)
 ​
         chull = convex_hull(self.s)
         mascara = check_if_is_inside(chull, xcoords, ycoords)
         xcoords, ycoords, z = xcoords[mascara], ycoords[mascara], z[mascara]
 ​
         pts = [[i, j, k] for i, j, k in zip(x_all, y_all, h_all)]
         pts = np.array(pts)
         pts_surf = [[i,j,k] for i, j, k in zip(xcoords, ycoords, z)]
         pts_surf = np.array(pts_surf)
        
         cloud1 = pv.PolyData(pts)
         cloud2 = pv.PolyData(pts_surf)
 ​
         self.surf_water = cloud1.delaunay_2d()
         self.surf_surface = cloud2.delaunay_2d()
         self.plotter.add_mesh(self.surf_water, name='water surface', color='blue')
         self.plotter.add_mesh(self.surf_surface, name='terrain surface', cmap='viridis', scalars=z)
         self.plotter.view_isometric()
 ​
 if __name__ == "__main__":
     import sys
     app = QtWidgets.QApplication(sys.argv)
     Dialog = QtWidgets.QDialog()
     ui = Ui_Dialog()
     ui.setupUi(Dialog)
     Dialog.show()
     sys.exit(app.exec_())