使用MATLAB-快速计算常用水文评价指标

%目标函数快速计算，用于查看 clear; clc; % 读取实测径流数据 Qo = load('observed.txt'); % 检查数据格式，假设数据是每行一个数值，如果数据是列向量则无需转置 if size(Qo, 1) == 1 % 如果Qo是行向量 Qo = Qo(:); % 转换为列向量 elseif size(Qo, 2) == 1 % 如果Qo已经是列向量 % Qo已经是列向量，无需转置 else error('实测径流数据文件格式不正确，应为每行一个数值。'); end % 读取模拟径流数据 Qs = load('simulated.txt'); % 检查数据格式，假设数据是每行一个数值，如果数据是列向量则无需转置 if size(Qs, 1) == 1 % 如果Qs是行向量 Qs = Qs(:); % 转换为列向量 elseif size(Qs, 2) == 1 % 如果Qs已经是列向量 % Qs已经是列向量，无需转置 else error('模拟径流数据文件格式不正确，应为每行一个数值。'); end % 调用函数计算NSE, KGE, R2 NSE = GetNSE(Qs, Qo); KGE = GetKGE(Qs, Qo); PBIAS=GetPBIAS(Qs, Qo); RMSE=GetRMSE(Qs, Qo); % 注意：这里需要修正R2计算函数中的变量名，使用a和b而不是Qs和Qo R2 = GetR2(Qs, Qo); % 应改为 R2 = calculateR2(Qs, Qo); % 显示结果 disp(['纳什效率系数（NSE）: ', num2str(NSE)]); disp(['Kling-Gupta效率系数（KGE）: ', num2str(KGE)]); disp(['决定系数（R²）: ', num2str(R2)]); disp(['PBIAS: ', num2str(PBIAS)]); disp(['RMSE: ', num2str(RMSE)]); % 将结果写入新的文本文件 %fileID = fopen('results.txt', 'w'); %if fileID == -1 % error('无法打开文件以写入结果。'); %end %fprintf(fileID, '纳什效率系数（NSE）: %f\n', NSE); %fprintf(fileID, 'Kling-Gupta效率系数（KGE）: %f\n', KGE); %fprintf(fileID, '决定系数（R²）: %f\n', R2); %fprintf(fileID, 'PBIAS: %f\n', PBIAS); %fclose(fileID); %% % 径流目标函数计算 clear; clc; % 读取实测径流数据 Qo = load('observed.txt'); % 检查数据格式，假设数据是每行一个数值，如果数据是列向量则无需转置 if size(Qo, 1) == 1 % 如果Qo是行向量 Qo = Qo(:); % 转换为列向量 elseif size(Qo, 2) == 1 % 如果Qo已经是列向量 % Qo已经是列向量，无需转置 else error('实测径流数据文件格式不正确，应为每行一个数值。'); end % 读取模拟径流数据 Qs = load('simulated.txt'); % 检查数据格式，假设数据是每行一个数值，如果数据是列向量则无需转置 if size(Qs, 1) == 1 % 如果Qs是行向量 Qs = Qs(:); % 转换为列向量 elseif size(Qs, 2) == 1 % 如果Qs已经是列向量 % Qs已经是列向量，无需转置 else error('模拟径流数据文件格式不正确，应为每行一个数值。'); end % 调用函数计算mNSE,PFAB %mNSE = GetmNSE(Qs, Qo); %PFAB = GetPFAB(Qs, Qo); PBIAS=GetPBIAS(Qs, Qo); %NRMSE=GetNRMSE(Qs, Qo); RMSE=GetRMSE(Qs, Qo); % 显示结果 %disp(['修改后的纳什效率系数（mNSE）: ', num2str(mNSE)]); %disp(['PFAB: ', num2str(PFAB)]); disp(['PBIAS: ', num2str(PBIAS)]); %disp(['NRMSE: ', num2str(NRMSE)]); disp(['RMSE: ', num2str(RMSE)]); % 将结果写入新的文本文件 fileID = fopen('results.txt', 'w'); if fileID == -1 error('无法打开文件以写入结果。'); end %fprintf(fileID, '修改后的纳什效率系数（mNSE）: %f\n', mNSE); %fprintf(fileID, 'PFAB: %f\n', PFAB); %fprintf(fileID, 'RMSE: %f\n', RMSE); fprintf(fileID, 'PBIAS: %f\n', PBIAS); fclose(fileID); % GetNSE函数定义（保持不变） function NSE = GetNSE(Qs, Qo) if length(Qs) == length(Qo) QoAve = mean(Qo); NSE = 1 - sum((Qs - Qo).^2) / sum((Qo - QoAve).^2); else error("实测径流和模拟径流长度不等"); end end % calculateKGE函数定义（保持不变） function KGE = GetKGE(Qs, Qo) QsAve = mean(Qs); QoAve = mean(Qo); COV = cov(Qs, Qo); CC = COV(1,2) / (std(Qs) * std(Qo)); BR = QsAve / QoAve; RV = (std(Qs) / QsAve) / (std(Qo) / QoAve); KGE = 1 - sqrt((CC-1)^2 + (BR-1)^2 + (RV-1)^2); end % calculateR2函数定义（修正变量名） function R2 = GetR2(Qs, Qo) % 修正为正确的输入变量名 if length(Qs) == length(Qo) QsAve = mean(Qs); QoAve = mean(Qo); tempX = sum( (Qo-QoAve).*(Qs-QsAve)); tempY = sqrt(sum( (Qo-QoAve).^2 ) * sum( (Qs-QsAve).^2)); R2 = (tempX / tempY).^2; else error("输入数据长度不等"); end end %function mNSE = GetmNSE(Qs, Qo) % if length(Qs) == length(Qo) % QoAve = mean(Qo); % mNSE = 1 - sum(abs(Qs - Qo)) / sum(abs(Qo - QoAve)); % else % error("实测径流和模拟径流长度不等"); % end %end %function PFAB = GetPFAB(Qs, Qo) % if length(Qs) == length(Qo) % PFAB = abs(sum(Qo - Qs) / sum(Qo)) * 100; % else % error("实测径流和模拟径流长度不等"); % end %end function PBIAS = GetPBIAS(Qs, Qo) PBIAS = sum(Qo - Qs) / sum(Qo) * 100; end %function NRMSE = GetNRMSE(Qs, Qo) % 计算均方误差 % mse = mean((Qo - Qs).^2); % 计算真实值的最大值和最小值 % y_range = max(Qo) - min(Qo); % 计算NRMSE % NRMSE = sqrt(mse)*100 / y_range; %end function RMSE = GetRMSE(Qs, Qo) % 计算预测值和真实值之间的误差平方和 squared_errors = (Qo - Qs).^2; % 计算平均误差平方 mean_squared_error = mean(squared_errors); % 计算均方根误差 RMSE = sqrt(mean_squared_error); end