Skip to content

Commit 2475f7d

Browse files
HolgerPollyNetHolgerPollyNet
authored
Merge pull request #295 from PollyNET/GHK2DEV
With this merge, polarization calibration and calculation is performed following latest ACTRIS standards, i.e. using the GHK nomenclature.
2 parents 8158162 + c78eed4 commit 2475f7d

19 files changed

Lines changed: 1945 additions & 573 deletions

‎config/pollyDefaults/template_polly_defaults.json‎

Lines changed: 13 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -3,14 +3,24 @@
33
"depolCaliConstStd532": 0.0,
44
"depolCaliConst355": 0.1201,
55
"depolCaliConstStd355": 0.0,
6+
"polCaliEta355": 33,
7+
"polCaliEtaStd355": 1,
8+
"polCaliEta532": 11,
9+
"polCaliEtaStd532": 1,
10+
"polCaliEta1064": 0.7,
11+
"polCaliEtaStd1064": 1,
612
"LC": [16893959541573.252000, 1, 1, 1, 28466210203696.203000, 1, 1, 67876166905915.734000, 1, 1, 1, 1],
713
"LCStd": [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
814
"overlapFile532": "template_polly_overlap_532.txt",
915
"overlapFile355": "template_polly_overlap_355.txt",
10-
"molDepol532": 0.0053,
16+
"molDepol532": 0.0044,
1117
"molDepolStd532": 0.001,
12-
"molDepol355": 0.0075,
18+
"molDepol355": 0.0077,
1319
"molDepolStd355": 0.0015,
20+
"molDepol1064": 0.0036,
1421
"wvconst": 11.1,
15-
"wvconstStd": 0
22+
"wvconstStd": 0,
23+
"volDepolerror355": 0.003,
24+
"volDepolerror532": 0.004,
25+
"volDepolerror1064": 0.005
1626
}

‎config/template_pollynet_processing_chain_config_links.xlsx‎

-22.5 KB
Binary file not shown.

‎lib/calibration/depolCaliGHK.m‎

Lines changed: 264 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,264 @@
1+
function [polCaliEta, polCaliEtaStd, polCaliStartTime, polCaliStopTime,cali_status, globalAttri ] = depolCaliGHK(signal_t, ...
2+
bg_t, signal_x, bg_x, time, polCaliPAngStartTime, ...
3+
polCaliPAngStopTime, polCaliNAngStartTime, ...
4+
polCaliNAngStopTime, K, caliHIndxRange, ...
5+
SNRmin, sigMax, rel_std_dplus, rel_std_dminus, segmentLen, smoothWin)
6+
% DEPOLCALI polarization calibration for PollyXT lidar system.
7+
%
8+
% USAGE:
9+
% [polCaliEta, polCaliEtaStd, depol_cal_time] = depolCaliGHK(signal_t,
10+
% bg_t, signal_x, bg_x, time, polCaliPAngStartTime,
11+
% polCaliPAngStopTime, polCaliNAngStartTime,
12+
% polCaliNAngStopTime, K, caliHIndxRange,
13+
% SNRmin, sigMax, rel_std_dplus, rel_std_dminus,
14+
% segmentLen, smoothWin, flagShowResults)
15+
%
16+
% INPUTS:
17+
% signal_t: matrix
18+
% background-removed photon count signal at total channel.
19+
% (nBins * nProfiles)
20+
% bg_t: matrix
21+
% background at total channel. (nBins * nProfiles)
22+
% signal_x: matrix
23+
% background-removed photon count signal at cross channel.
24+
% (nBins * nProfiles)
25+
% bg_x: matrix
26+
% background at cross channel. (nBins * nProfiles)
27+
% time: array
28+
% datenum array represents the measurement time of each profile.
29+
% polCaliPAngStartTime: array
30+
% datenum array represents the start time that the polarizer
31+
% rotates to the positive angle.
32+
% polCaliPAngStopTime: array
33+
% datenum array represents the stop time that the polarizer
34+
% rotates to the positive angle.
35+
% polCaliNAngStartTime array
36+
% datenum array represents the start time that the polarizer
37+
% rotates to the negative angle.
38+
% polCaliNAngStopTime: array
39+
% datenum array represents the end time that the polarizer
40+
% rotates to the negative angle.
41+
% K: float
42+
% K parameter from GHK to correct the calibration.
43+
% caliHIndxRange: 2-element array
44+
% range of height indexes at which the signal can be used for
45+
% polarization calibration.
46+
% SNRmin: array
47+
% minimum SNR for calibration.
48+
% sigMax: array
49+
% maximum signal that could be used in the calibration to prevent
50+
% pulse pileup effects. (Photon Count)
51+
% rel_std_dplus: float
52+
% maximum relative uncertainty of dplus that is allowed.
53+
% rel_std_dplus: float
54+
% maximum relative uncertainty of dminus that is allowed.
55+
% segmentLen: integer
56+
% segement length for testing the variability of the calibration results
57+
% to prevent of cloud contamintaion.
58+
% smoothWin: integer
59+
% width of the sliding window for smoothing the signal.
60+
% flagShowResults: logical
61+
% flag to control whether to save the intermediate results.
62+
%
63+
% OUTPUTS:
64+
% polCaliEta: array
65+
% eta from polarization calibration.
66+
% polCaliEtaStd: array
67+
% uncertainty of eta from polarizatrion calibration.
68+
% polCaliStartTime: array
69+
% start time for each successful calibration.
70+
% polCaliStopTime: array
71+
% stop time for each successful calibration.
72+
% cali_status: integer
73+
% 1 if calibration is successfull, 0 otherwise.
74+
% globalAttri: struct
75+
% all the information about the depol calibration.
76+
%
77+
% REFERENCE:
78+
% 1. Freudenthaler, V.: About the effects of polarising optics on lidar signals and the Delta90° calibration, Atmos. Meas. Tech., 9, 4181-4255, 10.5194/amt-9-4181-2016, 2016.
79+
%
80+
% HISTORY:
81+
% - 2018-07-25: First edition by Zhenping.
82+
% - 2019-06-08: If no depol cali, return empty array.
83+
% - 2019-09-06: Remove the part to replace the bins of low SNR with NaN, because it will lead to bias when doing smoothing.
84+
% - 2024-08-27: Transfrom to GHK parameters using only eta and not V* (polCaliFac) any more
85+
% - 2024-10-02: Add calibration status variable to output.
86+
%
87+
% .. Authors: - zhenping@tropos.de, haarig@tropos.de
88+
89+
%% parameters initialization
90+
polCaliEta = [];
91+
polCaliEtaStd = [];
92+
mean_dminus = [];
93+
mean_dplus = [];
94+
std_dminus = [];
95+
std_dplus = [];
96+
polCaliStartTime = [];
97+
polCaliStopTime = [];
98+
cali_status = 0;
99+
globalAttri = struct();
100+
globalAttri.sig_t_p = cell(0);
101+
globalAttri.sig_t_m = cell(0);
102+
globalAttri.sig_x_p = cell(0);
103+
globalAttri.sig_x_m = cell(0);
104+
globalAttri.caliHIndxRange = cell(0);
105+
globalAttri.indx_45m = cell(0);
106+
globalAttri.indx_45p = cell(0);
107+
globalAttri.dplus = cell(0);
108+
globalAttri.dminus = cell(0);
109+
globalAttri.segmentLen = cell(0);
110+
globalAttri.indx = cell(0);
111+
globalAttri.mean_dplus_tmp = cell(0);
112+
globalAttri.std_dplus_tmp = cell(0);
113+
globalAttri.mean_dminus_tmp = cell(0);
114+
globalAttri.std_dminus_tmp = cell(0);
115+
globalAttri.K = cell(0);
116+
globalAttri.segIndx = cell(0);
117+
globalAttri.caliTime = cell(0);
118+
119+
if isempty(signal_t) || isempty(signal_x)
120+
warning('No data for polarization calibration.');
121+
return;
122+
end
123+
124+
days = unique(fix(time));
125+
nDays = length(days);
126+
127+
for iDay = 1:nDays
128+
for iDepolCal = 1:length(polCaliNAngStartTime)
129+
indx_45p = find(time >= days(iDay) & time < (days(iDay) + 1) & ...
130+
time >= polCaliPAngStartTime(iDepolCal) & ...
131+
time <= polCaliPAngStopTime(iDepolCal));
132+
indx_45m = find(time >= days(iDay) & time < (days(iDay) + 1) & ...
133+
time >= polCaliNAngStartTime(iDepolCal) & ...
134+
time <= polCaliNAngStopTime(iDepolCal));
135+
136+
if (length(indx_45p) < 4) || (length(indx_45m) < 4)
137+
% if not enough depol cali profiles were found, break the loop
138+
break;
139+
end
140+
141+
thisCaliStartTime = min([polCaliPAngStartTime(iDepolCal), ...
142+
polCaliNAngStartTime(iDepolCal)]);
143+
thisCaliStopTime = max([polCaliPAngStopTime(iDepolCal), ...
144+
polCaliNAngStopTime(iDepolCal)]);
145+
146+
% neglect the first and last profile which could be unstable due to
147+
% the rotation of the polarizer
148+
indx_45m = indx_45m(2:end-1);
149+
indx_45p = indx_45p(2:end-1);
150+
151+
sig_t_p = nanmean(signal_t(:, indx_45p), 2);
152+
bg_t_p = nanmean(bg_t(:, indx_45p), 2);
153+
SNR_t_p = pollySNR(sig_t_p, bg_t_p);
154+
indxBad_t_p = (SNR_t_p <= SNRmin(1)) | (sig_t_p >= sigMax(1));
155+
156+
sig_t_m = nanmean(signal_t(:, indx_45m), 2);
157+
bg_t_m = nanmean(bg_t(:, indx_45m), 2);
158+
SNR_t_m = pollySNR(sig_t_m, bg_t_m);
159+
indxBad_t_m = (SNR_t_m <= SNRmin(2)) | (sig_t_m >= sigMax(2));
160+
161+
sig_x_p = nanmean(signal_x(:, indx_45p), 2);
162+
bg_x_p = nanmean(bg_x(:, indx_45p), 2);
163+
SNR_x_p = pollySNR(sig_x_p, bg_x_p);
164+
indxBad_x_p = (SNR_x_p <= SNRmin(3)) | (sig_x_p >= sigMax(3));
165+
166+
sig_x_m = nanmean(signal_x(:, indx_45m), 2);
167+
bg_x_m = nanmean(bg_x(:, indx_45m), 2);
168+
SNR_x_m = pollySNR(sig_x_m, bg_x_m);
169+
indxBad_x_m = (SNR_x_m <= SNRmin(4)) | (sig_x_m >= sigMax(4));
170+
171+
172+
dplus = smooth(sig_x_p, 'moving', smoothWin) ./ ...
173+
smooth(sig_t_p, 'moving', smoothWin);
174+
dminus = smooth(sig_x_m, 'moving', smoothWin) ./ ...
175+
smooth(sig_t_m, 'moving', smoothWin);
176+
dplus(isinf(dplus)) = NaN;
177+
dminus(isinf(dminus)) = NaN;
178+
dplus(indxBad_t_p | indxBad_x_p) = NaN;
179+
dminus(indxBad_t_m | indxBad_x_m) = NaN;
180+
dplus = dplus(caliHIndxRange(1):caliHIndxRange(2));
181+
dminus = dminus(caliHIndxRange(1):caliHIndxRange(2));
182+
183+
mean_dplus_tmp = [];
184+
std_dplus_tmp = [];
185+
mean_dminus_tmp = [];
186+
std_dminus_tmp = [];
187+
segIndx_tmp = [];
188+
% find the most stable region where the realtive std of the signal
189+
% is less than rel_std_dminus and rel_std_dplus
190+
for iReg = 1:(caliHIndxRange(2) - caliHIndxRange(1) - segmentLen)
191+
192+
if sum(~ isnan(dplus(iReg:(iReg + segmentLen)))) ...
193+
<= segmentLen/4 || ...
194+
sum(~ isnan(dminus(iReg:(iReg + segmentLen)))) <= segmentLen/4
195+
continue;
196+
end
197+
198+
this_mean_dplus = nanmean(dplus(iReg:(iReg + segmentLen)));
199+
this_std_dplus = nanstd(dplus(iReg:(iReg + segmentLen)));
200+
this_mean_dminus = nanmean(dminus(iReg:(iReg + segmentLen)));
201+
this_std_dminus = nanstd(dminus(iReg:(iReg + segmentLen)));
202+
203+
if abs(this_std_dminus / this_mean_dminus) <= rel_std_dminus && ...
204+
abs(this_std_dplus / this_mean_dplus) <= rel_std_dplus
205+
segIndx_tmp = cat(2, segIndx_tmp, iReg);
206+
mean_dplus_tmp = cat(2, mean_dplus_tmp, this_mean_dplus);
207+
mean_dminus_tmp = cat(2, mean_dminus_tmp, this_mean_dminus);
208+
std_dplus_tmp = cat(2, std_dplus_tmp, this_std_dplus);
209+
std_dminus_tmp = cat(2, std_dminus_tmp, this_std_dminus);
210+
end
211+
end
212+
213+
% if there is no stable calibration segments, start the next
214+
% calibration
215+
if isempty(mean_dplus_tmp)
216+
continue;
217+
end
218+
219+
% find the most stable calbiration region
220+
[~, segIndx] = min(sqrt((std_dplus_tmp./mean_dplus_tmp).^2 + ...
221+
(std_dminus_tmp./mean_dminus_tmp).^2));
222+
indx = segIndx_tmp(segIndx);
223+
polCaliStartTime = cat(2, polCaliStartTime, thisCaliStartTime);
224+
polCaliStopTime = cat(2, polCaliStopTime, thisCaliStopTime);
225+
mean_dplus = cat(2, mean_dplus, mean_dplus_tmp(segIndx));
226+
std_dplus = cat(2, std_dplus, std_dplus_tmp(segIndx));
227+
mean_dminus = cat(2, mean_dminus, mean_dminus_tmp(segIndx));
228+
std_dminus = cat(2, std_dminus, std_dminus_tmp(segIndx));
229+
230+
% save the intermediate results
231+
globalAttri.sig_t_p{end + 1} = sig_t_p;
232+
globalAttri.sig_t_m{end + 1} = sig_t_m;
233+
globalAttri.sig_x_p{end + 1} = sig_x_p;
234+
globalAttri.sig_x_m{end + 1} = sig_x_m;
235+
globalAttri.caliHIndxRange{end + 1} = caliHIndxRange;
236+
globalAttri.indx_45m{end + 1} = indx_45m;
237+
globalAttri.indx_45p{end + 1} = indx_45p;
238+
globalAttri.dplus{end + 1} = dplus;
239+
globalAttri.dminus{end + 1} = dminus;
240+
globalAttri.segmentLen{end + 1} = segmentLen;
241+
globalAttri.indx{end + 1} = indx;
242+
globalAttri.mean_dplus_tmp{end + 1} = mean_dplus_tmp;
243+
globalAttri.std_dplus_tmp{end + 1} = std_dplus_tmp;
244+
globalAttri.mean_dminus_tmp{end + 1} = mean_dminus_tmp;
245+
globalAttri.std_dminus_tmp{end + 1} = std_dminus_tmp;
246+
globalAttri.K{end + 1} = K;
247+
globalAttri.segIndx{end + 1} = segIndx;
248+
globalAttri.caliTime{end + 1} = mean([thisCaliStartTime, thisCaliStopTime]);
249+
250+
end
251+
end
252+
253+
if isempty(mean_dminus) || isempty(mean_dplus)
254+
print_msg('Plus or minus 45° calibration is missing.\n', 'flagTimestamp', true);
255+
cali_status =0; %Calibration was not successfull.
256+
return;
257+
end
258+
259+
% calculate the depol-calibration factor and unceratinty and correcting it
260+
% with K
261+
polCaliEta = 1/K .* sqrt(mean_dplus .* mean_dminus);
262+
polCaliEtaStd = 0.5 .* (mean_dplus .* std_dminus + mean_dminus .* std_dplus) ./ sqrt(mean_dplus .* mean_dminus);
263+
cali_status =1; % Calibration successfull.
264+
end

‎lib/calibration/pollyMDRGHK.m‎

Lines changed: 77 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,77 @@
1+
function [MDR, MDRStd, flagDeft] = pollyMDRGHK(sigT, bgT, sigC, bgC, flagT,flagC, eta, voldepol_error, minSNR, deftMDR, deftMDRStd, PollyConfig)
2+
% POLLYMDR etimate the molecular depolarization ratio according to the measurements at reference height.
3+
%
4+
% USAGE:
5+
% [MDR, MDRStd, flagDeft] = pollyMDRGHK(sigT, bgT, ...
6+
% sigC, bgC, flagT,flagC, eta, voldep_sys_uncertainty, ...
7+
% minSNR, deftMDR, deftMDRStd)
8+
%
9+
% INPUTS:
10+
% sigT: array
11+
% signal strength of the total channel at reference height.
12+
% [photon count]
13+
% bgT: array
14+
% background of the total channel at reference height. [photon count]
15+
% sigCross: array
16+
% signal strength of the cross channel at reference height.
17+
% [photon count]
18+
% bgCross: array
19+
% background of the cross channel at reference height. [photon count]
20+
% flagT:
21+
% flag the total channel for the respective wavelength.
22+
% flagC:
23+
% flag the cross channel for the respective wavelength.
24+
% eta: scalar
25+
% depolarzation calibration constant.
26+
% voldep_sys_uncertainty: scalar
27+
% systematic uncertainty of the volume depolarization ratio (in
28+
% future it should be given in the config file)
29+
% minSNR: float
30+
% the SNR constrain for the the signal strength at reference height.
31+
% Choose a strong constrain for ensuring a stable result,
32+
% like 50 or 100.
33+
% deftMDR: float
34+
% default molecular depolarization ratio.
35+
% deftMDRStd: float
36+
% default std of molecular depolarization ratio.
37+
% Polly.Config:
38+
% contains the GHK parameters for further calcualtions.
39+
%
40+
% OUTPUTS:
41+
% MDR: float
42+
% retrieved molecular depolarization ratio.
43+
% MDRStd: float
44+
% std of retrieved molecular depolarization ratio.
45+
% flagDeft: logical
46+
% flag to show whether using the default values. If true,
47+
% it means default MDR and MDRStd were used.
48+
%
49+
% HISTORY:
50+
% - 2021-05-31: first edition by Zhenping
51+
% - 2024-08-13: MH changed to GHK
52+
%
53+
% .. Authors: - zhenping@tropos.de, haarig@tropos.de
54+
55+
MDR = deftMDR;
56+
MDRStd = deftMDRStd;
57+
flagDeft = true;
58+
59+
snrTot = pollySNR(sum(sigT), sum(bgT));
60+
snrCro = pollySNR(sum(sigC), sum(bgC));
61+
62+
flagValidPointTot = (snrTot >= minSNR);
63+
flagValidPointCro = (snrCro >= minSNR);
64+
65+
if (~ flagValidPointCro) || (~ flagValidPointTot)
66+
fprintf(['Too noisy in the reference height to calculate ' ...
67+
'the molecular depol ratio.\n']);
68+
return;
69+
end
70+
71+
[MDR, MDRStd] = pollyVDRGHK(sum(sigT), sum(sigC), ...
72+
PollyConfig.G(flagT), PollyConfig.G(flagC), ...
73+
PollyConfig.H(flagT), PollyConfig.H(flagC), ...
74+
eta, voldepol_error(1), voldepol_error(2), voldepol_error(3), 1);
75+
flagDeft = false;
76+
77+
end

‎lib/calibration/pollyPolCali.m‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
function [polCaliEta, polCaliEtaStd, polCaliFac, polCaliFacStd, polCaliTime, polCaliAttri] = pollyPolCali(data, transRatio, varargin)
2-
% POLLYPOLCALI calibrate the PollyXT cross channels for 355 and 532 nm with ±45° method.
2+
% POLLYPOLCALI calibrate the PollyXT cross channels for 355, 532 and 1064 nm with Delta90° method.
33
%
44
% USAGE:
55
% [polCaliEta, polCaliEtaStd,polCaliFac, polCaliFacStd, polCaliTime, polCaliAttri] = pollyPolCali(data, transRatio)
@@ -12,7 +12,7 @@
1212
%
1313
% KEYWORDS:
1414
% wavelength: char
15-
% '355nm' or '532nm'.
15+
% '355nm' or '532nm' or '1064nm'
1616
% depolCaliMinBin: numeric
1717
% minimum search index for stable polarization calibration constants.
1818
% depolCaliMaxBin

0 commit comments

Comments
 (0)