1.8 湍流谱代码解读

Kaimal 模型

代码

import numpy as np
from scipy import signal
from scipy.optimize import curve_fit
from scipy.signal import welch
import matplotlib.pyplot as plt


def KaiFit(u, meanU, dt, TI_min, Nsegms, OveLapPer, component='u', averageing_method='mean', fig=False):
    """

    Args:
        u:
        meanU:
        dt:
        TI_min:
        Nsegms:
        OveLapPer:
        component:
        averageing_method:
        fig:

    Returns:

    """
    # u: velocity vector(time-series) of the subgroup
    # meanU: mean value of vector u
    # TI_min: Time Interval of the subgroup
    # Nsegms: Number of Welch's segments
    # OveLapPer: percent of Welch's Overlaping, 0.5 means 50%
    N = np.size(u)
    fs = 1/dt
    # Length of each segments
    Segms_L = round(TI_min*60*fs/Nsegms)
    # Length of overlapping interval
    OveLap_L = round(TI_min*60*fs*OveLapPer/Nsegms)
    u1 = signal.detrend(u, type='linear')
    # Transformed original frequency and PSD
    f0, Su0 = welch(u1, fs, nfft=N, nperseg=Segms_L, noverlap=OveLap_L, detrend=False)
    # Nomalized PSD
    Su0 = Su0 / np.trapz(Su0, f0)

    Nfreq = 60
    
    index0 = f0 > 0
    index1 = f0 <= 1
    index = index0.tolist() and index1.tolist()
    # Original frequency & PSD being cut
    Su0 = Su0[index][1:]
    f0 = f0[index][1:]
        
    f3 = np.logspace(np.log10(f0[0]), np.log10(f0[-1]), num=Nfreq, endpoint=True)
    Bin = np.digitize(f0, f3, right=True)
    bin_no = np.arange(0, Nfreq)
    
    # Representative frequency & PSD by sub-interval mean or median
    Su = []
    f = []
    stdY = []
    for i in bin_no:
        if averageing_method == 'mean': 
            Su.append(np.mean(Su0[Bin == i]))
            f.append(np.mean(f0[Bin == i]))
        elif averageing_method == 'median':
            Su.append(np.median(Su0[Bin == i]))
            f.append(np.median(f0[Bin == i]))
            stdY.append(np.std(Su0[Bin == i]))
    
    index = ~np.isnan(Su)
    # Representative frequency & PSD 
    Su_res = np.array(Su)[index.tolist()]
    f_res = np.array(f)[index.tolist()]
    # Su_res = Su0
    # f_res = f0
    if component == 'u' :
        def Kaifunc(f,L):
            fr = L * f /meanU
            S = 4. * fr / (1.+6.*fr)**(5./3.)
            S = S / f
            varU = np.trapz(S, f)
            S = S / varU  # PSD: Su/var(u)
            return S
    elif component == 'v':
        def Kaifunc(f, L):
            fr = L * f / meanU
            S = 4. * fr / (1.+6.*fr)**(5./3.)
            S = S /f 
            varU = np.trapz(S, f)
            S = S / varU  # PSD: Sv/var(v)
            return S
    elif component =='w':
        def Kaifunc(f, L):
            fr = L * f / meanU
            S = 4. * fr / (1.+6.*fr)**(5./3.)
            S = S / f
            varU = np.trapz(S, f)
            S = S / varU  # PSD: Sw/var(w)
            return S
        
    # PSD Fitting & Lengthscale estimating
    popt, pcov = curve_fit(Kaifunc, f_res, Su_res, p0=50)
    L = popt
    # f_new =  np.logspace(np.log10(f_res[0]), np.log10(f_res[-1]), num=Nfreq, endpoint=True)
    # Real frequencies in Hz
    f_new = f_res
    # Fitted PSD
    Su_new = np.array(Kaifunc(f_new, popt))

    # Plot [Su/var(u)] on [f in Hz]
    if fig == True:
        plot1 = plt.loglog(f_res, Su_res, 'o', markerfacecolor='none',color ='b',label='original values')
        plot2 = plt.loglog(f_new, Su_new, 'r', label='fit values_Kai')
        plt.xlabel('f [Hz]')
        if component == 'u':
            plt.title('Normlized Su_Kaimal_CF')
            plt.ylabel('Su/var(u)')
        elif component == 'v':
            plt.title('Normlized Sv_Kaimal_CF')
            plt.ylabel('Sv/var(v)')
        elif component == 'w':
            plt.title('Normlized Sw_Kaimal_CF')
            plt.ylabel('Sw/var(w)')       
        plt.legend()
        plt.show()
    # Output [Su/var(u)] on [f in Hz]
    return L, f_res, Su_res, Su_new

解读

这段代码是用于处理和分析风速数据的，主要用于计算和拟合功率谱密度（Power Spectral Density, PSD）。这个函数 KaiFit 的作用是对风速时间序列进行分析，计算其功率谱密度，并基于Kaimal谱模型进行拟合。下面我将解释代码中的各个部分：

1. 参数说明:

   • u: 风速向量（时间序列）。
   • meanU: 风速向量 u 的平均值。
   • dt: 时间间隔。
   • TI_min: 子组的时间间隔。
   • Nsegms: Welch方法的段数。
   • OveLapPer: Welch重叠百分比，例如0.5表示50%。
   • component: 表示风速的成分（‘u’, ‘v’, 或 ‘w’）。
   • averageing_method: 平均方法，可以是’mean’或’median’。
   • fig: 是否绘制图形。

2. 处理步骤:

   • 计算样本数、采样频率、每个段的长度以及重叠间隔的长度。
   • 对风速向量进行去趋势处理。
   • 使用Welch方法计算原始频率和功率谱密度（PSD）。
   • 将PSD归一化。
   • 对频率进行筛选和划分。
   • 根据平均方法计算每个频率段的代表频率和PSD。
   • 定义Kaimal谱模型函数 Kaifunc，用于拟合PSD。
   • 使用曲线拟合工具 curve_fit 对数据进行拟合，估算长度尺度。
   • 根据需要绘制图形。

3. 输出:

   • 返回长度尺度、频率和原始/拟合后的PSD值。

这个函数对于分析环境科学或工程领域中的风速数据很有用，尤其是在需要根据实际数据拟合风速模型时。它能够提供关于风速变化模式的详细信息，这对于例如风能评估等应用场合非常重要。