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

亚图跨际

已于 2024-04-15 15:39:50 修改

阅读量1k

点赞数 15

分类专栏： Python 技术工程交叉知识文章标签： python 流体统计建模线性非线性数学

于 2024-04-15 15:38:00 首次发布

本文链接：https://blog.csdn.net/jiyotin/article/details/137783335

版权

Python 同时被 3 个专栏收录

359 篇文章 24 订阅

订阅专栏

交叉知识

88 篇文章 5 订阅

订阅专栏

技术工程

29 篇文章 0 订阅

订阅专栏

🎯要点

流体数据统计模型：
1. 🎯曲线拟合和回归模型预测水流：🖊建立流量水位回归模型，模型拟合，绘制结果，🖊以为水位相对湿度的函数，建立多线性回归模型，绘制三维结果。🖊以水位，水位和相对湿度乘积的函数，建立非线性回归模型，绘制三维结果。
2. 🎯平稳性、趋势、季节性和周期性数据序列分析：🖊时间戳的水流数据分解为趋势和季节性分量，使用增强迪基富勒测试水流数据平稳性，绘制二维趋势图，🖊建立水流自回归、自回归移动平均线、自回归综合移动平均线和简单指数平滑模型，使用增强迪基富勒测试测试平稳性，均方根误差量化预测误差，绘制二维趋势图评估有效性。
3. 🎯显着性测试两处水文流量差异方法：🖊单因素和双因素方差分析，🖊t-测试，F-测试，柯尔莫哥洛夫-斯米尔诺夫检验，曼惠特尼测试。
4. 🎯估计不确定性：🖊非特定分布区间下，估计可能水文下降的范围，🖊自举置信区间估计，生成可替换多重采样水文，🖊对称置信区间估计一维水流，🖊分位数非参数置信区间估计和预测水文范围，🖊分位数回归推理和绘图检测水文超出指定分位数的异常值，🖊蒙特卡罗不确定性传播模拟曼宁方程，估计流体动力明渠中的流速，绘制二维图示。
流体模型：
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_())