電気学会論文誌 C（電子・情報・システム部門誌）

IEEJ Transactions on Electronics, Information and Systems
Vol.140 No.1 pp.40–48
DOI: 10.1541/ieejeiss.140.40

Paper

Three-Dimensional Data Compression and Fast High-Quality
Reconstruction for Phased Array Weather Radar
Ryosuke Kawami∗
Akira Hirabayashi∗a)
Hiroshi Kikuchi∗∗∗Student Member,
Non-member,
Member,Daichi Kitahara∗
Eiichi Yoshikawa∗∗
Tomoo Ushio∗∗∗∗Non-member
Non-member
Senior Member(Manuscript received April 11, 2019, revised Aug. 16, 2019)

This paper proposes a fast high-quality three-dimensional (3D) compressed sensing for a phased array weather radar(PAWR), which is capable of spatially and temporally high-resolution observation of the atmosphere. Because of thehigh-resolution, the PAWR generates huge observation data of approximately 500 megabytes every thirty seconds. To
transfer this huge data in a public internet line for real time weather forecast, an eﬃcient data compression technology
is required. The proposed method compresses the PAWR data by randomly transferring several measurements only inthe troposphere, and then reconstructs the missing measurements for each small 3D tensor data by minimizing a costfunction based on a prior knowledge on weather phenomena. The minimizer of the cost function can be quickly computed by using a convex optimization algorithm with Nesterov’s acceleration technique. Numerical simulations using
real PAWR data show the eﬀectiveness of the proposed method compared to conventional two-dimensional methods.
Keywords: Phased array weather radar, data compression, compressed sensing, convex optimization, Nesterov’s acceleration1.

pects of weather phenomena, and two bytes are used to describe each aspect. As a result, huge data of 490.9 megabytes
is generated every thirty seconds, and a data transfer without any compression technique needs a transfer rate of 130.9
megabits per second (Mbps). Thus the current PAWR system
exploits a private internet line with some cost, but it makes
real time weather forecast on the network diﬃcult.
To transfer the PAWR data in a public internet line, an efﬁcient compression technique is necessary as well as a highquality reconstruction algorithm. A standard public internet
line of 100 Mbps at most can only maintain 30 to 40 Mbps inaverage. By assuming that 35 Mbps is available, the requiredcompression rate becomes 26.73% approximately. Thus weset 25% compression rate as our goal in this paper.
To this end, we exploit random selection for compression
scheme as used in several existing papers (5)–(9) . In the PAWR
system, the observed data is processed as shown in the left
column of Fig. 1. The observed data are mixtures of the backscattered signals from various elevation angles within the antenna beamwidth. A digital beamforming technique (10)–(12) divides the mixed signals into the backscattered signals of theindividual elevations. Then, the demixed signals are used tocompute thirteen weather parameters including the reﬂection
intensity, whose data size is 37.8 megabytes for one rotationof the PAWR. The computational cost to obtain each parameter is not ignorable. In a scenario without compression, theentire 37.8 megabytes data are transferred to a remote server
and then analyzed. On the other hand, we exploit random selection after the beamforming as shown in the right columnof Fig. 1. Then, we can reduce not only the data size of thereﬂection intensity to 9.44 megabytes but also the cost for its
computation. Since the entire reﬂection intensity data are not

Introduction

Occurrence of extreme weather events is getting more frequent (1) . Among such events, there exist thunderstorms withheavy rain caused by quick growth of cumulonimbus clouds.
Thunderstorms only take tens of minutes from their generation to decay, while classical weather radars equipped with aparabolic antenna take ﬁve or ten minutes for sparse scan of
the whole sky. Thus the temporal and spatial resolution is not
suﬃcient for observation of these thunderstorms.
A phased array weather radar (PAWR) (2)–(4) was developedin order to overcome the above limitation. That developed in
Osaka University achieved the temporal resolution of thirty
seconds for scanning a hemisphere of a radius sixty kilometers. Furthermore, the spatial resolution along the range axis
was improved to one hundred meters, which means that the
number of sampling points in the range axis increased to 600.
Those in the azimuth and the elevation axes increased to 300
and 110, respectively. The measurements include thirteen asa) Correspondence to: Akira Hirabayashi. E-mail: akirahrb@
media.ritsumei.ac.jp∗
Graduate School of Information Science and Engineering,
Ritsumeikan University
1-1-1 Noji-higashi, Kusatsu, Shiga 525-8577, Japan∗∗
Aeronautical Technology Directorate, Japan Aerospace Exploration Agency
6-13-1 Osawa, Mitaka, Tokyo 181-0015, Japan∗∗∗
Center for Space Science and Radio Engineering, The University of Electro-Communications
1-5-1 Chofugaoka, Chofu, Tokyo 182-8585, Japan∗∗∗∗
Division of Electrical, Electronic and Infocommunications
Engineering, Osaka University
2-1 Yamadaoka, Suita, Osaka 565-0871, Japan© 2020 The Institute of Electrical Engineers of Japan.

40

3D Compression and Fast High-Quality Reconstruction for PAWR（Ryosuke Kawami et al.）
PAWR⇓
Observed mixture signals
⇓

PAWR
⇓
Observed mixture signals⇓

Digital beamforming

Digital beamforming⇓
Demixed reﬂected signals⇓⇓
Demixed reﬂected signals⇓⇓⇓

Random selection⇓

Full parameter computation

Reduced parameter computation⇓
Reﬂection intensity X (37.8 MB)⇓⇓
Reﬂection intensity Y (9.44 MB)⇓

Data transfer⇓

Data transfer⇓⇓⇓
Transferred intensity X⇓

Reconstruction⇓
Reconstructed intensity X∗⇓

Data analysis

Data analysis

the compressed tensor by minimizing a convex cost function
based on the following prior knowledge on the PAWR data.
First, all measurements outside of the troposphere are zeros
because there is no weather events. Second, the reconstructeddata should be consistent with the compressed one. Third, the
3D total variation is relatively small because of the local similarity in the spatial domain. Fourth, the ℓ1 -norm of the 3D
discrete cosine transform (DCT) coeﬃcients is also small because of the sparsity in the frequency domain. Since the costfunction is convex, the minimizer can be eﬃciently found byusing one of the convex optimization techniques (21)–(23) .
In the previous paper (9) , we further improved the convergence speed by Nesterov’s acceleration technique (24)–(27) , but
one issue that the cost function is not strongly convex† was
not considered. For a strongly convex cost function, the convergence speed of the alternating direction method of multipliers (ADMM) can be improved by applying Nesterov’s acceleration to dual variables (26) . However, since the proposedcost function is not strongly convex, the previous algorithm (9)cannot guarantee the convergence to the optimal solution. Although Nesterov’s acceleration with a restart rule is proposedto theoretically guarantee the convergence (26) , this strategy isslow because it takes much time, in every iteration, to judge
whether the restart should be done or not. Instead, in this paper, we apply the acceleration technique only in the ﬁrst several iterations, and after that the standard ADMM iterations
are used. Numerical simulations using real PAWR data showthat the proposed method outperforms the conventional 2D
reconstruction algorithms (5) (7) (8) with less computational time.
The rest of this paper is organized as follows. Section 2 formulates our PAWR data compression scheme. In Section 3.1,
we deﬁne a cost function for 3D reconstruction and propose
a blockwise reconstruction algorithm based on a convex optimization technique. Section 3.2 describes Nesterov’s acceleration for the reconstruction algorithm. Section 4 shows theeﬀectiveness of the proposed algorithm by numerical simulations using real data. Finally, Section 5 concludes this paper.

Fig. 1. Flowcharts of original and our data processing.

computed, we cannot exploit a Huﬀman code, even though it
can theoretically achieve the lowest compression rate.
The random selection compression scheme matches withreconstruction techniques in compressed sensing (CS) (13)–(15) .
CS was ﬁrst applied to radars for point targets (16)–(18) and thenapplied to weather radars (5)–(9) (19) (20) . For example, Mishra etal. used a low-rank matrix approximation to reconstruct twodimensional (2D) slices, in the entire three-dimensional (3D)
PAWR data, from compressed slices (5) . Shimamura et al. proposed a one-dimensional (1D) reconstruction algorithm that
promotes the sparsity in a wavelet domain (6) . In our previous
work, we proposed a two-dimensional (2D) reconstruction algorithm which promotes the local similarity and the sparsity
of 2D slice data in the spatial and the frequency domains, respectively (7) (8) . In this method, matrices are extracted as slice
data for ﬁxed elevations, in the same manner as the methodof Mishra et al. Such matrices are generated only after eachthirty seconds rotation of the PAWR. Hence the PAWR system has to transfer these matrices while saving new measurements in the next rotation. As a result, our previous technique
is not eﬃcient with respect to memory usage although the reconstruction quality is better than the other methods (5) (6) .
In this paper, not only to reduce the memory usage but also
to further improve the reconstruction quality, we propose a
3D reconstruction method. The PAWR sequentially captures
2D slice data for ﬁxed azimuth directions, and then creates
the entire 3D PAWR data by combining all 2D slice data. If
we apply our previous 2D reconstruction method (7) (8) to sucha slice in each azimuth direction, we can transfer the PAWR
data with low memory usage because we do not need to wait
the thirty seconds rotation. However, the reconstruction quality of this simple strategy is not good because the local similarity along the elevation axis is lower than that along the azimuth axis. To achieve both of low memory usage and highreconstruction accuracy, we have proposed a 3D reconstruction method (9) , and this paper proposes its improved version.
In the 3D reconstruction method (9) , we combine several contiguous slices into a 3D tensor and compress it. We randomly
select and transfer some measurements of the 3D tensor onlyin the troposphere. Then, we reconstruct the 3D tensor from

2.

Data Compression by Random Selection

In this paper, we focus on the reﬂection intensity among
thirteen weather parameters observed by the PAWR (2)–(4) . The
other parameters can be handled in similar ways. The reﬂection intensities are observed as a 3D tensor X ∈ RNR ×NA ×NE ,and we denote by x[nR , nA , nE ] ∈ R the (nR , nA , nE )-entry ofX, where 1 ≤ nR ≤ NR := 600, 1 ≤ nA ≤ NA := 300, and1 ≤ nE ≤ NE := 110, i.e., the number of measurements is
N := NR NA NE = 19,800,000.
The PAWR can observe the reﬂection intensities within ahemisphere of a radius sixty kilometers. This hemisphere includes the out of the troposphere and x[nR , nA , nE ] is 0 there.
We use this knowledge for both sampling and reconstruction.
Let Htro be the height of the troposphere, where the curvatureof the earth is not taken into account. In order to compress theoriginal data X, we randomly select M (< N) measurements
within regions in which the height is lower than Htro and thenconvert (N − M) unselected measurements into 0. This compression process is denoted by A : RNR ×NA ×NE → RNR ×NA ×NE ,† A function f : RN → R ∪ {∞} is called convex if f (λx + (1 − λ)y) ≤λ f (x) + (1 − λ) f (y) for all x, y ∈ RN and all λ ∈ (0, 1), and strongly convex
if f (λx + (1 − λ)y) < λ f (x) + (1 − λ) f (y) for all x, y ∈ RN and all λ ∈ (0, 1).

41

IEEJ Trans. EIS, Vol.140, No.1, 2020

3D Compression and Fast High-Quality Reconstruction for PAWR（Ryosuke Kawami et al.）constructed by solving a convex optimization problem
Xl∗ :=

argminXl ∈S s.t. Al (Xl )=Yl

TV(Xl ) + λ ∥C(Xl )∥1 . · · · · · · ·(1)NR ×s×NE
E
In (1), λ >C : RNR ×s×N
denotes the 3D DCT,
P0,
P NE → R
NR P s∥Xl ∥1 := nR =1 nA =1 nE =1 |xl [nR , nA , nE ]|, and the total variation is deﬁned byTV(Xl ) :=

NX
NE
s X
R −1 X|xl [nR + 1, nA , nE ] − xl [nR , nA , nE ]|nR =1 nA =1 nE =1(a) Original data X[ : , 1, : ] [dBZ]+(b) Compressed data Y[ : , 1, : ] [dBZ]

Fig. 2. Original data and compressed data.

|xl [nR , nA + 1, nE ] − xl [nR , nA , nE ]|nR =1 nA =1 nE =1+

the compressed data is expressed as Y := A(X), and the compression ratio is deﬁned by α := M/N ∈ (0, 1).
An example of a pair of the original observed data X andthe compressed data Y is shown in Fig. 2, where slices of the
1st azimuth direction are extracted from the entire 3D tensors.
The original data X was acquired on March 30th, 2014, by the
PAWR equipped in Suita Campus of Osaka University, Japan,
and the compressed data Y retains only 25% measurements
among X. Now, we must recover the missing measurementsfrom the selected remaining ones and their positions. The
next section proposes a high-quality reconstruction method.
3.

NR X
NE
s−1 X
X

NR X
s NX
E −1
X

|xl [nR , nA , nE + 1] − xl [nR , nA , nE ]|.

nR =1 nA =1 nE =1

In (1), the cost function is convex and two constraint sets arealso convex† . Therefore, this problem is a convex optimization problem and can be solved by the simultaneous direction
method of multipliers (SDMM) (22) , which is a special case of
ADMM (21) , as follows.
b b
Let xl := vec(Xl ) ∈ RN (N:= sNR NE ) be a vectorized version of a 3D tensor Xl ∈ RNR ×s×NE , and the inverse mappingb bof vec is denoted by ten, i.e., Xl = ten(xl ). Let Al ∈ RN×N be a
diagonal matrix which converts the unselected measurementsof xl into 0 while retaining the selected ones. Deﬁne
0 ... 0 −1 1
 0 −1 1 . . . 0 

(N −1)×NR,
DNR :=  .
..  ∈ R R.. .. ..
 ..
.
.
. . 

0 . . . 0 −1 1

Three-Dimensional Reconstruction Algorithm
Based on Blockwise Convex Optimization

The proposed method in this paper is an improved versionof our previous 3D reconstruction method (9) . In Section 3.1,
we summarize our previous method (9) , where the data ﬁdelity
term is modiﬁed from the cost into constraint. In Section 3.2,
we apply Nesterov’s acceleration while considering the nonstrong convexity which was ignored in the previous paper (9) .
3.1 Proposed 3D Blockwise Reconstruction
In thePAWR system, a 2D matrix data X[ : , nA , : ] ∈ RNR ×NE as inFig. 2(a) is acquired at the same time. By rotationally acquiring X[ : , nA , : ] (nA = 1, 2, . . . , NA ) in sequence and combining them, the entire 3D tensor data X ∈ RNR ×NA ×NE is generated. In this section, we propose to divide the entire tensor X
into small tensors and to reconstruct each tensor in sequence.
At ﬁrst, with some factor s of NA , we divide the original
data X and the compressed one Y into NA /s tensors Xl :=(X[ : , (l − 1)s + 1, : ], X[ : , (l − 1)s + 2, : ], . . . , X[ : , ls, : ]) ∈RNR ×s×NE and Yl = Al (Xl ) ∈ RNR ×s×NE (l = 1, 2, . . . , NA /s).
We suppose that each compression process Al : RNR ×s×NE →RNR ×s×NE has the same compression ratio M/N, i.e., we randomly select sM/NA measurements from each Xl within regions in which the height is lower than Htro .
To reconstruct the complete measurements Xl from the randomly selected ones Yl , we exploit four characteristics of Xl .
First, the measurements outside of the troposphere are 0 as
mentioned before. Let us denote the set of all Xl ∈ RNR ×s×NE
satisfying this condition by S. Second, the observation of Xl
should be consistent with Yl . Third, the reﬂection intensities
are locally similar because rain falling areas exist continuously. Fourth, we suppose that Xl can be sparsely describedin the frequency domain. On the basis of the above characteristics, we reconstruct X as X∗ := (X1∗ , X2∗ , . . . , XN∗ A /s ) ∈RNR ×NA ×NE , where each block Xl∗ (l = 1, 2, . . . , NA /s) is re-and also deﬁne D s ∈ R(s−1)×s and DNE ∈ R(NE −1)×NE in thesame manner. The identity matrix and the 1D DCT matrix ofsize K × K are denoted by IK and C K , respectively. By deﬁning matrices, with the use of the Kronecker product ⊗, as
L1 := INE ⊗ I s ⊗ INR ,L2 := INE ⊗ I s ⊗ DNR ,




L
:=
I
⊗
D
⊗
I
,
L4 := DNE ⊗ I s ⊗ INR ,

3
NE
s
NR



L := C ⊗ C ⊗ C , L := I ⊗ I ⊗ I ,
5
NE
s
NR
6
NE
s
NR
and deﬁning g1 (y1 ) := 0 if Al y1 = vec(Yl ), g1 (y1 ) := ∞ if
Al y1 , vec(Yl ), gi (yi ) := ∥yi ∥1 (i = 2, 3, 4), g5 (y5 ) := λ ∥y5 ∥1 ,g6 (y6 ) := 0 if ten(y6 ) ∈ S, and g6 (y6 ) := ∞ if ten(y6 ) < S, avectorized version of Xl∗ in (1) is expressed as
6X∗
gi (Li xl ).
xl = argmin
xl ∈RNb i=1· · · · · · · · · · · · · · · · · · · · · · · · (2)The vector x∗l in (2) can be computed by SDMM as shownin Algorithm 1, and the original tensor data is reconstructedPb b
by Xl∗ = ten(x∗l ). In Algorithm 1, Q := 6i=1 LiT Li ∈ RN×N
−1
is invertible and the analytic form of Q can be derived in asimilar way as in our previous method (7) (8) without any prodb × N.
b Moreover, for any γ >uct of huge matrices of size N0, each proximity operator proxγgi (wi ) := argminyi gi (yi ) +12(22) (23).
2γ ∥yi − wi ∥2 (i = 1, 2, . . . , 6) can be easily computed† Since the range resolution of the PAWR is 100 [m] and the elevation resπ[rad], the height corresponding to xl [nR , nA , nE ] is given byolution is 220nE π
100 nR sin( 220 ) [m]. Furthermore, the height of the troposphere Htro [m] is
ﬁxed, and hence S is expressed as S := {Xl ∈ RNR ×s×NE | xl [nR , nA , nE ] = 0
Eπfor all (nR , nA , nE ) s.t. 100 nR sin( n220
) > Htro }. This set is clearly convex.

42

IEEJ Trans. EIS, Vol.140, No.1, 2020

3D Compression and Fast High-Quality Reconstruction for PAWR（Ryosuke Kawami et al.）Algorithm 1 3D PAWR data reconstruction
Input: Yl ∈ RNR ×s×NE , γ ∈ (0, ∞), and an integer Imax
1: yi,0 ← Li vec(Yl ) and zi,0 ← 0(i = 1, 2, . . . , 6)
2: for k = 1, 2, . . . , Imax do
P3:
xk ← Q−1 6i=1 LiT (yi,k−1 − zi,k−1 )4:yi,k ← proxγgi (Li xk + zi,k−1 ) (i = 1, 2, . . . , 6)5:zi,k ← zi,k−1 + Li xk − yi,k (i = 1, 2, . . . , 6)
6: end for
Output: x∗l ← xk

4.

Numerical Simulations

To show the eﬀectiveness of the proposed method, we conducted simulations using real PAWR data shown in Fig. 3(a).
This data was acquired on March 30th, 2014, by the PAWR
equipped in Suita Campus of Osaka University, Japan, and isreferred as Data 1. The following simulations were executedby Matlab on iMac (OS 10.10, Intel Core i5, 2.7 GHz, 8 GB).
The compressed data Y was generated with compression ratio
α = 0.25, where the height parameter Htro was set to 15 [km].
Each block Xl was generated by combining s = 4 contiguous
slices. This means that the PAWR system only has to store
4 slices for sending each small block while the conventional
2D reconstruction methods (5) (7) (8) have to store the entire 300
slices corresponding to one rotation of the PAWR. After sending one block, memory can be overwritten by the succeeding
block. As a result, we can achieve 300/4 = 75 times lower
memory usage than those of the conventional methods (5) (7) (8) .
Algorithm 1 reconstructed the 3D tensors Xl from Yl withthe use of λ = 4 and γ = 1.7, and Algorithm 2 reconstructed
Xl with the use of λ = 4, γ1 = 2.2, γ2 = 1.3 and P = 5. These
values were set so that the reconstruction quality and the convergence rate would be the best. We compared the proposed
method with the 2D methods of Mishra et al. (5) and ours (7) (8) ,which are called Conventional 1 and Conventional 2. We also
compared the results of the proposed method for s = 300, i.e.,
the results of the non-blockwise reconstruction. Note that theshapes of the reconstructed blocks are diﬀerent depending on
the reconstruction methods, e.g., sampling and reconstructionin Conventional 1 and Conventional 2 are done slice-by-slice
for each ﬁxed elevation angle. Therefore, the measurements
were not the same as those for the proposed method, but compression rate was ﬁxed by α = 0.25 for all methods.
We evaluated each reconstruction method by the normalized error 100 ∥X∗ − X∥F /∥X∥F and the averaged computational time as shown inPTableP1, where
P the Frobenius norm
is deﬁned by ∥X∥F := ( nNRR=1 nNAA=1 nNEE=1 |x[nR , nA , nE ]|2 )1/2 ,and the numbers of iterations Imax were set to 650, 100, 30,
and 25 for Conventional 1, Conventional 2, Algorithm 1, andAlgorithm 2, respectively. From Table 1, the normalized errors for Data 1 by Conventional 1 and Conventional 2 were
16.94% and 13.61%, respectively. Those by Algorithm 1 andAlgorithm 2 were 9.75% and 9.78% for the entire reconstruction (s = 300), and were 9.88% and 9.90% for the blockwisereconstruction (s = 4). These results show that the proposed
3D reconstruction method outperforms the conventional 2D
methods. Moreover, we can see that the deterioration due tothe blockwise optimization was less than 0.5%, which is acceptable because this diﬀerence does not cause any change inthe estimation of rainfall (6) . Although Nesterov’s acceleration
caused the slight deterioration of the accuracy, this is also acceptable. Therefore, the proposed 3D CS technique achieves
both of low memory usage and high reconstruction accuracy.
Table 2 shows the normalized errors of Algorithm 2 with different s for Data 1. From Table 2, although the reconstruction
accuracy tends to become better for larger s, we judged thats = 4 is the best value among all factors of NA = 300 in termsof both the memory usage and the reconstruction accuracy.
Figure 3 shows the reconstruction results of the 10th elevation angle. The original data X and the compressed one Y

3.2 Nesterov’s Acceleration for Algorithm 1
In thissection, we improve the convergence speed of Algorithm 1 by
Nesterov’s acceleration technique (24)–(27) . Goldstein et al. proposed Fast ADMM as an accelerated variant of ADMM, andthe convergence rate can be improved from O(1/k) to O(1/k2 )for a strongly convex optimization problem (26) . Fortunately,
SDMM is just a special case of ADMM, and hence we apply
Nesterov’s technique to dual variables following the paper (26) .
Since the function in (2) is not strongly convex, we should
use Fast ADMM with a restart rule for stable computations(see Algorithm 8 in the paper (26) ). However, in our case, thisstrategy cannot accelerate Algorithm 1 because it takes muchtime, in every iteration, to judge whether the restart should bedone or not. Through our numerical experiments, we foundthat the restart hardly happens in the ﬁrst several iterations,e.g., ﬁve iterations for the problem in (2). Therefore, we apply Nesterov’s technique only in the ﬁrst P iterations without
checking the restart condition, and after that the same iterations as in Algorithm 1 are continued. Algorithm 2 summarizes this idea, where we use diﬀerent parameters γ1 and γ2for computing proxγgi to achieve the best convergence speed.
Note that Algorithms 1 and 2 are suitable for parallel computing. Hence, GPU implementations of these algorithms, inwhich the computation for each block Xl is assigned to eachcore, would further accelerate the reconstruction process.

Algorithm 2 Fast 3D PAWR data reconstruction
Input: Yl ∈ RNR ×s×NE , γ1 , γ2 ∈ (0, ∞), and integers P < Imax
1: yi,0 ← Li vec(Yl ) and zi,0 ← 0(i = 1, 2, . . . , 6)(i=
1, 2, . . . , 6)
2: yi, 1 ← yi,0 and zi, 1 ← zi,022
3: α0 ← 1
4: for k = 1, 2, . . . , P doP5:
xk ← Q−1 6i=1 LiT (yi,k− 12 − zi,k− 12 )6:
yi,k ← proxγ1 gi (Li xk + zi,k− 21 ) (i = 1, 2, . . . , 6)
7:zi,k ← zi,k− 12 + Li xk − yi,k (i = 1, 2, . . . , 6)√
1+ 1+4α2k−1
8:
αk ←
2
−1
9:
yi,k+ 12 ← yi,k + αk−1(i = 1, 2, . . . , 6)
αk (yi,k − yi,k−1 )−1(i = 1, 2, . . . , 6)
zi,k+ 12 ← zi,k + αk−1
αk (zi,k − zi,k−1 )11: end forγ2
12: zi,P ← γ zi,P(i = 1, 2, . . . , 6)1
13: for k = P + 1, P + 2, . . . , Imax do
P
14:
xk ← Q−1 6i=1 LiT (yi,k−1 − zi,k−1 )15:yi,k ← proxγ2 gi (Li xk + zi,k−1 ) (i = 1, 2, . . . , 6)16:zi,k ← zi,k−1 + Li xk − yi,k (i = 1, 2, . . . , 6)17: end for
Output: x∗l ← xk
10:

43

IEEJ Trans. EIS, Vol.140, No.1, 2020

3D Compression and Fast High-Quality Reconstruction for PAWR（Ryosuke Kawami et al.）(a) Original observed data X[ : , : , 10] [dBZ](b) Magniﬁcation of the red rectangle in (a)(c) 25% randomly selected data Y[ : , : , 10] [dBZ](d) Conventional 1 (normalized error: 16.35%)(e) Conventional 2 (normalized error: 11.93%)(f) Algorithm 2 with s = 4 (normalized error: 8.77%)(g) Magniﬁcation of the red rectangle in (d)(h) Magniﬁcation of the red rectangle in (e)(i) Magniﬁcation of the red rectangle in (f)Fig. 3. Simulation results for the reﬂection intensity X [dBZ] observed on March 30th, 2014.

are shown in Figs. 3(a) and (c). The reconstructed one X∗ by
Conventional 1, Conventional 2, and Algorithm 2 with s = 4
are shown in Figs. 3(d), (e), and (f). Magniﬁcated images of
X and X∗ are shown in Figs. 3(b), (g), (h), and (i). The normalized errors for this slice by Conventional 1 and Conventional 2 were 16.35% and 11.93%, respectively. On the other
hand, that by Algorithm 2 was 8.77%. From Fig. 3, we canconﬁrm high reconstruction quality of the proposed method.
Figure 4 shows the scatter plots of the simulation resultsfor Data 1. Figures 4(a), (b), and (c) show the results by Con-ventional 1, Conventional 2, and Algorithm 2 using s = 4,
respectively. We can see that the scatter plot of the proposed
method is narrower than those of the 2D methods. In particular, the plot of the proposed method more than 30dBZ seems
narrower. Table 3 shows the distributions of the normalized
errors in detail. Table 3 shows that the proposed 3D method
outperformed the other methods for regions especially morethan 30dBZ, 40dBZ and 50dBZ. These ranges indicate that itis heavily raining. Thus the proposed method may be eﬀective for such heavy rain regions. Figure 4 also implies that the
44

IEEJ Trans. EIS, Vol.140, No.1, 2020

3D Compression and Fast High-Quality Reconstruction for PAWR（Ryosuke Kawami et al.）Table 1. Normalized error and computational time of each reconstruction method.
Data 1 (2014/3/30)
Data 2 (2014/1/8)
Normalized Error [%]Data 3 (2013/6/19)
Data 4 (2013/5/10)
Averaged Computational Time [sec](Computational Time per Slice or Block)

Conventional 1
16.94
20.51
16.65
30.35
2,180
(19.82 per Slice)

Conventional 2
13.61
17.21
13.37
21.03
678
(6.16 per Slice)

Algorithm 1 (s = 300)
9.75
12.49
9.56
15.15
304
—

Algorithm 2 (s = 300)
9.78
12.48
9.62
15.33
223
—

Algorithm 1 (s = 4)
9.88
12.69
9.68
15.61
314
(4.19 per Block)

Algorithm 2 (s = 4)
9.90
12.65
9.72
15.82
242
(3.23 per Block)

30

50

60

70

50

60

70

50

60

70

Table 2. Normalized errors depending on s for Data 1.
s
Error

2
9.99

3
9.93

4
9.90

5
9.90

6
9.91

10
9.87

20
9.85

30
9.83

60
9.81

100
9.82

150
9.80

70

Reconstructed Values [dBZ]

60

Table 3. Normalized error distributions for Data 1.
Normalized Error [%]≥ 20dBZ≥ 30dBZ≥ 40dBZ≥ 50dBZ

Conventional 1
16.94
12.93
11.00
12.77
16.88

Conventional 2
13.61
9.38
7.42
7.68
8.03

Algorithm 2 (s = 4)
9.90
6.75
5.32
5.06
4.93

50
40
30
20
10
0

true PAWR data do not take the reﬂection intensities less than
15dBZ. This is because 15dBZ is too small to be valid measurements, and hence the system is designed to discard suchvalues. The proposed method does not use this knowledge,
and if do, the reconstruction quality can be further improved.
Table 1 also shows the computational time of each method.
From Table 1, we ﬁnd that both Algorithms 1 and 2 achieved
better reconstruction accuracy with less computational timethan Conventional 1 and Conventional 2. In addition, Algorithm 2 is approximately 1.3 times faster than Algorithm 1,
and hence Nesterov’s technique is eﬀective for acceleration.
We further conducted simulations for other data acquiredon January 8th, 2014, June 19th, 2013, and May 10th, 2013,
which are referred as Data 2, Data 3, and Data 4, respectively.
The normalized errors for each data are shown in Table 1 inthe same manner as Data 1. Further, Figs. 5, 6, and 7 show thereconstruction results of the 10th elevation angle for Data 2,
Data 3, and Data 4, respectively. From Table 1 and these ﬁgures, we can see that the proposed method robustly achieves
better reconstruction results than the conventional methods.
In the above simulations, the number of the blocks Xl is 75.
Yet, to complete the reconstruction for all thirteen parametersof the PAWR, the proposed method requires 3.23 × 75 × 13 ≈3,150 seconds. Since the reconstruction should be completed
within 30 seconds, 105 times acceleration is required approximately. To further accelerate Algorithm 2, implementation
by some compiler language like C is eﬀective. In addition,
since the reconstruction for each block can be performed in aparallel way, the computation based on GPU is also eﬀective.
If each block can be reconstructed by a separate computer inparallel, 105/75 ≈ 1.4 times acceleration is required to transfer all parameters observed by the PAWR data in real time.
5.

-10
-10

0

10

20

40

True Values [dBZ](a) Conventional 1

70

Reconstructed Values [dBZ]

60
50
40
30
20
10
0
-10
-10

0

10

20

30

40

True Values [dBZ](b) Conventional 2

70

Reconstructed Values [dBZ]

60
50
40
30
20
10
0

Conclusion

-10
-10

0

10

20

30

40

True Values [dBZ]

This paper proposed a fast high-quality 3D CS technique
for transferring the PAWR data in a public internet line. First,
we summarized our compression and reconstruction schemes
based on random sampling in the troposphere and blockwise
convex optimization, respectively. In the reconstruction part,
we deﬁned a new cost function that expresses a prior knowl-(c) Algorithm 2 with s = 4

Fig. 4. Scatter plot of each method for Data 1.

edge such as the local similarity and sparsity of the 3D PAWR
data in the spatial and frequency domains. Because the cost
45

IEEJ Trans. EIS, Vol.140, No.1, 2020

3D Compression and Fast High-Quality Reconstruction for PAWR（Ryosuke Kawami et al.）(a) Original observed data X[ : , : , 10] [dBZ](c) Conventional 1 (normalized error: 17.67%)(b) 25% randomly selected data Y[ : , : , 10] [dBZ](d) Conventional 2 (normalized error: 14.23%)(e) Algorithm 2 with s = 4 (normalized error: 10.68%)

Fig. 5. Simulation results for the reﬂection intensity X [dBZ] observed on January 8th, 2014.

function is convex, an eﬃcient reconstruction algorithm wasderived by SDMM. Then, for improvement of the computational speed, we applied Nesterov’s acceleration technique toSDMM. Numerical experiments using real data showed thatthe proposed method achieves low memory usage, high accuracy, and fast reconstruction compared with 2D methods.
Since the reconstruction can be performed for each block in aparallel way, our method can be further accelerated by GPU.
Acknowledgments
This work was supported by JSPS KAKENHI Grant Number 17H02069. The authors are grateful to Prof. Takashi Ijiri
in Shibaura Institute of Technology for fruitful discussions.

(5)

(6)

(7)

(8)

(9)

References
(1)

(2)

(3)

(4)

Committee on Extreme Weather Events and Climate Change Attribution:
Attribution of Extreme Weather Events in the Context of Climate Change,
Washington D.C.: The National Academy Press (2016)
F. Mizutani, M. Wada, T. Ushio, E. Yoshikawa, S. Satoh, and T. Iguchi: “Development of active phased array weather radar”, in Proc. 35th Conf. Radar
Meteor., Pittsburgh, PA: 4 pages (2011)
E. Yoshikawa, T. Ushio, Z. Kawasaki, S. Yoshida, T. Morimoto, F. Mizutani,and M. Wada: “MMSE beam forming on fast-scanning phased array weather
radar”, IEEE Trans. Geosci. Remote Sens., Vol.51, No.5, pp.3077-3088(2013)
F. Mizutani, T. Ushio, E. Yoshikawa, S. Shimamura, H. Kikuchi, M. Wada,
S. Satoh, and T. Iguchi: “Fast-scanning phased-array weather radar with angular imaging technique”, IEEE Trans. Geosci. Remote Sens., Vol.56, No.5,( 10 )( 11 )( 12 )

46

pp.2664-2673 (2018)
K. V. Mishra, A. Kruger, and W. F. Krajewski: “Compressed sensing applied
to weather radar”, in Proc. IEEE Int. Geosci. Remote Sens. Symp., Quebec
City, QC, Canada: pp.1832-1835 (2014)
S. Shimamura, H. Kikuchi, T. Matsuda, G. Kim, E. Yoshikawa, Y. Nakamura,and T. Ushio: “Large-volume data compression using compressed sensing
for meteorological radar”, Electron. Comm. Jpn., Vol.99, No.10, pp.103-111(2016)
R. Kawami, A. Hirabayashi, N. Tanaka, M. Shibata, T. Ijiri, S. Shimaura,
H. Kikuchi, G. Kim, and T. Ushio: “2-dimensional high-quality reconstruction of compressive measurements of phased array weather radar”, in Proc.
Asia-Pac. Signal Inf. Process. Assoc. Annual Summit Conf. (APSIPA ASC),
Jeju, South Korea: 7 pages (2016)
R. Kawami, A. Hirabayashi, N. Tanaka, M. Shibata, T. Ijiri, S. Shimaura,
H. Kikuchi, G. Kim, and T. Ushio: “Simple compression technique forphased array weather radar and 2-dimensional high-quality reconstruction”,
IEEJ Trans. Electron. Inf. Syst., Vol.137, No.7, pp.864-870 (2017)
R. Kawami, H. Kataoka, D. Kitahara, A. Hirabayashi, T. Ijiri, S. Shimamura,
H. Kikuchi, and T. Ushio, “Fast high-quality three-dimensional reconstruction from compressive observation of phased array weather radar,” in Proc.
Asia-Pac. Signal Inf. Process. Assoc. Annual Summit Conf. (APSIPA ASC),
Kuala Lumpur, Malaysia: pp.44-49 (2017).
E. Kudeki and F. Sürücü, “Radar interferometric imaging of ﬁeld-aligned
plasma irregularities in the equatorial electrojet,” Geophys. Res. Lett., Vol.18,
No.1, pp.41-44 (1991)
E. Yoshikawa, T. Ushio, Z. Kawasaki, S. Yoshida, T. Morimoto, F. Mizutani, and M. Wada: “MMSE Beam Forming on Fast-Scanning Phased Array
Weather Radar,” IEEE Trans. Geosci. Remote Sens., Vol. 51, No. 5, pp.30773088 (2013)
D. Kitahara, M. Nakahara, A. Hirabayashi, E. Yoshikawa, H. Kikuchi, andT. Ushio, “Nonlinear beamforming via convex optimization for phased arrayweather radar,” in Proc. Asia-Pac. Signal Inf. Process. Assoc. Annual Summit
Conf. (APSIPA ASC), Honolulu, HI: pp.1831-1835 (2018)

IEEJ Trans. EIS, Vol.140, No.1, 2020

3D Compression and Fast High-Quality Reconstruction for PAWR（Ryosuke Kawami et al.）(a) Original observed data X[ : , : , 10] [dBZ](c) Conventional 1 (normalized error: 17.54%)(b) 25% randomly selected data Y[ : , : , 10] [dBZ](d) Conventional 2 (normalized error: 13.53%)(e) Algorithm 2 with s = 4 (normalized error: 9.85%)

Fig. 6. Simulation results for the reﬂection intensity X [dBZ] observed on June 19th, 2013.
( 13 )( 14 )( 15 )( 16 )( 17 )( 18 )( 19 )( 20 )( 21 )( 22 )( 23 )( 24 )( 25 )( 26 )

D.L. Donoho: “Compressed sensing”, IEEE Trans. Inf. Theory, Vol.52, No.4,
pp.1289-1306 (2006)
E.J. Candès and M.B. Wakin: “An introduction to compressive sampling”,
IEEE Signal Process. Mag., Vol.25, No.2, pp.21-30 (2008)
Y.C. Eldar and G. Kutyniok, Eds.: Compressed Sensing: Theory and Applications, Cambridge, MA: Cambridge University Press (2012)
R. Baraniuk and P. Steeghs: “Compressive radar imaging”, in Proc. IEEE
Radar Conf., Boston, MA: pp.128-133 (2007)
M.A. Herman and T. Astronomer: “High-resolution radar via compressedsensing”, IEEE Trans. Signal Process., Vol.57, No.6, pp.2275-2284 (2009)
J.H.G. Ender: “On compressive sensing applied to radar”, Signal Process.,
Vol.90, No.5, pp.1402-1414 (2010)
P. Meischner, Ed.: Weather Radar: Principles and Advanced Applications,
New York, NY: Springer (2004)
T.Y. Yu, L. Ding, and S. Ozturk: “Application of compressive sensing torefractivity retrieval with a network of radars”, in Proc. IEEE Radar Conf.,
Kansas City, MO: pp.757-761 (2011)
D. Gabay and B. Mercier: “A dual algorithm for the solution of nonlinear variational problems via ﬁnite element approximation”, Comput. Math.
Appl., Vol.2, No.1, pp.17-40 (1976)
P.L. Combettes and J.C. Pesquet, Eds.: Fixed-Point Algorithms for Inverse
Problems in Science and Engineering, New York, NY: Springer (2011)
L. Condat: “A generic proximal algorithm for convex optimization—
Application to total variation minimization”, IEEE Signal Process. Lett.,
Vol.21, No.8, pp.985-989 (2014)
Y.E. Nesterov: “A method for solving the the convex programming problem
with convergence rate O(1/k2 )”, Sov. Math. Dokl., Vol.27, No.2, pp.372-376(1983)
A. Beck and M. Teboulle: “A fast iterative shrinkage-thresholding algorithm
for linear inverse problems”, SIAM J. Imag. Sci., Vol.2, No.1, pp.182-202(2009)
T. Goldstein, B. O’Donoghue, S. Setzer, and R. Baraniuk: “Fast alternating
direction optimization methods”, SIAM J. Imag. Sci., Vol.7, No.3, pp.1588-( 27 )

1623 (2014)
M. Yamagishi, D. Kitahara, and I. Yamada: “A fast dual iterative algorithm
for convexly constrained spline smoothing”, in Proc. IEEE Int. Conf. Acoust.
Speech Signal Process. (ICASSP), Shanghai, China: pp.4538-4542 (2016)

Ryosuke Kawami (Student Member) received the B.E. and M.E. degrees in computer science from Ritsumeikan University, Shiga, Japan, in 2017 and 2019, respectively. His
research interests include compressed sensing and its
applications. He received the Young Researcher Paper Award of IEEJ Transactions on Electronics, Information and Systems in 2018.

Daichi Kitahara (Non-member) received the B.E. degree in computerscience in 2012, and the M.E. and Ph.D. degrees incommunications and computer engineering in 2014
and 2017 from the Tokyo Institute of Technology,
Tokyo, Japan, respectively. From 2014 to 2017, he
was a recipient of the Research Fellowship of theJapan Society for the Promotion of Science (JSPS).
Since 2017, he has been an Assistant Professor with
the College of Information Science and Engineering,
Ritsumeikan University, Shiga, Japan. His current research interests are in signal processing and its applications, inverse problem,
optimization theory, and multivariate spline theory.

47

IEEJ Trans. EIS, Vol.140, No.1, 2020

3D Compression and Fast High-Quality Reconstruction for PAWR（Ryosuke Kawami et al.）(a) Original observed data X[ : , : , 10] [dBZ](c) Conventional 1 (normalized error: 27.16%)(b) 25% randomly selected data Y[ : , : , 10] [dBZ](d) Conventional 2 (normalized error: 18.55%)(e) Algorithm 2 with s = 4 (normalized error: 14.09%)

Fig. 7. Simulation results for the reﬂection intensity X [dBZ] observed on May 10th, 2013.
Akira Hirabayashi (Non-member) received his D.E. degree in computer science in 1999, from the Tokyo Institute ofTechnology, Tokyo, Japan. Currently, he is a Professor of the College of Information Science andEngineering, Ritsumeikan University, Shiga, Japan.
His research interests include compressed sensing,sampling theory, and signal and image processing.
Prof. Hirabayashi is also a member of the Institute ofElectrical and Electronics Engineers (IEEE) and the
Acoustical Society of Japan (ASJ).

Hiroshi Kikuchi (Member) received the B.S. degree from the Department of Engineering, Doshisha University, Kyoto, Japan, in 2008, and the M.S. and Ph.D. degrees
from the Division of Electrical, Electronic, and Information Engineering, Osaka University, Suita, Japan,
in 2010 and 2013, respectively. He joined the Division of Electrical, Electronic, and Information Engineering, Osaka University as a Specially Appointed
Researcher in 2013. In 2017, he was a Research Assistant Professor at Tokyo Metropolitan University. In
2018, he joined The University of Electro-Communications, where he is currently an Assistant Professor. His research specialties are the remote sensing
for an atmospheric electricity with space-borne platforms, the weather radar
remote sensing, and a development of the radar system.

Eiichi Yoshikawa (Non-member) received the B.E. degree in aerospace engineering from Osaka Prefecture University,
Sakai, Japan, in 2005, and the M.E. and Ph.D. degrees from Osaka University, Suita, Japan, in 2008
and 2011, respectively. In 2011, he was a PostDoctoral Researcher with Osaka University and Colorado State University, Fort Collins, CO, USA, andalso a Post-Doctoral Fellow of the Research Fellowship for Young Scientists sponsored by the Japan Society for the Promotion of Science. In 2012, he joined
the Japan Aerospace Exploration Agency (JAXA), Mitaka, Japan, where heis currently an Associate Senior Researcher. He is a Visiting Associate Professor with Tokyo Metropolitan University, Tokyo, Japan. His research interests include weather radar remote sensing, radar-based analyses, and applications for general and aviation weather.

Tomoo Ushio (Senior Member) received the B.S., M.S., and Ph.D. degrees in electrical engineering from Osaka University,
Suita, Japan, in 1993, 1995, and 1998, respectively.
He was with the Global Hydrology and Climate Center, Huntsville, AL, USA, as a Post-Doctoral Researcher from 1998 to 2000. In 2000, he joined theDepartment of Aerospace Engineering, Osaka Prefecture University, Sakai, Japan, as an Assistant Professor. In 2006, he joined the Division of Electrical,
Electronic and Information Engineering, Osaka University, as an Associate Professor. In 2017, he joined the Department of
Aeronautics and Astronautics, Tokyo Metropolitan University, as a Professor. Since 2019, he has been a Professor in Osaka University. His research
specialties are radar-based remote sensing, passive and active remote sensing
of atmosphere from space-borne platforms, and atmospheric electricity.

48

IEEJ Trans. EIS, Vol.140, No.1, 2020