Incorporating Satellite Precipitation Estimates into a Radar-Gauge Multi-Sensor Precipitation Estimation Algorithm

Yuxiang He 1,2,* ID , Yu Zhang 3, Robert Kuligowski 4, Robert Cifelli 5 and David Kitzmiller 1 1 Office of Water Prediction (OWP), National Weather Service (NWS), NOAA, Silver Spring, MD 20910, USA; david.kitzmiller@noaa.gov 2 University Corporation for Atmospheric Research (UCAR), Boulder, CO 80307, USA 3 University of Texas at Arlington, Arlington, Texas 76019, USA; yu.zhang@uta.edu 4 NOAA/NESDIS/Center for Satellite Applications and Research, College Park, MD 20740, USA; Bob.Kuligowski@noaa.gov 5 Physical Sciences Division, NOAA/Earth System Research Laboratory, Boulder, CO 80305, USA; rob.cifelli@noaa.gov * Correspondence: Yuxiang.He@noaa.gov; Tel.: +1-301-683-3567


Introduction
Precipitation estimation is critical in hydrology and weather forecasting over a wide range of spatiotemporal scales.For example, in the National Oceanic and Atmospheric Administration (NOAA)/National Weather Service's (NWS) flash flood and river forecast operations, quantitative precipitation estimates (QPEs) are one of the major inputs to the models and serve as a basis for watches and warnings [1,2].At present, operational QPEs are produced at NWS River Forecast Centers (RFCs) based primarily on rain gauge reports and the observations of Weather Surveillance Radars-1988 Doppler (WSR-88D) through the Multi-sensor Precipitation Estimator (MPE) [3].MPE consists of an interactive interface and a suite of automated algorithms, including those for correcting mean-field and range-dependent errors in radar, and one for blending radar and point gauge data.
Over the past two decades, a number of satellite-based QPE (SQPE) algorithms have been developed for applications in many areas of hydrology.Satellite-based precipitation estimates have both advantages and disadvantages compared to ground-based estimates.On the negative side, SQPE is limited by overpass frequency for microwave-based retrievals and a lack of robust correlation between cloud top brightness temperature and surface rainfall for geostationary platforms.However, satellite platforms also have inherent advantages compared to ground-based sensors because their views are not subject to terrain-based blockages or discontinuities along political boundaries due to lack of data, instrumentation differences or data embargo.In particular, The Geostationary Operational Environmental Satellite-R Series (GOES-R) products will improve precipitation estimation quality utilizing Lightning Imaging Sensor (LIS) and/or the higher temporal resolution that can fill in gaps between microwave overpasses.As a result, SQPEs are particularly useful in areas with few ground gauges, under radar gaps, and along national borders.Some of the most widely used algorithms retrieve rain rate using infrared (IR) brightness temperature observations from geostationary platforms, which have short latency and wide spatial-temporal coverage, and more accurate estimates from microwave sensors onboard lower-orbit satellites.These include the Climate Prediction C-Morphing (CMORPH) [4,5], TRMM Multi-Sensor Precipitation Analysis (TMPA) and Global Precipitation Measurement Integrated Multi-satellite Retrievals for GPM (IMERG) [6,7]; and the Self-Calibrating Multivariate Precipitation Retrieval (SCAMPR) [8][9][10][11].Recent advances in precipitation retrieval and calibration algorithms have led to substantial improvements to the latency, resolution and accuracy of SQPEs.For example, the recently launched Geostationary Operational Environmental Satellite (GOES)-16 satellite tripled the number of spectral bands offered by the previous GOES, which has the potential to further enhance the detection of precipitation.These enhancements offer new opportunities of using SQPEs to address the quality issue of radar QPE arising from factors such as overshooting, beam widening, and partial beam blockage.
In order to address inadequacies in the coverage of radar and gauges, the MPE has the capability to display the SQPE as a separate field for qualitative comparison, or the SQPE field can be cut and pasted into any other MPE products through the Graphical User Interface.Also the satellite QPE based on the Hydro-Estimator (H-E) algorithm [12] has been introduced to NWS operations [13], which however was limited by SQPE quality and the uneven distribution of errors in radar QPE.Further enhancement has been sought.There are many new efforts and attempts to integrate radar mosaic, satellite, and rain gauge network [14,15], but these works focus on a few cases and it is not clear how the algorithms perform over a wide range of precipitation regimes and geographic areas.As forecasting missions require seamless precipitation maps, a new module was developed for the MPE to bias-correct SQPE and to use the results to fill the radar gaps.The resulting products have been used operationally in a few NWS RFCs, most notably the West Gulf RFC (WGRFC), whose service area covers the drainage of the Rio Grande and requires SQPE for filling the radar-gauge gaps across the Mexico-US border.Though this MPE satellite-radar-gauge (SRG) module has been found to be overall satisfactory by forecasters, it is limited by design for filling radar coverage gaps and unable to fully harness the complementary skill of SQPE within the areas of the effective radar coverage.In addition, simple gap-filling operations often create discontinuities along the borders of the effective radar coverage areas due to dissimilarities in the storm configurations as seen by satellite-borne sensors and ground-based radar.
In this work, we develop an enhanced SRG fusion module that goes beyond simple gap filling, and seeks to blend the radar, gauge and satellite QPE over regions within the effective radar coverage and over areas immediately outside the effective radar coverage.In this way, the algorithm more effectively takes advantage of the different sensor inputs when and where they are available.This module consists of a preprocessor that employs local bias correction to mitigate systematic biases in the SQPE, and an advanced double two-way co-Kriging mechanism to objectively blend radar, gauge and satellite using an RQI-like (RQI: Radar Quality Index) mechanism to determine weighting of satellite data over locations.This new module is written in object-oriented C++ programming language and can run in standalone mode.It can be easily adapted for real-time integration of radar, gauge and satellite observations.The remainder of the paper is structured as follows.Section 2 reviews the existing MPE SQPE module and describes the new module and its underlying mechanisms.Section 3 describes the evaluation methodology for assessing the merits of the enhanced module.Section 4 presents the experimental design and data.Section 5 discusses the findings and discussions.Section 6 is the summary and conclusions of this work.

Multi-Sensor Precipitation Estimator (MPE)
The NWS MPE was conceived in the early 2000s to provide a full range of capabilities for radar QPE correction and multi-sensor blending [1,2,16].The earlier versions of MPE ingested radar-only precipitation estimates, including the Digital Precipitation Array (DPA) generated by the WSR-88D Precipitation Processing System (PPS), rain gauge records, and externally defined radar climatology to create bias-corrected multi-radar and multi-sensor QPEs.The current version is able to ingest Multi-Radar Multi-Sensor (MRMS) [2] radar-only estimates, and incorporates basic quality-assurance procedures.The primary algorithms include (1) mean-field bias (MFB) correction to correct the errors caused by radar calibration and uncertain Z-R relationship, (2) local bias (LB) correction to correct spatially varying bias, which is caused by rainfall inhomogeneity, range-degradation, and vertical profile of reflectivity, and (3) multi-sensor blending using the objective analysis techniques of Single Optimal Estimation (SOE) or Double Optimal Estimation (DOE) [17][18][19][20] to account for intermittency in statistical modeling of the spatial variability of precipitation.MFB determines a uniform gauge-radar ratio using spatial average during each time step based on recursive estimation and exponential smoothing.The ratio was calculated through dividing the average rainfall estimates measured at all of the valid rain gauges by the average radar rainfall estimates at the corresponding radar pixels over those gauges.LB calculation generates one bias value for each grid bin within a given coverage area; this bias may be spatially inhomogeneous since it is a function of range, azimuth, rainfall type, and other location-dependent factors.In MPE-SOE/DOE, weighting factors are calculated based on a covariance vector and matrix based on distances among gauges and radar values in the influencing radius for each bin.Heavy weight is placed on the gauge value for bins immediately surrounding the gauge, while the weight of the radar values increases as the distance from the gauge increases.Areas not covered by gauge or radar are filled in by the estimates from nearby gauge and radar values using the objective analysis technique to merge them by calculating their optimal weights and influence radius.SOE uses ordinary Kriging to directly estimate rainfall amount; DOE uses indicator (which is a Kriging analysis performed on a binary-transformed sample population [21]) and ordinary Kriging to estimate probability of rainfall and rainfall amount under the condition that it is raining [19].The operational version of the MPE system is configured to run in the Advanced Weather Interactive Prediction System (AWIPS) of the NWS.
Since 2006, SQPE based on the GOES-based H-E has also been ingested into MPE to produce satellite-radar-gauge mosaicked products for routine operational hydrologic forecasting.This SQPE first undergoes local bias correction using the MPE local bias module.The bias-corrected SQPE are then used to fill the gaps in the bias-corrected multi-radar mosaic of hourly precipitation totals, and the result is blended with point gauge data using the MPE SOE algorithm.This blended product has been used in real-time operations at several NWS RFCs.

Enhanced SRG Fusion Module
In order to go beyond simple gap-filling and fully utilize the complementary strength of SQPE, we devised a new, double two-way blending module that is implemented in a standalone C++ framework so that it can be executed in a number of computing environments including the AWIPS.
The new system improves over the existing one in four key ways: (1) Pre-processes the data to correct SQPE biases; (2) introduces the Radar Range Index (RRI) to consider the radar data quality due to beam broadening and the change in beam height with distance from the radar;  This module has been designed to be easily integrated into the operational MPE system.Three primary input data sources (radar DPA, gauge, and satellite) are required; in the future, lighting data might be added to identify severe convective weather.The major functions and algorithms of the system are similar to the operational MPE except for the enhanced MFM.It includes (A) multi-radar mosaicking, (B) mean field bias correction, (C) local bias correction, and (D) double two-way multisensor blending.The enhanced MFM includes a procedure to identify whether radar echoes are inside or outside the area where the radar beam is below the melting layer.Three different radar mosaic algorithms were implemented for a common location observed from multiple radars mosaic, This module has been designed to be easily integrated into the operational MPE system.Three primary input data sources (radar DPA, gauge, and satellite) are required; in the future, lighting data might be added to identify severe convective weather.The major functions and algorithms of the system are similar to the operational MPE except for the enhanced MFM.It includes (A) multi-radar mosaicking, (B) mean field bias correction, (C) local bias correction, and (D) double two-way multi-sensor blending.The enhanced MFM includes a procedure to identify whether radar echoes are inside or outside the area where the radar beam is below the melting layer.Three different radar mosaic algorithms were implemented for a common location observed from multiple radars mosaic, they are average, maximum, and the lowest elevation of multi-radar rainfall intensities.In this investigation, the rain rate on the lowest radar beam was used, and the information of corresponding radar beam height was saved.There are several ways of identifying the melting layer, such as from atmospheric sounding, from vertical profile of reflectivity, and from numerical weather prediction model forecast and analysis.For simplicity, in this study, the climatological melting layer height on the warm season was used to identify below or above melting layer based on the beam height product.If the echoes are outside this area, it determines how far away they are from the area where the radar beam is below the melting layer and computes the radar range index (RRI).RRI can be considered as a simplified version of radar quality index (RQI) which integrates additional information (e.g., blockages, depth of the bright band layer) to determine radar quality [22].The reason we employ RRI is that the method can be applied to existing DPA products, for which RQI is not available.
After computing the RRI, the modified double optimal estimation (MDOE) is used to blend radar and satellite QPE.Before this step, local bias correction is used to correct satellite data (Figure 1E).MDOE, the local bias correction, and the RRI method are described in the following section.

Local Bias Correction of Satellite Quantitative Precipitation
Before blending the satellite rainfall with the radar-gauge MPE fusion products as described above, the satellite rainfall needs to be corrected with respect to in-situ gauge data.This is done using a local bias correction algorithm in [17].The correction problem is formulated as a space-time estimation procedure in the radius of influence of the gauge.The correction is computed as soon as the corresponding satellite rainfall QPEs and bin-averaged gauge rain rates become available (sometimes more than one rain gauge was located in a HRAP grid box, the average procedure was implemented in the offline-MPE), and the correction is calculated sub-optimally via exponential smoothing to minimize processing time.

RRI Calculation
Satellite rainfall QPEs are more spatially coherent and temporally consistent than radar QPEs, so they are potentially useful in the region of poor radar data quality.Generally, they are less accurate than radar-derived precipitation, and therefore just considered as a useful supplement and not a replacement for radar.The RRI was developed to indicate the quality of radar QPEs as a function of beam height, with values calculated for each radar bin ranging from 0 (bad quality) to 1 (good quality).The RRI is based on concepts from the Italian Regional Environmental Protection Agency (ARPA-SMR) [23], the German Aerospace Center (DLR) [22,24].It is based primarily on radar beam height and width, as shown in Figure 2. The first circle O 1 is the radar coverage and the second circle O 2 represents the area of influence of the satellite QPE; their area of intersection is used to determine the RRI.To calculate it is actually a "crop circle" problem.So, , and a = 1 . D is the distance between radar site 1 and satellite pixel location 2. R 1 and R 2 are the radar QPE coverage radius and satellite QPE influence radius, respectively.Figure 2b shows the resulting RRI values as a function of distance from the radar, where R 1 and R 2 are set at 65 and 12 Hydrologic Rainfall Analysis Project (HRAP) pixels (approximately 4 km each), respectively.This is because WSR-88D precipitation products are on the Hydrologic Rainfall Analysis Project (HRAP) coordinate system.The digital precipitation array product for each radar is on a 131 × 131 HRAP pixel grid mesh with the radar situated in the center.The radius of coverage is therefore 65 HRAP pixels.The default value of R 2 is 12 HRAP pixels is determined using the radar data influence radius currently used in radar-gauge merging.This value is set based on the conditional correlograms of rainfall amount averaged over eight directions (refer to Figure 2 in the paper of Seo and Smith [25]).Sensitivity analysis suggests that R 2 varies is within a narrow range surrounding the 12 HRAP pixels for a specific geographical region, and therefore a fixed, rather than a spatially and temporally varying R 2 is used.
Remote Sens. 2018, 10, 106 6 of 20 gauge merging.This value is set based on the conditional correlograms of rainfall amount averaged over eight directions (refer to Figure 2 in the paper of Seo and Smith [25]).Sensitivity analysis suggests that R2 varies is within a narrow range surrounding the 12 HRAP pixels for a specific geographical region, and therefore a fixed, rather than a spatially and temporally varying R2 is used.
(a) (b) Figure 2. Illustration of the RRI calculation.On the left, (a) shows radar and satellite relative location, O1 is the radar location, O2 is the location for which the rain rate will be estimated, R1 is the radar coverage radius, and R2 is the range of influence of satellite QPE, it is updated cut-off lag distance (HRAP) with pre-defined correlation coefficient.On the right, (b) shows an example of calculated radar range index as a function of distance from radar, the solid line is the Radar Range Index and the dashed line is the satellite supplemental index, the distance D between O1 and O2 on the left is same as the x-axis on the right.

Double Optimal Estimation (DOE) for Radar and Satellite
The problem of optimally estimating the spatial distribution of rainfall using satellite (S) and radar (R) rainfall may be stated as determining the rainfall  0 at a location  0 using the radar and satellite rainfall data within the estimation domain .Assuming isotropy in the correlation structure of rainfall, we define the estimation domain  to be a circle centered at the point of estimation  0 whose radius equals the spatial correlation scale of rainfall.Our goal is to estimate the expected amount of rainfall at this ungauged location  0 , given the radar estimates of rainfall directly over and/or near  0 ,  0 through   , and the satellite estimates near  0 ,  0 through   , [ 0 | 0 , … ,   ,  0 , … ,   ], where the upper-and lowercase letters signify random variables and their observed values, respectively.Mathematical notation we used followsSeo [18].
Simply following the original DOE in the work of Seo [18] would weigh the radar and satellite data equally within the domain A, since the weighting is based only on the distance between each radar or satellite observation and the location  0 .However, the errors in the satellite rain rates are generally higher than radar estimates.Furthermore, the quality of the radar QPE varies in space and time due to a number of factors; e.g., the range dependence [22].In order to generate an optimal merged precipitation field, a weighting coefficient is introduced to allocate the contribution for each dataset in domain ; in our study, it is determined by the radar range index (RRI) and indicated by   for radar and by   for satellite.Here   is a function of the distance of  0 from the radar site

Double Optimal Estimation (DOE) for Radar and Satellite
The problem of optimally estimating the spatial distribution of rainfall using satellite (S) and radar (R) rainfall may be stated as determining the rainfall Z 0 at a location u 0 using the radar and satellite rainfall data within the estimation domain A. Assuming isotropy in the correlation structure of rainfall, we define the estimation domain A to be a circle centered at the point of estimation u 0 whose radius equals the spatial correlation scale of rainfall.Our goal is to estimate the expected amount of rainfall at this ungauged location u 0 , given the radar estimates of rainfall directly over and/or near u 0 , z R0 through z Rn , and the satellite estimates near u 0 , z S0 through z Sm , E[Z 0 |Z R0 , . . ., Z Rn , Z S0 , . . ., Z Sm ], where the upper-and lowercase letters signify random variables and their observed values, respectively.Mathematical notation we used followsSeo [18].
Similar to radar-gauge estimation, the first step in DOE to fuse the radar and satellite data is to rewrite the conditional expectation and variance, E[Z 0 |z R , z S ] and Var[Z 0 |z R , z S ], as follows: where Pr[Z 0 > 0|z R , z S ] is the conditional probability, and the next step is then to estimate Interested readers are referred to Seo [18] for the full derivations.
Simply following the original DOE in the work of Seo [18] would weigh the radar and satellite data equally within the domain A, since the weighting is based only on the distance between each radar or satellite observation and the location u 0 .However, the errors in the satellite rain rates are generally higher than radar estimates.Furthermore, the quality of the radar QPE varies in space and time due to a number of factors; e.g., the range dependence [22].In order to generate an optimal merged precipitation field, a weighting coefficient is introduced to allocate the contribution for each dataset in domain A; in our study, it is determined by the radar range index (RRI) and indicated by σ R for radar and by σ S for satellite.Here σ R is a function of the distance of u 0 from the radar site and the radar effective sampling volume above u 0 , its value is determined by the RRI, i.e., σ R = RRI, and σ S + σ R = 1; the satellite contribution is 1 − σ R for each data point in the estimation domain A. If the radar beam is below the melting layer at u 0 , σ S = 0; if u 0 is far outside the region where the radar beam is below the melting layer, we assume radar data no contribution to this region, and errors in satellite QPE is constant everywhere, so σ S = 1.0.The procedure for calculating RRI was given in Section 2.4.
Since the weight coefficient does not impact the conditional probability estimation, the conditional probability Pr[Z 0 > 0|z R , z s ] is estimated via the indicator kriging estimator based on conditional-probability approximation of the conditional expectation of indicator random variable, which is similar to [1].However, since SQPEs have a very high false alarm ratio, in the radar-satellite overlay region the conditional probability Pr 1) and ( 2) is influenced by the new parameter RRI which is estimated using the following linear estimator: where, E[Z 0 |z 0 >0] is the expected positive rainfall estimate at u 0 , E[Z Ri |z 0 >0] is the expected radar rainfall at u i given that it rained at u 0 , and E[Z Sj |z 0 >0] is the expected satellite rainfall at u j given that it rained at u 0 .The optimal weights, Γ Ri and Γ Si are obtained by minimizing the conditional error variance ].The optimal weights not only are function of radar and satellite data, but also function of σ S , σ R , σ S and σ R are determined by RRI: where The values of E[Z 0 |z R , z S , z 0 >0] and Var[Z 0 |z R , z S , z 0 >0] are estimated as in [16], but have different weighting coefficients which are assigned to the radar and satellite dataset covariance matrix/vector according to the value of RRI in the domain A.

Evaluation Method
The hourly rainfall estimates on the HRAP grid are evaluated to understand the accuracy and limitations of the new SRG-MSF algorithm, as shown in the following sections.The surface rain gauges are considered to be the ground truth for the HRAP grid boxes in which they are located; if there are multiple gauges in one grid box, the values of all the gauges within the HRAP grid box are averaged.Several continuous verification and categorical metrics are used in this study, as described below.

Continuous Verification Statistics
The percentage bias (PBIAS) [9], mean absolute error (MAE) Root Mean Square Error (RMSE), and correlation coefficient (r) are defined as follows: where E i is the MPE estimate, O i is the gauge observation, and N is the number of pair numbers of matched gauge and MPE estimates.

Categorical Verification Statistics
Categorical verification statistics measure the correspondence between the estimated and observed occurrences of precipitation ≥0.254 mm/h (minimum value recorded by rain gauge is 0.01 inch corresponding to 0.254 mm/h).They can be derived from a 2 × 2 contingency table which is constructed using dichotomous (rain-no rain) values from MPE estimates and independent gauge observations [26].In this table, the entries in each "cell" represent the number of hits, misses, false alarms and correct rejections.From the table values can be computed for bias score (BIAS), Probability of Detection (POD), False Alarm Ratio (FAR), and Critical Success Index (CSI) as follows: BI AS = hits + f alse alarms hits + misses POD = hits hits + misses FAR = f alse alarms hits + f alse alarms CSI = hits hits + misses + f alse alarms For BIAS, the valid range is 0 to infinity, the perfect score is 1.For POD and CSI, valid ranges are 0 (bad) to 1 (good), for FAR, its valid ranges is from 0 (good) to 1.0 (bad).

Experimental Design and Data
The experimental region is focused around Texas, where the climate varies from semi-arid in the west to humid in the east, with fewer gauges available in the west than in the east.As shown in Figure 3, the primary source of hourly rain gauge data is the Hydro-meteorological Automated Data System (HADS), which was reprocessed by Kim, et al. [27].It is based on the historical archive of the real-time data collected and distributed by the HADS Program in the NWS Office of Hydrologic Development (OHD), now the Office of Water Prediction (OWP).Data loss and quality issues with HADS and possible biases with respect to neighboring COOP gauges have previously been investigated [27].The available gauge sites are shown in Figure 3 as yellow dots.In order to make the validation independent, the total available gauge sites (~800) were split into two groups-one for MPE bias correction, the other for validation.
The DPA products were operationally produced by WSR-88D and stored in the Level III data suite.The DPA are on the HRAP grid (~3.7 km grid space) and cover the area within 230 km of the radar site, in each individual DPA product, reflectivity was converted to rain rate using one of several Z-R relationships, the choice being made at the supporting Weather Forecast Office based on the weather situation [28].The locations of the 22 NEXRAD radars used in this study are indicated by red stars in Figure 3. SCaMPR is used in SRG-MSF as satellite QPE, the SCaMPR used in this study have been projected on the HRAP coordinate system.SCaMPR rain rates are derived from IR brightness temperatures and derived quantities after being calibrated against microwave rain rates [8,10,11] in an effort to capture the accuracy of microwave rain rates along with the rapid refresh of GOES data.The SCaMPR, HADS, and radar DPA were collected hourly, so hour by hour fusion procedure was performed in this study.
Remote Sens. 2018, 10, 106 9 of 20 weather situation [28].The locations of the 22 NEXRAD radars used in this study are indicated by red stars in Figure 3. SCaMPR is used in SRG-MSF as satellite QPE, the SCaMPR used in this study have been projected on the HRAP coordinate system.SCaMPR rain rates are derived from IR brightness temperatures and derived quantities after being calibrated against microwave rain rates [8,10,11] in an effort to capture the accuracy of microwave rain rates along with the rapid refresh of GOES data.The SCaMPR, HADS, and radar DPA were collected hourly, so hour by hour fusion procedure was performed in this study.

Results
Precipitation estimations over the study area from satellite only (SQPE), satellite-gauge local bias correction (LSQPE), radar-gauge multi-sensor fusion (RG-MSF), and the new satellite-radar-gauge multi-sensor fusion (SRG-MSF) are compared and analyzed in this section, starting with a case study and then with an analysis of the five years of warm seasons (April-October) from 2003 to 2007.

Case Study
The case study is a summertime heavy rain event on 27 June 2007 that was part of an unusually wet weather pattern that began in May and lasted into early August.During that time, a trough of low pressure in the middle and upper atmosphere over Texas was surrounded by ridges of high pressure to the east and west.Atmospheric disturbances tracking south from the Plains states interacted with a rich supply of moisture from the Gulf of Mexico.On 26 June, a complex of thunderstorms developed in the unstable atmosphere over central Texas.The storms entrained a very warm, moist flow of air off the Gulf of Mexico (low-level jet).The intersection of the jet with outflow in the upper atmosphere caused the storm complex to stall.Rainfall intensity increased as individual storms repeatedly moved over the same areas.The highest rain total measured during the event was 19.06 inches (484 mm) at the Hydro-met station near Marble Falls, Texas.Approximately 18 inches (457 mm) of rain fell in a six-hour period, exceeding the 500-year return frequency of central Texas.This case was chosen in order to illustrate errors rising from assumed Z-R relationships in estimating radar rain rates.

Results
Precipitation estimations over the study area from satellite only (SQPE), satellite-gauge local bias correction (LSQPE), radar-gauge multi-sensor fusion (RG-MSF), and the new satellite-radar-gauge multi-sensor fusion (SRG-MSF) are compared and analyzed in this section, starting with a case study and then with an analysis of the five years of warm seasons (April-October) from 2003 to 2007.

Case Study
The case study is a summertime heavy rain event on 27 June 2007 that was part of an unusually wet weather pattern that began in May and lasted into early August.During that time, a trough of low pressure in the middle and upper atmosphere over Texas was surrounded by ridges of high pressure to the east and west.Atmospheric disturbances tracking south from the Plains states interacted with a rich supply of moisture from the Gulf of Mexico.On 26 June, a complex of thunderstorms developed in the unstable atmosphere over central Texas.The storms entrained a very warm, moist flow of air off the Gulf of Mexico (low-level jet).The intersection of the jet with outflow in the upper atmosphere caused the storm complex to stall.Rainfall intensity increased as individual storms repeatedly moved over the same areas.The highest rain total measured during the event was 19.06 inches (484 mm) at the Hydro-met station near Marble Falls, Texas.Approximately 18 inches (457 mm) of rain fell in a six-hour period, exceeding the 500-year return frequency of central Texas.This case was chosen in order to illustrate errors rising from assumed Z-R relationships in estimating radar rain rates.
Figure 4 shows the daily accumulation during 00:00:00 to 23:59:59 on 27 June 2007 with the four different algorithms in the MPE software package, (a) SQPE, (b) LSQPE, (c) RG-MSF, and (d) SRG-MSF.The rain map of SQPE is very different from the radar QPE distribution in terms of resolution and distribution pattern (Figure 4a,b vs. Figure 4c,d) because satellite QPE is retrieved from cloud-top properties tuned using infrequent microwave-based precipitation estimates.Figure 4b shows the SQPE after the local bias correction, which has a greater similarity with RG-MSF than with the unadjusted SQPE with respect to the main precipitation core distribution.The locations of the mesoscale convective centers in the LSQPE in Figure 4b corresponded reasonably well with the maxima in the RG-MSF and SRG-MSF in Figure 4c,d, the heavy precipitation in Figure 4a is obvious overestimated; however, even after local bias correction the LSQPE still places too much precipitation in northeast Texas and the Texas panhandle compared to the radar QPE products in Figure 4c,d.The LSQPE product also does not contain the spatial detail of the RG-MSF (Figure 4d).However, the LSQPE did fill the gaps and/or extended rainfall region on the radar-gauge MSF in the western portion of Texas (left bottom region in Figure 4c,d) and improved the quality of the estimates where the radar beam was above the freezing level.In addition, the rainfall distribution smoothly extends from radar coverage to satellite-only region.The discontinuity in the LSQPE over western Texas that also carries over into MSF, it appears to be a problem, actually the artifacts around radar site KMAF appeared in the uncorrected satellite QPE map.After local bias correction, the discontinuity became more conspicuous because the rain gauge distribution is sparse around this region.In some hours there was no gauge record (missing records or no gauge site) during the event on the far western region.Therefore, satellite data quality did not improve much.More gauge observations are needed in this region to make the local bias correction more effective.4c,d.The LSQPE product also does not contain the spatial detail of the RG-MSF (Figure 4d).However, the LSQPE did fill the gaps and/or extended rainfall region on the radar-gauge MSF in the western portion of Texas (left bottom region in Figure 4c,d) and improved the quality of the estimates where the radar beam was above the freezing level.In addition, the rainfall distribution smoothly extends from radar coverage to satellite-only region.The discontinuity in the LSQPE over western Texas that also carries over into MSF, it appears to be a problem, actually the artifacts around radar site KMAF appeared in the uncorrected satellite QPE map.After local bias correction, the discontinuity became more conspicuous because the rain gauge distribution is sparse around this region.In some hours there was no gauge record (missing records or no gauge site) during the event on the far western region.Therefore, satellite data quality did not improve much.More gauge observations are needed in this region to make the local bias correction more effective.The impact of the various stages of multi-sensor merging can be evaluated objectively using partition biases into hit, false alarm, miss and correct null.Thus, an analysis of biases within the sample of hourly precipitation gauge observations and remote sensor estimates were carried out (Table 1).Biases were calculated for various criteria.These included the cases in which both the estimate and gauge were greater than 0.254 (hit cases), either was greater than 0.254 (false alarm and missed event cases), and where the estimate or gauge value were greater than zero (total cases).The hit and total bias are very similar in Table 1.For total biases (either gauge or MPE estimate was nonzero), the MPE estimates had biases between 0.9 and 1.1 except for the uncorrected satellite (SQPE).It apparently overestimates rain rate, which had a highest bias of 1.54.Applying the LB correction to the SQPE (for LSQPE), the systematical bias is corrected, the bias was reduced to 0.90, and also reduced the RMSE by 13% (Table 1 Total).While the RG-MSF is significantly better than the SQPE even after the local bias correction, blending the LSQPE with the RG-MSF still produces a slight improvement in correlation coefficient (CC) and reduces the RMSE, producing the best estimates for this case (ρ = 0.62 and RMSE = 4.99).
As might be expected, given their generally lower skill, the satellite estimates produced appreciably more missed events and false alarms than did the radar-gauge and radar-gauge-satellite estimates.The number of hits increases and the numbers of false alarm and misses decrease when radar data is blended in RG-MSF and SRG-MSF.Both SQPE and LSQPE had approximately 40% more false alarms and over 100% more missed events than did either RG-MSF or SRG-MSF.Also, SQPE and LSQPE mean precipitation were nearly 400% higher in their false alarm cases than were the RG-MSF and SRG-MSF.

Statistical Verification
In addition to the case study, multi-sensor estimates were derived throughout five consecutive warm seasons (April through October 2003-2007) and evaluated using the continuous and categorical verification metrics.The verification dataset was divided into two subsets of gauge reports: one is the area which is within the radar coverage and locate above the melting level (Region I) as well, and another is outside the radar coverage, but close enough for the SRG-MSF to still be influenced by radar (Region II).The continuous and categorical verification statistics for the 1-hour totals of the four products for Region I and Region II are compared in the following section respectively.
A total of about 392 rainfall gauges and 125,016 hourly reports are considered here, in which about 40,000 h were located in the region I, and 10,000 h located in the region II.It should be noted that the statistical sample contains sites and hours at which the verifying gauge or MPE products reported measurable precipitation and therefore False Alarm Ratio (FAR) is a fairly high due to the extended light rain effect which is lower than rain gauge minimum records (<0.254 mm/h).

Region I
Table 2 shows the results for Region I.As expected, satellite QPE performs worse than radar, though the local bias correction significantly improves every continuous statistical measure except correlation coefficient, which has no substantial change.The RG-MSF is significantly better than LSQPE, with the notable exception of percent bias.The SRG-MSF performs best of all, with a 10% reduction in RMSE from adding the satellite component.However, the correlation coefficient is essentially unchanged from the RG-MSF and a significant wet bias remains from the radar component even after the addition of satellite QPE.The categorical verification statistics (which focus on rain/no rain identification without any information about rainfall amount) are consistent with the continuous validation statistics, the RG-MSF and SRG-MSF have a significant wet bias, which appears to indicate that the radar estimates have a wet bias in areas where rainfall is correctly identified and the algorithm extended the light rain under the gauge measurable records.To look at the differences in more detail, the continuous and categorical metrics for each of the warm months (April-October) are shown in Figures 5 and 6 respectively.The month by month variability in the following comparison reflects the most dominant type of events sampled in each month (e.g., deep convective or stratiform precipitation).Figure 5a indicates that most of the wet bias of the SQPE occurs in August and October, and that in the other months it generally has less of a wet bias than the RG-MSF.The local bias correction reduces the SQPE percentage bias to near zero in every month except for October in those cases, the adjustment causes the LSQPE to be a little bit dry.This indicates that the local bias correction algorithm is very useful and robust for satellite QPE bias correction in most months.Meanwhile, the RG-MSF has a wet bias every month, with peaks in June and August (Figure 5a), possible reason of overestimating rain rate in radar is likely due to inappropriate Z-R for the types of storms that are sampled in this region.For all months, there is less bias in the SRG-MSF than for the RG-MSF (Figure 5a), which shows that the value in using satellite rainfall information is helpful to improve rain rates estimation in the new algorithm.Comparison of MAE and RMSE statistics for the four estimates (Figure 5b,c) reveals that LSQPE has the largest MAE and RMSE almost for all months, even after the local bias correction; however, the SRG-MSF still has a lower MAE and RMSE than the SG-MSF.Since the addition of satellite rain rates to the RG-MSF has a largely neutral impact on correlation coefficient (Figure 5d       Figure 6 contains a similar comparison of the categorical metrics.They show that the tendency of radar to overestimate the spatial coverage of rainfall is persistent throughout the warm season, but satellite bias shows variation with severe overestimation in May and severe underestimation in October (Figure 6a).While the satellite significantly under-detects rainfall relative to radar (Figure 6b), it generally does produce significantly more false alarms than radar (Figure 6c).However, the addition of satellite, while improving detection slightly (Figure 6b) also increases false alarms somewhat (Figure 6c), resulting in a generally neutral impact on CSI (Figure 6d and Table 3).Carefully reviewed SRG-MSF products and compared with RG-MSF, the increased false alarms and other categorical metrics underperformers due to light rain, which is found on the edge of rain regions, its influence over percentage bias is negligible, but the addition of satellites in Region I the wet bias in the radar data is corrected along with the improvement of percentage bias.Figure 6 contains a similar comparison of the categorical metrics.They show that the tendency of radar to overestimate the spatial coverage of rainfall is persistent throughout the warm season, but satellite bias shows variation with severe overestimation in May and severe underestimation in October (Figure 6a).While the satellite significantly under-detects rainfall relative to radar (Figure 6b), it generally does produce significantly more false alarms than radar (Figure 6c).However, the addition of satellite, while improving detection slightly (Figure 6b) also increases false alarms somewhat (Figure 6c), resulting in a generally neutral impact on CSI (Figure 6d and Table 3).Carefully reviewed SRG-MSF products and compared with RG-MSF, the increased false alarms and other categorical metrics underperformers due to light rain, which is found on the edge of rain regions, its influence over percentage bias is negligible, but the addition of satellites in Region I the wet bias in the radar data is corrected along with the improvement of percentage bias.Region II is the area outside the radar coverage area (so there is no RG-MSF), but close enough to radar coverage edge, so that the SRG-MSF algorithm still can merge satellite QPE with radar component which is within the radar influence range.Table 4 shows the continuous statistics of Region II corresponding to the time period in Table 2.As in Region I, the LB correction overcorrects the upward bias in the SQPE, and integrating the radar data exacerbates the dry bias, though it does reduce the MAE and RMSE.The correlation coefficient of the SQPE is degraded slightly by the LB correction but is identical to the SQPE value for the SRG-MSF.The categorical validation statistics in Table 5 show a downward rainfall detection bias (POD) for the SQPE that is similar to that within the radar coverage umbrella (Table 3), and as in Region I it is made slightly worse by the LB correction.However, in Region II the radar is unable to correct the negative bias to the degree that it does in Region I, which would be expected since the radar has less weight in the SRG-MSF in Region II.The slight drying by the LB correction reduces the POD of the satellite somewhat, and yet merging with radar produces the best POD.Since the FAR is essentially very close for all three fields, the improved POD leads to the CSI of the SRG-MSF being the best of the three fields.Figure 7 shows the continuous statistical metrics for Region II, analogous to Figure 5 for Region I.The general patterns of behavior are similar, with the notable exception that the LB correction more frequently turns a wet bias in the SPE into a dry bias in the LSPE (Figure 7a), and once again the weaker radar adjustment in the SRG-MSF in Region II is unable to address this (and in some months the dry bias becomes even worse when radar data are included).The LB generally reduces the MAE of the LSPE relative to the SPE (with the interesting exception of April), and the SRG-MSF consistently has the lowest MAE and RMSE of the three (Figure 7b,c).The LB also appears to degrade the correlation significantly in April but has a more neutral impact in other months, while after blending with radar, the SRG-MSF product generally has better correlation than LB does (Figure 7d and Table 4).
The categorical metrics in Figure 8 for Region II are analogous to Figure 6 for Region I, and they again show similar results.The LB makes a dry rain detection bias in the SQPE slightly worse, and the addition of radar data is less able to remove this bias in Region II than in Region I because it has relatively less influence in the RG-MSF.This is driven by changes in the POD, since the FAR is not affected significantly by the additional blending in most months.The result is that the CSI is still better for the SRG-MSF than for the LSQPE, but the improvement is smaller outside the radar coverage umbrella than within it.
Remote Sens. 2018, 10, 106 17 of 20 by the additional blending in most months.The result is that the CSI is still better for the SRG-MSF than for the LSQPE, but the improvement is smaller outside the radar coverage umbrella than within it. (a)

Discussion
Ground-based radar, an active microwave sensor, provides more accurate precipitation estimates than can infrared-based precipitation methods, since the radar detects some physical

Discussion
Ground-based radar, an active microwave sensor, provides more accurate precipitation estimates than can infrared-based precipitation methods, since the radar detects some physical properties of hydrometeors, rather than temperature properties of cloud tops.However, satellite platforms provide spatially continuous monitoring, while ground radar cannot.All of the results shown in the previous section clearly indicate the superiority of the radar QPE to the satellite QPE, but subject to the range limitations of the radar estimates.Finally, rain gauge observations are crucial to correcting both satellite and radar QPE from time-dependent biases.Our results have shown that in situations with less-than-perfect radar coverage, satellite QPE is a useful supplement to the overall merged QPE product.
The heavy rainfall event case study shown here illustrates the capability of the MPE algorithm to fuse multisensory precipitation.Although there was modest enhancement of skill through the incorporation of satellite QPE (SRG-MSF vs. RG-MSF) for this case, it apparent that satellite QPE field functioned to smoothly fill the gaps and improved correlation considerably within the radar effective coverage with no negative impact on radar/gauge QPE in areas with below-melting-layer radar coverage.
In an evaluation using warm-season data for a five-year period, results indicate that the enhanced multi-sensor module satellite-radar-gauge QPE improved RMSE, MAE and PB relative to radar-gauge only, even in area close to radar units.Thus, the satellite information proved useful not only for filling radar coverage gaps, but areas where radar detects only hydrometeors above the melting layer.For the region just outside effective radar coverage, the infusion of information from nearby radar QPE improved on the satellite QPE through blending the local bias corrected satellite QPE value and the distance-weighted value from the radar-gauge QPE analysis.

Summary and Conclusions
In this paper, we describe an enhanced module that objectively blends QPEs from gauge, radar and satellite on the basis of radar range, availability of radar data, and distance of gauge and available radar data to the location of interest.This MFM exploits the complementary accuracy of satellite, radar and gauge QPE, and expands the use of GOES-R satellite QPE within radar coverage where quality of radar product is considered low.The basic mechanism of the MFM is the DOE algorithm from the original MPE, which is a form of co-Kriging.The MFM not only utilizes satellite QPE within radar range, but also allows radar QPE to be used immediately outside of radar range so to improve the accuracy of blended product and mitigate discontinuities along the radar boundary due to differences in radar and satellite QPEs.
To evaluate and demonstrate the value of the new algorithm SRG-MSF, we compared the various rain rate products with HADS gauge data over Texas for a test case and for a 5-year warm-season longitudinal study.The major findings are as follows: (1) The satellite-radar-gauge SRG MSF improves over using satellite-gauge or radar-gauge alone, mainly because the estimates of radar and satellite are corrected consulting with gauge before merging them together, this procedure makes SRG-MSF less biased than the radar and more highly correlated with gauges.(2) Within the radar coverage region (Region I) the merging of satellite data with radar and gauges improves the data quality during the warm season, apparently by compensating for degradation of the radar data in regions where the beam is at or above the melting level (See, Figure 5).Outside the radar coverage region (Region II) nearby radar rain rates can be used via SOE and DOE to create a merged field that improves over the LB-adjusted satellite rain rates.(3) Introducing the parameter RRI to express the degradation of radar data quality with distance from the radar (beam height) provides an effective strategy for merging satellite and radar data in the framework of MPE.The smooth transition, as well as improvement in PB and other evaluation metrics between radar and satellite in the final product map reveal that the new fusing algorithm is reliable and robust.(4) Although the MPE and the enhanced satellite data fusion software package are reliable, the quality of the merged fields is directly influenced by the quality of each component dataset; and some errors in the individual component data set (particularly satellite) will propagate into the merged fields.Consequently, pre-processing the component data to reduce its errors (e.g., using local bias correction based on gauge values) is necessary to minimize the final product errors, e.g., this study demonstrates the value of LB correction of the satellite QPEs used in this framework.
It should be noted that some of the errors identified during the evaluation are because of the mismatch in spatial scale between the point gauge values and the areal averages represented by the multi-sensor products; i.e., the right tail of the distribution of the point values will tend to be higher than the tail of the distribution of the corresponding areal averages.Furthermore, since the MSF algorithm is trying to find the optimal weights for gauge and radar to blend them by minimizing the conditional error variance [18], the fields will tend to be smoothed and thus will tend to underestimate the small-scale heaviest rain rates.
The algorithm in this work is a procedure to blend different sources of information which represent the same physical quantity.It is based on the idea that the merging result of multiple source of rainfall rates can provide a better quantity of interest than any single source.This is particularly the case when the errors from the different methods are complementary, in the sense that radar rain rates tend to perform well where satellite rain rates are weak.Thus, the SRG-MSF algorithm has potential benefits to estimate rain rates in/above the radar melting layer where radar data quality is very low.Future work might integrate the lightning detection information to identify different precipitation types to improve heavy rain rate estimates in MPE to achieve more accurate rain-rate distribution.
(3) uses DOE to seamlessly blend the satellite and radar data; (4) identifies the areas of low radar data quality.The schematic diagram of this module is shown in Figure 1.It depicts the data flow and algorithms at each processing step.The left side of the diagram shows the structure of the operational MPE.The right side shows the new multi-sensor fusion module (MFM) and major steps involved in merging radar, satellite, and gauge data.Remote Sens. 2018, 10, 106 4 of 20 The new system improves over the existing one in four key ways: (1) Pre-processes the data to correct SQPE biases; (2) introduces the Radar Range Index (RRI) to consider the radar data quality due to beam broadening and the change in beam height with distance from the radar; (3) uses DOE to seamlessly blend the satellite and radar data; (4) identifies the areas of low radar data quality.The schematic diagram of this module is shown in Figure 1.It depicts the data flow and algorithms at each processing step.The left side of the diagram shows the structure of the operational MPE.The right side shows the new multi-sensor fusion module (MFM) and major steps involved in merging radar, satellite, and gauge data.

Figure 1 .
Figure 1.The new Multi-Sensor Fusion System, showing the existing MPE on the left, and new developed Multi-Sensor Fusion System on the right.Major features include the Radar Range Index (RRI) to quantify the radar data quality based on beam broadening and increase in beam height with distance from the radar, the Modified-Sensor Estimator which uses Double Optimal Estimation (DOE) to blend the satellite and radar data, and pre-processing to identify regions of low radar data quality region and to correct SQPE biases, here A: Radar mosaic (RM), B: Mean field bias correction (MB), C: Local bias correction (LB), D: Radar-gauge fusion (RG-MSF), E: Satellite local bias correction (LSQPE), H: Satellite-radar-gauge multi-sensor fusion (SRG-MSF).

Figure 1 .
Figure 1.The new Multi-Sensor Fusion System, showing the existing MPE on the left, and new developed Multi-Sensor Fusion System on the right.Major features include the Radar Range Index (RRI) to quantify the radar data quality based on beam broadening and increase in beam height with distance from the radar, the Modified-Sensor Estimator which uses Double Optimal Estimation (DOE) to blend the satellite and radar data, and pre-processing to identify regions of low radar data quality region and to correct SQPE biases, here A: Radar mosaic (RM), B: Mean field bias correction (MB), C: Local bias correction (LB), D: Radar-gauge fusion (RG-MSF), E: Satellite local bias correction (LSQPE), H: Satellite-radar-gauge multi-sensor fusion (SRG-MSF).

Figure 2 .
Figure2.Illustration of the RRI calculation.On the left, (a) shows radar and satellite relative location, O 1 is the radar location, O 2 is the location for which the rain rate will be estimated, R 1 is the radar coverage radius, and R 2 is the range of influence of satellite QPE, it is updated cut-off lag distance (HRAP) with pre-defined correlation coefficient.On the right, (b) shows an example of calculated radar range index as a function of distance from radar, the solid line is the Radar Range Index and the dashed line is the satellite supplemental index, the distance D between O 1 and O 2 on the left is same as the x-axis on the right.

Figure 3 .
Figure 3. Test area of offline MPE.The red stars are radar sites and the yellow points are HADS gauge locations.There are approximately 800 gauges in the study area.Half of them are used for MPE bias correction, the rest of them serve as independent validation gauge.

Figure 3 .
Figure 3. Test area of offline MPE.The red stars are radar sites and the yellow points are HADS gauge locations.There are approximately 800 gauges in the study area.Half of them are used for MPE bias correction, the rest of them serve as independent validation gauge.
Remote Sens. 2018, 10, 106 10 of 20 MSF.The rain map of SQPE is very different from the radar QPE distribution in terms of resolution and distribution pattern (Figure 4a,b vs. Figure 4c,d) because satellite QPE is retrieved from cloudtop properties tuned using infrequent microwave-based precipitation estimates.Figure 4b shows the SQPE after the local bias correction, which has a greater similarity with RG-MSF than with the unadjusted SQPE with respect to the main precipitation core distribution.The locations of the mesoscale convective centers in the LSQPE in Figure 4b corresponded reasonably well with the maxima in the RG-MSF and SRG-MSF in Figure 4c,d, the heavy precipitation in Figure 4a is obvious overestimated; however, even after local bias correction the LSQPE still places too much precipitation in northeast Texas and the Texas panhandle compared to the radar QPE products in Figure

Figure 5 .Figure 5 .
Figure 5. Monthly average continuous statistical metrics from four different algorithms in the MPE software package vs. gauges over Region I for April through October of 2003-2007: (a) Percentage bias, (b) mean absolute error, (c) root mean square error, and (d) correlation coefficient.

Figure 7 .Figure 8 .
Figure 7. Monthly average continuous statistical metrics from three different algorithms in the MPE software package vs. gauges over Region II for April through October of 2003-2007: (a) Percentage bias, (b) mean absolute error, (c) root mean square error, and (d) correlation coefficient.

Figure 8 .
Figure 8. Monthly average categorical statistical metrics from three different algorithms in the MPE software package vs. gauges over Region II for April through October of 2003-2007 (Bias score, POD, FAR, and CSI).
and Q SS are the n × n, n × m, and m × m conditional matrices, whose ij th entries are given by Cov[Z Ri , Z Rj |z 0 >0], Cov[Z Ri , Z Sj |z 0 > 0], and Cov[Z Si , Z Sj |z 0 > 0], respectively, and Q 0R and Q 0S are the 1 × n and 1 × m conditional covariance vectors, whose i th and j th entries are given by Cov[Z 0 , Z Ri |z 0 >0] and Cov[Z 0 , Z Sj |z 0 >0], respectively.The estimation variance is given by:

Table 1 .
Area average and bias for estimated hourly rainfall from 26 June 2017 to 28 June 2007 vs. the amounts of the independent gauges from different algorithm in the MPE software package, (a) SQPE, (b) LSQPE, (c) RG-MSF, and (d) SRG-MSF, last two rows are correlation coefficient and RMSE, where Q MPE and Q Gauge are hourly rain rate of offline-MPE products and independent gauge, respectively.

Table 2 .
Percentage bias (PB%), mean average error (MAE), root mean square error (RMSE), and correlation (r) between hourly rainfall from four different algorithms in the MPE software package vs. rain gauge for all of Region I during 2003 to 2007 from April to October.
and Table2.) it appears that the reduction in MAE and RMSE is entirely due to bias reduction.

Table 3 .
Same as Table2, but area bias (BIAS), Probability of Detection (POD), False Alarm Ration (FAR), and Critical Success Index (CSI) for rain/no rain discrimination (the cut off threshold is 0.254 mm/h).

Table 4 .
Same as Table 2, but for Region II.

Table 5 .
Same as Table3, but for Region II.