qaa.png

Rrs.txt

0.0012	0.00169	0.00329	0.00404	0.00748	0.00346
0.00097	0.00119	0.00184	0.00229	0.00425	0.00161
0.00091	0.00102	0.00151	0.0019	0.0028	0.00179
0.00229	0.00346	0.00537	0.00649	0.00907	0.00638
0.00163	0.0017	0.00255	0.00309	0.00515	0.00268
0.0008	0.0008	0.00121	0.00147	0.00269	0.00122
0.00109	0.00136	0.0022	0.00268	0.00389	0.00182
0.00136	0.00213	0.00381	0.0047	0.00655	0.00282
0.0093	0.00799	0.00727	0.00444	0.00208	0.00017
0.0031	0.00247	0.00269	0.00234	0.0017	0.00013
0.00206	0.00158	0.00175	0.0016	0.0015	0.00029
0.00359	0.003	0.00341	0.00292	0.00235	0.0004

aw.txt

0.0031	0.0048	0.0132	0.0302	0.0595	0.439

bbw.txt

0.0034	0.0025	0.00158	0.00133	0.0009	0.00034

bands.txt

412	443	490	510	555	670

qaa_line.m

function [adg, aph] = qaa_line(Rrs, aw, bbw, bands)
%{
%一行
Rrs = [0.0012 0.00169 0.00329 0.00404 0.00748 0.00346];
aw = [0.0031 0.0048	0.0132 0.0302 0.0595 0.439];
bbw = [0.0034 0.0025 0.00158 0.00133 0.0009 0.00034];
bands = [412 443 490 510 555 670];
[adg, aph] = qaa_line(Rrs, aw, bbw, bands);

%多行
Rrs = load('Rrs.txt'); % 遥感反射率
aw = load('aw.txt'); % 纯水的吸收系数
bw = load('bbw.txt');  % 纯水的后向散射
bands = load('bands.txt'); % 波段波长数值
n = size(Rrs, 1);
adg = Rrs;
aph = Rrs;
for i = 1:n
    [adg(i, :), aph(i, :)] = qaa_line(Rrs(i, :), aw, bbw, bands);
end
%}

i_412 = 1;
i_443 = 2;
i_490 = 3;
i_55x = 5;
i_670 = 6;

%step 1，计算rrs
rrs = Rrs ./ (0.52+1.7*Rrs);

%step 2，计算u
g0 = 0.089;
g1 = 0.1245;
u = ( -g0 + sqrt(g0*g0 + 4*g1*rrs) ) / (2 * g1);

%step 3
Rrs_670 = Rrs(i_670);
rrs_443 = rrs(i_443);
rrs_490 = rrs(i_490);
rrs_55x = rrs(i_55x);
rrs_670 = rrs(i_670);
if Rrs_670 < 0.0015
    p1 = rrs_443 + rrs_490;
    p2 = rrs_55x + 5 * rrs_670 / rrs_490 * rrs_670;
    x = log10(p1/p2);  % 改为log10
    
    i_lamda = i_55x;
    h0 = -1.146;
    h1 = -1.366;
    h2 = -0.469;
    a_0 = aw(i_lamda) + 10^(h0+h1*x+h2*x*x);
else
    i_lamda = i_670;
    a_0 = aw(i_lamda) + 0.39*(rrs_670/(rrs_443+rrs_490))^1.14;
end

%step 4
bbp_0 = u(i_lamda)*a_0/(1-u(i_lamda))-bbw(i_lamda);

%step 5
e = exp( -0.9 * rrs_443 / rrs_55x );
g = 2.0*(1-1.2*e);

%step 6，计算bbp
bbp = bbp_0 * (bands(i_lamda)./bands).^g;

%step 7，计算a
a = (1-u) .* (bbw+bbp) ./ u;

%step 8
Zeta = 0.74+0.2/(0.8+rrs_443/rrs_55x);
S = 0.015+0.002/(0.6+rrs_443/rrs_55x);
Xi = exp(S*(442.5-415.5));

%step 9
a_412 = a(i_412);
a_443 = a(i_443);
aw_412 = aw(i_412);
aw_443 = aw(i_443);
adg_443 = ((a_412-Zeta*a_443)-(aw_412-Zeta*aw_443)) / (Xi-Zeta);
adg = adg_443*exp(-S*(bands-443));
aph = a - adg - aw;