使用Perez模型计算墙面辐照度时除零错误及pvlib异常问题求解

一、DHI/DNI同时为零的除零错误解决 #

问题原因 #

Perez模型计算天空清晰度参数ε时，核心公式为 ε = DHI / (DNI + DHI)，当DNI和DHI同时为0时，分母为0触发除零错误。这类情况通常出现在夜间或完全阴天无太阳辐射的时段。

解决方案 #

1. 数据预处理：筛选无效时段

   • 结合太阳位置，当太阳天顶角>90°（太阳在地平线以下）时，直接将GHI、DNI、DHI设为0，无需计算辐照度
   • 对太阳在地平线以上但DNI+DHI=0的时段，添加极小值（如1e-8）到分母避免除零，或直接标记为0（此时无有效辐射）

2. 代码层面的处理
   在调用get_total_irradiance前，对辐射数据做清洗：

   # 筛选太阳在地平线以上的时段
   day_time = solar_position['zenith'] < 90
   # 对夜间时段直接设为0
   ghi.loc[~day_time] = 0
   dni.loc[~day_time] = 0
   dhi.loc[~day_time] = 0
   # 处理白天DNI+DHI=0的情况，添加极小值避免除零
   zero_sum_mask = (dni + dhi) == 0
   dni.loc[zero_sum_mask & day_time] = 1e-8
   dhi.loc[zero_sum_mask & day_time] = 1e-8
   

二、数据不连续与数值异常问题解决 #

可能原因 #

1. 气象数据本身存在异常值（如DNI为负、DHI>GHI等不符合物理逻辑的数值）
2. 太阳位置与辐射数据的时间戳不匹配
3. 未处理太阳在地平线以下的无效计算（此时Perez模型可能输出异常值）

解决方案 #

1. 校验并清洗气象数据
   添加数据校验逻辑，修正不符合物理规则的数值：

   # 确保DNI非负
   dni = dni.clip(lower=0)
   # 确保DHI不超过GHI（物理上DHI是GHI的漫射部分，不能大于总辐射）
   dhi = dhi.clip(upper=ghi, lower=0)
   # 确保GHI非负
   ghi = ghi.clip(lower=0)
   

2. 强制夜间辐照度为0
   太阳在地平线以下时，平面接收的总辐照度应为0，直接覆盖计算结果：

   results['poa_global'] = total_irradiance['poa_global']
   results.loc[~day_time, 'poa_global'] = 0
   

3. 确认时间戳匹配
   确保weather_data的索引与timestamps完全对齐，避免因时间错位导致的异常跳变。

修改后的完整代码 #

import pandas as pd
import pvlib
from pvlib.location import Location
from pvlib.irradiance import get_total_irradiance

# 多伦多地理位置
latitude = 43.7
longitude = -79.4
location = Location(latitude, longitude, 'America/Toronto')

# 全年逐小时时间序列
timestamps = pd.date_range(start='2023-01-01 00:00:00', 
                           end='2023-12-31 23:59:59', 
                           freq='H', tz='America/Toronto')

# 读取气象数据
weather_data = pd.read_csv('WeatherfileML_tor_1.0.csv')
assert len(weather_data) == 8760, "数据需包含全年8760小时的逐小时记录"
weather_data.index = timestamps

# 提取并清洗辐射数据
ghi = weather_data['GHI'].clip(lower=0)
dni = weather_data['DNI'].clip(lower=0)
dhi = weather_data['DHI'].clip(lower=0, upper=ghi)

# 计算太阳位置
solar_position = location.get_solarposition(weather_data.index)
day_time = solar_position['zenith'] < 90

# 处理夜间时段与白天零辐射情况
# 夜间直接设为0
ghi.loc[~day_time] = 0
dni.loc[~day_time] = 0
dhi.loc[~day_time] = 0
# 白天DNI+DHI=0时添加极小值避免除零
zero_sum_mask = (dni + dhi) == 0
dni.loc[zero_sum_mask & day_time] = 1e-8
dhi.loc[zero_sum_mask & day_time] = 1e-8

# 计算外太空辐射
dni_extra = pvlib.irradiance.get_extra_radiation(weather_data.index)

# 计算不同朝向墙面的总辐照度
surface_tilt = 90  # 垂直墙面
# 四个朝向对应的方位角：北=0，东=90，南=180，西=270
for azimuth, direction in zip([0, 90, 180, 270], ['north', 'east', 'south', 'west']):
    total_irradiance = get_total_irradiance(
        surface_tilt=surface_tilt,
        surface_azimuth=azimuth,
        solar_zenith=solar_position['zenith'],
        solar_azimuth=solar_position['azimuth'],
        dni=dni,
        ghi=ghi,
        dhi=dhi,
        dni_extra=dni_extra,
        model='perez'
    )
    weather_data[f'poa_global_{direction}'] = total_irradiance['poa_global']
    # 强制夜间辐照度为0
    weather_data.loc[~day_time, f'poa_global_{direction}'] = 0

# 保存结果
weather_data.to_csv('solar_radiation_toronto_2023.csv')
print(weather_data[['poa_global_north', 'poa_global_east', 'poa_global_south', 'poa_global_west']].head())

额外注意事项 #

• 若气象数据存在缺失值，建议用线性插值或相邻时段均值填充，避免异常跳变
• 可通过pvlib.irradiance.check_irradiance函数自动校验辐射数据的合理性
• 垂直墙面的辐照度计算需注意方位角的定义：pvlib中方位角以正北为0，顺时针递增（东=90，南=180，西=270）

内容的提问来源于stack exchange，提问作者Aynkaran Aymmugan