本文参考文献主要为章文波教授发表的论文《不同类型雨量资料估算降雨侵蚀力》，这也是目前很多中文论文计算降雨侵蚀力的主要参考文献之一。

大家在使用时，我认为应注意几个问题：（1）不同数据源可能带来的误差影响；（2）根据自己的研究与数据选择合适的公式；（3）基于日降雨量的公式相对最准确，但国内外还有很多人使用不同的公式计算降雨侵蚀力。

这里主要介绍两个公式，一个是基于月降水，一个是基于日降水。

1、基于月降水的降雨侵蚀力计算：

分别使用月平均雨量与逐月雨量，如下图所示，逐月雨量精度相对较高，推荐使用逐月雨量。 青藏高原数据中心等网站也提供了逐月降雨数据等。

2、使用日降水数据

基于站点的日降水数据计算降雨侵蚀力，精度超过基于月降水的降雨侵蚀力。

α与β的取值可以用参数估计的方法，值分别为1.1515与1.6273；另外也可以根据公式8与9对参数α与β进行拟合。公式中，每年每个站点的参数并不完全相同，α与β为多年多个站点的平均值。

------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

一些代码，希望大家指正：

首先将TXT的站点数据通过代码整理到了excel，

path = 'E:\R_Factor\A_Pre_Txt\';
maindir = dir('E:\R_Factor\A_Pre_Txt\*.txt');for i =1:1:21for j=1:1:12inname = (i-1)*12+j;pathori = strcat(path,maindir(inname).name);%2010年1月的要单独读取if inname == 121importdata = importtxt201001(pathori);elseimportdata = import_txt(pathori);%读取txt文件，只读取前十列endif j == 1outdata = importdata;%第一个月的情况下elseoutdata =  [outdata;importdata];%使用分号纵向合并，后面11个月追加到一个月的矩阵里endend   %输出结果outname = string(1999+i);%设置输出年份名outpath = strcat('E:\R_Factor\B_PRE_CSV\\',outname,'.csv');csvwrite(outpath,outdata); disp(outpath);
end

下面这个是使用matlab自动生成的，可能不同的版本啥的不太一样

function import_txt = import_txt(filename, startRow, endRow)
%IMPORTFILE 将文本文件中的数值数据作为矩阵导入。
%   SURFCLICHNMULDAYPRE13011200001 = IMPORTFILE(FILENAME) 读取文本文件 FILENAME
%   中默认选定范围的数据。
%
%   SURFCLICHNMULDAYPRE13011200001 = IMPORTFILE(FILENAME, STARTROW, ENDROW)
%   读取文本文件 FILENAME 的 STARTROW 行到 ENDROW 行中的数据。
%
% Example:
%   SURFCLICHNMULDAYPRE13011200001 = importfile('SURF_CLI_CHN_MUL_DAY-PRE-13011-200001.TXT', 1, 25823);
%
%    另请参阅 TEXTSCAN。% 由 MATLAB 自动生成于 2022/10/28 18:56:52%% 初始化变量。
if nargin<=2startRow = 1;endRow = inf;
end%% 每个文本行的格式:
%   列1: 双精度值 (%f)
%	列2: 双精度值 (%f)
%   列3: 双精度值 (%f)
%	列4: 双精度值 (%f)
%   列5: 双精度值 (%f)
%	列6: 双精度值 (%f)
%   列7: 双精度值 (%f)
%	列8: 双精度值 (%f)
%   列9: 双精度值 (%f)
%	列10: 双精度值 (%f)
% 有关详细信息，请参阅 TEXTSCAN 文档。
formatSpec = '%5f%5f%6f%7f%5f%3f%3f%7f%7f%7f%[^\n\r]';%% 打开文本文件。
fileID = fopen(filename,'r');%% 根据格式读取数据列。
% 该调用基于生成此代码所用的文件的结构。如果其他文件出现错误，请尝试通过导入工具重新生成代码。
dataArray = textscan(fileID, formatSpec, endRow(1)-startRow(1)+1, 'Delimiter', '', 'WhiteSpace', '', 'TextType', 'string', 'EmptyValue', NaN, 'HeaderLines', startRow(1)-1, 'ReturnOnError', false, 'EndOfLine', '\r\n');
for block=2:length(startRow)frewind(fileID);dataArrayBlock = textscan(fileID, formatSpec, endRow(block)-startRow(block)+1, 'Delimiter', '', 'WhiteSpace', '', 'TextType', 'string', 'EmptyValue', NaN, 'HeaderLines', startRow(block)-1, 'ReturnOnError', false, 'EndOfLine', '\r\n');for col=1:length(dataArray)dataArray{col} = [dataArray{col};dataArrayBlock{col}];end
end%% 关闭文本文件。
fclose(fileID);%% 对无法导入的数据进行的后处理。
% 在导入过程中未应用无法导入的数据的规则，因此不包括后处理代码。要生成适用于无法导入的数据的代码，请在文件中选择无法导入的元胞，然后重新生成脚本。%% 创建输出变量
import_txt = [dataArray{1:end-1}];

接下来还要去掉异常：

clc;
clear;
path = 'E:\R_Factor\B_PRE_CSV\';
maindir = dir('E:\R_Factor\B_PRE_CSV\*.csv');for i = 20  %这个方法对2020年不适用，2020年没有一个站点正常，只能手动处理pathori = strcat(path,maindir(i).name);%导入数据%使用load函数导入数据，导入前四列（站点、经纬度和海拔）第六七列（月日）第十列（降雨量）alldata=load(pathori);%导入所有数据后，删除第五 第八九列alldata(:,9) = [];alldata(:,8) = [];alldata(:,5) = [];data=sortrows(alldata);%对数据进行排序，先按照第一列进行排序%删除data的第7列降雨量数据异常值，并且转为mm单位for k =1:length(data)if data(k,7)>30000data(k,7)=0;elsedata(k,7)=data(k,7)/10;endend%添加侵蚀性降水列,第8列for j=1:length(data)if data(j,7) >= 12data(j,8)=data(j,7);elsedata(j,8) = 0;endendzhandianhao = data(:,1);%读取站点号chuxiancishu = tabulate(zhandianhao);%统计站点出现的次数zhandianqingkuang = sortrows(chuxiancishu,2,'descend');%站点情况排序%计算平闰年,计算天数，天数应该等于站点出现的次数year = 1999+i;if(mod(year,4)==0 && mod(year,100)~=0 || mod(year,400)==0)daynumber = 366;elsedaynumber = 365;end%删除站点出现次数不等于天数的站点，这样很可能会少很多站点，但是比较快捷for k = 1:length(chuxiancishu)if chuxiancishu(k,2) ~= daynumber for m = 1:length(data)if data(m,1)==chuxiancishu(k,1)data(m,:) = 0;%将异常的行全部赋值为0endendend  end%删掉data中全部为0的行data(all(data==0,2),:)=[];%输出outname = string(1999+i);outpath=strcat('E:\R_Factor\C_PRE_removeabnormal\\',outname,'.xlsx');xlswrite(outpath,data); disp(outpath);%计算一下是否是整除365或366length(data)/365length(data)/366end    

然后计算R，但是好像是当时海拔不对，没有修改还说怎么回事来的，后来是修改了

clc;
clear;
path = 'E:\R_Factor\C_PRE_removeabnormal\';
maindir = dir('E:\R_Factor\C_PRE_removeabnormal\*.xlsx');%首先读取excel，判断平年与闰年，并计算站点数，保存到矩阵中
for i = 1:21pathori = strcat(path,maindir(i).name);%导入数据%使用xlsread函数导入数据data=xlsread(pathori);year = 1999+i;if(mod(year,4)==0 && mod(year,100)~=0 || mod(year,400)==0)daynumber(i,1) = 366;zhandianshu(i,1) = length(data)/366;elsedaynumber(i,1) = 365;zhandianshu(i,1) = length(data)/365;end
end%计算R因子值
for z = 1:21pathori = strcat(path,maindir(z).name);%导入数据%使用xlsread函数导入数据zong=xlsread(pathori);%使用站点数进行循环ZHDS = zhandianshu(z,1)-1;if daynumber(z,1) == 365 for i=0:1:ZHDS%这里的值是气象站点数减一nqsxjs(i+1,1)=0;%初始化年侵蚀性降水qsxjyts(i+1,1)=0;% 侵蚀性降雨天数%使用矩阵R存储年侵蚀R(i+1,1)= 0;%初始化矩阵for j=1:1:365if(zong(i*365+j,8)>0)%导入站点数据，命名为zong，365天。nqsxjs(i+1,1)=nqsxjs(i+1,1)+zong(i*365+j,8);qsxjyts(i+1,1)=qsxjyts(i+1,1)+1;endendif(nqsxjs(i+1,1)== 0)%该站点一年是否产生侵蚀性降雨，没产生等于0rpjjy(i+1,1)=0;bt(i+1,1)=0;aef(i+1,1)=0;R(i+1,1)=0;elserpjjy(i+1,1)=nqsxjs(i+1,1)/qsxjyts(i+1,1);%日平均侵蚀性降水量bt(i+1,1)=0.8363+18.177/rpjjy(i+1,1)+24.455/nqsxjs(i+1,1);%计算β与αaef(i+1,1)=21.586*(1/(bt(i+1,1)^7.1891));for m=1:1:365%计算年R因子R(i+1,1)=R(i+1,1)+(zong(i*365+m,8)^bt(i+1,1));    endR(i+1,1) = aef(i+1,1)* R(i+1,1);endend%为年度R因子值添加经纬度和海拔信息。for k=1:1:zhandianshu(z,1)zhandian(k,1)=zong(365*k,1);%站点名称weidufen(k,1)=zong(365*k,2);%提取纬度数据并转换为度a=floor(weidufen/100);%纬度的度right=(weidufen-100*a)/60;%纬度的分转换为0.xx度weidu=a+right;zhandian(k,2)=weidu(k,1);jingdufen(k,1)=zong(365*k,3);%提取经度数据并转换为度b=floor(jingdufen/100);%经度的度rightb=(jingdufen-100*b)/60;%经度的分转换为0.xx度jingdu=b+rightb;%xx.xx度格式zhandian(k,3)=jingdu(k,1);zhandian(k,4)=zong(365*k,4)/10;%海拔zhandian(k,5)=R(k,1);%年R值 月R值end   else%366天的情况for i=0:1:ZHDS %这里的值是气象站点数减一nqsxjs(i+1,1)=0;%初始化年侵蚀性降水qsxjyts(i+1,1)=0;% 侵蚀性降雨天数%使用矩阵R存储年侵蚀R(i+1,1)= 0;%初始化矩阵for j=1:1:366if(zong(i*366+j,8)>0) %导入站点数据，命名为zong，366天。nqsxjs(i+1,1)=nqsxjs(i+1,1)+zong(i*366+j,8);qsxjyts(i+1,1)=qsxjyts(i+1,1)+1;endendif(nqsxjs(i+1,1)== 0)%该站点一年是否产生侵蚀性降雨，没产生等于0rpjjy(i+1,1)=0;bt(i+1,1)=0;aef(i+1,1)=0;R(i+1,1)=0;elserpjjy(i+1,1)=nqsxjs(i+1,1)/qsxjyts(i+1,1);%日平均侵蚀性降水量bt(i+1,1)=0.8363+18.177/rpjjy(i+1,1)+24.455/nqsxjs(i+1,1);%计算β与αaef(i+1,1)=21.586*(1/(bt(i+1,1)^7.1891));for m=1:1:365%计算年R因子R(i+1,1)=R(i+1,1)+(zong(i*366+m,8)^bt(i+1,1));    endR(i+1,1) = aef(i+1,1)* R(i+1,1);endend%为年度R因子值添加经纬度和海拔信息。for k=1:1:zhandianshu(z,1)zhandian(k,1)=zong(366*k,1);%站点名称weidufen(k,1)=zong(366*k,2);%提取纬度数据并转换为度a=floor(weidufen/100);%纬度的度right=(weidufen-100*a)/60;%纬度的分转换为0.xx度weidu=a+right;zhandian(k,2)=weidu(k,1);jingdufen(k,1)=zong(366*k,3);%提取经度数据并转换为度b=floor(jingdufen/100);%经度的度rightb=(jingdufen-100*b)/60;%经度的分转换为0.xx度jingdu=b+rightb;%xx.xx度格式zhandian(k,3)=jingdu(k,1);zhandian(k,4)=zong(366*k,4)/10;%海拔zhandian(k,5)=R(k,1);%年R值 月R值end   end%输出outname = string(1999+z);outpath=strcat('E:\R_Factor\D_Rfactor\\',outname,'.xls');%输出表头biaotou=[{'sta','lon','lat','ele','Rfactor'};num2cell(zhandian)];xlswrite(outpath,biaotou); disp(outpath);  
end