Simplified Reflectivity Model and Its Application in Millimeter-Wave Radiometry

Millimeter-wave radiative sensing has been used in several close-range applications, such as security checks, military detection, terrain modeling, etc. Obtaining target information through radiometry is the goal of the above applications. However, there was no method to recover stable material information and complete geometric information through radiometry. The common problem of physically model-based information recovery is the solution to underdetermined equations. In this article, we analyze the multipolarization brightness temperature model and propose the equivalent permittivity (EP) to characterize the material information and simplify the reflectivity model. With the simplified reflectivity model, we solve the underdetermined problem in information inversion and realize the material classification and surface normal vector (SNV) reconstruction. Simulations and experiments are conducted to analyze the error and suitable range of the simplified reflectivity model and to verify the validity and accuracy of our methods of material classification and normal vector reconstruction. The results show that our classification method is suitable for most objects at an incident angle of 20°∼60°, and our SNV estimation method is effective at all incident angles. The possible applications include liquid composition analysis, target detection, road perception, and three-dimensional reconstruction.


I. INTRODUCTION
A LL substances at an absolute temperature radiate electromagnetic energy, this phenomenon is called thermal radiation. Millimeter-wave (MMW) radiometry is the measurement of electromagnetic radiation in the MMW band, and the sensors are referred to as radiometers [1]. With the ability to penetrate through a variety of common substances (e.g., atmosphere, fog, clouds, and to some extent rain, etc.) independently of the time of day, MMW radiometry has been used in astronomy [2], [3] and remote sensing [4], [5], [6], [7]. Since the 1990 s, with the development of MMW devices [8], MMW radiometry has been gradually applied to the information acquisition of close-range (compared with remote sensing) targets, such as the detection of guns, explosives, and other dangerous goods in security inspection [9], [10], the acquisition of road type, slope, and obstacles [11], environmental mapping and monitoring [12], [13], fire fighting [14] and military applications [15], [16]. Comparing the MMW radiation imaging results in remote sensing and close-range applications, we can see that, on the one hand, as the observation distance changes from at least hundreds of kilometers to hundreds of meters or even a few meters, there may be multiple clear targets in one image, such as guns, gasoline, human skin in the security inspection image, and water surface, cement ground, land, and various buildings in the ground image; on the other hand, the minimum resolution unit of the image can reach the centimeter or even millimeter level, so the objects in the image will show an obvious three-dimensional (3-D) structure. The above changes have also brought changes in the demand for information extraction. In remote sensing, we usually pay attention to the inversion of physical and chemical parameters of specific targets, such as atmospheric composition, seawater salinity, soil humidity [5], [6], [7], etc., but in close-range applications, we need to distinguish or identify different objects, so it is particularly important to obtain the material information and geometric information of objects. The polarization of electromagnetic radiation is governed by the geometrical configuration of the medium, and by its dielectric properties [17]. MMW radiometry with different polarizations (called polarimetry) has been successfully used for obtaining the parameters of the wind field, snow, ice, and hurricane [18], [19], [20], [21]. As radiometric resolution and spatial resolution increase, MMW polarimetry is applied for closerange objects. Liao et al. [22] developed a dual-polarization radiometer and imaged some objects, such as cars and forests. The experimental results showed that different objects had different characteristics of polarization difference. Welson et al. [23] also developed a dual-polarization radiometer, its imaging results showed similar polarization characteristics to [22]. Kin et al. [24] developed a multilinear polarization radiometer and imaged a metal ball, metal cup, and ceramic cup. The experimental results showed that there were obvious differences between these linear polarization imagings, and a higher contrast can be obtained by image processing. Cheng et al. [25] developed a multipolarization MMW radiometer and imaged the outdoor parking lot and the dormitory buildings. All the above studies show that there are obvious differences between the polarization brightness temperature (TB) of different objects, and the This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ polarization characteristics are dominated by material and surface orientation. That is, it is possible to infer some material and 3-D information by MMW polarimetry.
In the infrared and visible bands, the material classification, complex refractive index estimation, and 3-D reconstruction have been realized through polarimetry with the irradiation source or some prior information [26], [27], [28], [29], [30], [31], [32], [33], and the corresponding models and methods have been used for reference in the MMW band. Inspired by Partial Polarization (PP) [26], Degree of Polarization (DoP) [27], [29], Degree of Linear Polarization (DoLP) [28], [30], [31], and other feature parameters [32], [33] in the infrared and visible bands, some TB-derived feature parameters, such as Linear Polarization Ratio (LPR) [34], Passive Degree of Polarization (PDoP) [35], and Linear Polarization Difference Ratio (LPDR) [36], are proposed to material classification in MMW polarimetry. However, the above feature parameters are governed by the incident angle, so they no longer work when the observed angle is not fixed. The complex permittivity is only related to the material properties, so it is the best classification criteria. However, on the one side, we have realized the estimation of complex permittivity in [37], but this work requires too much prior information and multiangle measurements, which is not conducive to practical applications; on the other side, retrieving complex permittivity from feature parameter is an underdetermined problem. In the acquisition of 3-D information, Cheng et al. [38] realized the estimation of surface azimuth referring to Angle of Polarization (AOP) in infrared and visible bands. However, the surface normal vector (SNV) is determined by both azimuth and incident angle, and the estimation of incident angle is an underdetermined problem. The cause of the underdetermination problem is that there are too many parameters in the model on which the information inversion is based. There are only two solutions to the problem: one is to increase the number of measured parameters, but it will cause application limitations; the other is to reduce the parameter dimension of the model. In this article, the latter is chosen. Based on the analysis of the influence of complex permittivity on TB, the EP is proposed to replace the complex permittivity, and the reflectivity model is simplified. The EP can be used to characterize the material information and as the criteria for classification. Based on the simplified model, the estimation of incident angle can be realized without increasing the amount of measurement, so the reconstruction of SNV can be completed. The contributions of our work are summarized as follows.
1) Although having eliminated the influence of surrounding radiation, all material feature parameters in MMW radiometry are still related to the incident angle, so they are not stable parameters. In this article, we have created a new parameter EP, which can be retrieved from feature parameters, to eliminate both the influence of incident angle and surroundings. 2) EP can be used as the criteria in material classification.
To the best of our knowledge, this is the first work to realize material classification at different incident angles by MMW radiometry. The resolution and applicable range of EP have been analyzed by simulation and experiment. 3) The reflectivity model can be simplified by replacing the complex permittivity with EP. Based on the simplified model and multipolarization measurement, we have proposed a new reconstruction method of SNV. This is the first work of SNV reconstruction by monocular MMW radiometry, which does not require the active illumination source and the complex permittivity of the object. This article is organized as follows. Section II presents the fundamental theory of MMW radiation and then proposes the EP and the simplified reflectivity model. Section III presents the application of EP in material classification and analyzes the resolution and applicable range. Section IV proposes a new estimation method of SNV and verifies its feasibility and error through simulations. Finally, Section V concludes this article.

A. Fundamental Theory
MMW radiations of natural and artificial objects are in general partially polarized. Arbitrary polarization states can be decomposed into several linear polarization and circular polarization. Existing theories and experiments have shown that there is little circular polarization difference in most scenarios [39], [40], so we mostly focus on the linear polarization characteristics of radiation.
A pictorial illustration of MMW radiometry is shown in Fig. 1: a linear polarization radiometer is measuring an object with a horizontal surface. θ is the incident angle. The object radiation of different linear polarization can be measured by rotating the radiometer around the observation axis. α is the linear polarization angle (LP-angles), which determines the linear polarization state of radiation. According to the Rayleigh-Jeans approximation, the radiation intensity can be described with TB in the MMW band. Arbitrary linear polarization can be expressed as a linear combination of horizontal polarization and vertical polarization [41] where T h and T v are the horizontal polarization TB and vertical polarization TB, respectively. As shown in Fig. 1(a), the radiation received by the radiometer is partly due to the self-emission of the object and partly due to the reflected radiation originating from the surroundings, as expressed by (2) where = − i is the complex permittivity of object; T obj is the physical temperature of object; T inc is the unpolarized TB of surroundings; e h and e v are the emissivity of horizontal polarization and vertical polarization, respectively; r h and r v are the reflectivity of horizontal polarization and vertical polarization, respectively. Permeability is neglected in (2) because most objects are nonmagnetic. According to the Kirchhoff Laws and the conservation of energy, the emissivity can be calculated from where the reflectivity is given by the Fresnel's law [1] There are three parameters in (4), i.e., , , and θ. The next section will propose an equivalent reflectivity model based on EP, which reduces the dimension of parameters to two.

B. Simplified Reflectivity Model
In [37], we have revealed that the materials whose complex permittivities are located on the same "C-curve" have the same TB, i.e., have the same reflectivity. Let us give a clear definition to "C-curve": the curve formed on − plane by the complex permittivities with the same reflectivity at a certain incident angle and a certain polarization. That is, the "C-curve" through point ( 1 , 2 ) can be expressed as where p = h or v stands for horizontal polarization and vertical polarization, respectively. "C-curve" changes with polarization and incident angle, but we find that the change has special characteristics. To give an example, Fig. 2 shows the "C-curve" of water ( = 8.10 − 13.56i [1]) at different polarizations and incident angles. We can see that, the "C-curve" of horizontal polarization hardly changes with the incident angle, while that of vertical polarization is almost unchanged when θ < 60 • , but changes dramatically near the Brewster angle. In the complex permittivities on "C-curve," there is a special one, that is the intersection point of "C-curve" and -axis because it is a real quantity. It means that each complex permittivity on the same "C-curve" can be replaced by a real quantity e with no loss of accuracy r p (θ, , ) = r p (θ, e , 0).
We call e as the EP. As shown in Fig. 2, for a given complex permittivity, since its "C-curve" changes with polarization and incident angle, its EP is also related to polarization and incident angle. Fortunately, similar to "C-curve," EP is almost the same except near the Brewster angle at vertical polarization. Therefore, the reflectivity can be simplified with EP as and We call r h and r v as the simplified reflectivity, and [0 • , θ 0 ] is the suitable range of r v . Given the maximum tolerable error e MAX , the suitable range meets From the complex permittivity to the EP, it is mathematically a mapping from 2-D space to 1-D space. This mapping is shown in Fig. 3, where five "C-curve"s and the corresponding EPs are plotted. Although the EP given in Fig. 3 is at p = h and θ = 45 • , it is also applicable to other incident angles and polarizations except near the Brewster angle at vertical polarization as analyzed above. We have selected five complex permittivities on the same "C-curve" (E1 through E5, marked as a red dot in Fig. 3) and calculated their theoretical reflectivities as shown in Fig. 4. E5 is the EP of E1 through E4. We can see that, as already analyzed, the equivalent reflectivity has almost no error for horizontal polarization at all incident angles. However, for vertical polarization, the equivalent reflectivity is accurate  only for a certain range of incident angle, and the error of equivalent reflectivity increase with the incident angle until it reaches a maximum value at the Brewster angle, and then the error decreases.
For a lossless media, the incident angle, where r v goes to zero is called the Brewster angle [1]. Lossless means that the imaginary part of complex permittivity is 0. Therefore, from the strict definition, only E5 has a Brewster angle. However, most materials are lossy media, i.e., the imaginary part of complex permittivity is nonzero, so r v cannot reach zero at any incident angle. As an extension of the above definition, the angle where rv goes to the minimum (not always zero) is also called as Brewster angle. From the definition of the Brewster angle, the imaginary part of complex permittivity will affect the vertical polarization characteristics of reflectivity at the Brewster angle. We use EP (a real number) to replace the complex permittivity, which actually sacrifices the accuracy of the reflectivity model near the Brewster angle at vertical polarization, so as to ensure the accuracy of the model at other angles of vertical polarization and at horizontal polarization to the greatest extent. Therefore, the Brewster angle determines the center of the unsuitable range of incident angles. Or, to be exact, the center is between the Brewster angles of the complex permittivity and its corresponding EP, and the range depends on the distance of EP and the complex permittivity. For

A. Underdetermined Problem
In the previous section, we have treated the direct model, as shown in Fig. 5(a), where the inputs are TB-relevant parameters, and the output is the polarized TB. The following treats the inverse process, which uses TB, usually measured at different incident angles and polarizations, to produce the estimations of the interest parameter. Complex permittivities are the best choices to denote the material characteristic for object classification. However, the estimation of complex permittivity requires multiangle TB and T inc measurements, which is not conducive to practical applications [37]. So, some TB-derived feature parameters are proposed for material classification, which is easier to measure. For example, the measuring equations of LPR and PDoP are and which states that LPR and PDoP of an object can be obtained by measuring its T h , T v , and T obj at a certain incident angle. Putting (2) into (10) and (11), we can get the definitions of LPR and PDoP From (12) and (13), we can find that LPR and PDoP are functions of three quantities, i.e., , , and θ. and indicate that these feature parameters contain material information, so they are proposed for material classification. However, θ shows that these feature parameters are governed by the incident angle, so they no longer work when the observed angle is not fixed.
To realize the classification under variable incident angles, we hope to further retrieve the complex permittivity from these feature parameters based on (12) and (13). Since θ is usually known to the observer and LPR and PDoP can be measured by (10) and (11), there are still two unknown quantities ( and ) in (12) and (13). To estimate = − i , we need at least two independent equations. However, (12) and (13) are not independent of each other on account of The number of independent equations is less than the unknowns, so the estimation of complex permittivity from feature parameters is an underdetermined problem. In the next sections, we will use the simplified reflectivity model to reduce the dimension of unknowns, with which we solve the underdetermined problem.

B. Method for Material Classification
As shown in Fig. 3, using the simplified reflectivity model, complex permittivities in the two-dimension space are mapped to the 1-D space. That is, the number of unknowns in (12) and (13) is reduced to one PDoP Fig. 6 shows the theoretical calculation of LPR and PDoP for E1 through E5 (E5 is the EP of E1-E4). We can see that, within the suitable range of r v , (15) and (16) can work very accurately. Therefore, we no longer seek to invert the complex permittivity from (12) or (13), but instead invert the EP from (15) or (16) without an underdetermined problem. While reducing the dimension, EP retains the characteristics of complex permittivity to the greatest extent and is not affected by the incident angle and polarization within the suitable range, so it is suitable as the criteria basis for classification. Based on the above analysis, we propose a new material classification method with a three-step process. First, measuring the feature parameter; second, estimating the EP from the feature parameter based on (15) or (16); last, classifying different materials based on the EP estimations. As shown in Fig. 5(b), this is an inversion process from TB to object information. The core step of the above classification method is the estimation of EP. The accuracy of estimation will be analyzed in the later section. Here, we focus on the number of solutions to equations (15) and (16). Fig. 8 shows the PDoP as a function of EP. We can see that for a given incident angle, the PDoP increases with EP until it reaches a maximum value, and then the PDoP decreases. Therefore, there are two solutions of EP when PDoP and θ are given, and the two solutions are distributed on the two sides of maximum value. The same is true when estimating EP from LPR. We use water as an example to analyze the characteristics of the two solutions. According to the complex permittivity model of water [1], the complex permittivity of pure water at 20 • C and 94 GHz is 8.10 − 13.56i. According to (13), the PDoP of water is 0.3819 when θ = 50 • . Regardless of the measurement error, PDoP = 0.3819 and θ = 50 • are substituted into (15), and the two solutions are e = 0.620 and e = 21.74, respectively. The PDoP of both two solutions as a function of incident angle is shown in Fig. 7. It can be seen that 21.74 is the correct solution for EP, whose PDoP is very similar to that of water within a certain range; 0.620 is the wrong solution, which has the same PDoP value as the water only at θ = 50 • . The characteristic of these two solutions is that the Brewster angle of one solution is less than 50 • and that of the other one is greater than 50 • . This is because the PDoP as a function of incident angle is not monotonic. For the same PDoP value, it may appear before or after the Brewster angle.
As shown in Fig. 8, the smaller one of the two solutions is all in the range of e < 3 when θ ≤ 60 • . We can see from Fig. 3 that the EP of the most complex permittivity is greater than 3. In other words, when θ ≤ 60 • , the larger of the two solutions is the correct solution for the vast majority of cases. Therefore, we take the larger of the two solutions as the estimation of EP. This strategy is always effective when θ ≤ 60 • and e > 3. However, when θ ≥ 60 • , both solutions may be correct, so this strategy cannot be guaranteed to choose the correct solution. For the object of e < 3, the above strategy may also be effective, but the effective range will become smaller. The relationship between EP and the corresponding effective range can be found in Fig. 8: for example, for objects with e = 1, the effective range is θ ≤ 45 • . Based on the above analysis, the recommended incident angle range of our EP-based classification method is θ ≤ 60 • , within, which most objects can be accurately classified even measured from different directions.
For metals, r h = r v = 1, so the PDoP of metals is 0. As discussed in the appendix, considering the measurement error, we set the threshold value of PDoP to 0.01 to distinguish metals and nonmetals first, and then use the derived EP to classify nonmetals. To choose (15) or (16) to estimate EP, we need to consider the following two points. First, as shown in Fig. 6, LPR increases sharply near the Brewster angle, so it is very sensitive to measurement error. (16) should be selected in this case. Second, in other cases, the calculation cost of (16) is larger than (15), since its function form is more complex, so (15) should be chosen. Therefore, we set an empirical threshold of LPR = 50 to distinguish between these two cases. For LPR < 50, we choose (15), otherwise, choose (16). Based on the above analysis, our method for material classification is summarized as Algorithm 1.

Algorithm 1: Algorithm for Material Classification.
M and N stand for metal and nonmetal, respectively 1: else It is nonmetal 6: Calculate

C. Simulation
This section analyzes the influence of measurement error on the estimation of EP and further discusses the applicable range of incident angles of our method. LPR and PDoP are obtained by measuring T h , T v , and T obj . The sensitivity, that is the standard deviation (SD) of measured TB, is used to express the precision of measurement. It decides the SD of LPR and PDoP (ΔLPR and ΔPDoP), and then affects the SD of EP estimations (Δ e ). Fig. 9 shows the relationship between ΔPDoP and Δ e . We can see that, Δ e is proportional to ΔPDoP and the reciprocal of ΔPDoP.
The influence of error is analyzed by 500 independent Monte Carlo trials. Four materials, i.e., water, skin, concrete, and soil, are involved in the simulation. Table I shows their complex permittivities in 94 GHz obtained from published papers and corresponding EP with a suitable range calculated from the definition. We assume the infinite thickness of materials, which is usually reasonable. For example, Table I also gives the corresponding skin depth [1] of four materials. When the material thickness is 2 3 times the skin depth, it can be regarded as infinite thickness. In the simulation, T obj was set to 298 K; T inc was set to a set of measured data, as shown in Fig. 10. Then T h and T v at different incident angles were generated by (2). White Gaussian noises were added to T h , T v , and T obj to simulate the measurement error. Finally, the EP was estimated by Steps 1∼10 of Algorithm 1. Table II shows the mean and SD of EP estimations from 500 independent Monte Carlo trials. Two levels of noise with SD (σ N ) equaling 0.5 K and 1 K, respectively, and five different incident angles are included in the simulation. As analyzed in Section III-B, there are two solutions when estimating EP. We have discarded the smaller one when θ ≤ 60 • but reserve both two solutions when θ = 75 • . We can see that, for water and skin, θ = 75 • is out of the suitable range of EP. So, both the two solutions are wrong estimations of EP. On the other hand, for concrete and plaster, the suitable range is θ ∈ [0 • , 90 • ], so one of these two solutions is correct. However, the strategy proposed to eliminate the wrong solution in the previous section is effectless: the lower one is the correct solution. Therefore, our method is not recommended at θ > 60 • whether considers the suitable range of EP or eliminates the wrong solution correctly. Fig. 11 shows the theoretical ΔLPR and ΔPDoP calculated by (37) and (38). We can see that, when θ ≤ 60 • , SD of LPR and PDoP meets: soil > concrete>skin > water. But from Table II, SD of e meets: water > skin>concrete > soil. This is because, as shown in Fig. 8, ∂ e /∂PDoP increases with EP and plays a dominant role in (17). Similarly, we can find that SD of LPR and PDoP increases with incident angle; but in Table II, SD of e decreases with incident angle. This is also because, for a given EP, ∂ e /∂PDoP decreases with incident angle (see Fig. 8).
Overall, the higher the incident angle is and the lower EP is, the higher the precision of EP is.
When θ ∈ [30 • , 60 • ], the mean of EP is almost stable with the incident angle and the SD of EP estimations is within a Fig. 11. Theoretical calculation of ΔLPR and ΔPDoP with ΔT = ΔT obj = 0.5 K for water, skin, concrete, and soil. small range. These two points are conducive to the use of EP for material classification. When measurement error is doubled from 0.5 to 1 K, the SD of EP estimates is also doubled, not exponentially. This proves the robustness of the EP estimation. However, when θ = 15 • , the result is different. If no noise is added, we can estimate EP with very high precision. But as shown in Table II, when the noise of 0.5 K is added, the SD of EP estimations becomes very large; when noise is doubled, the SD of EP increases exponentially, and the mean of EP estimations deviates from the theoretical value. This is because, when θ = 15 • , ∂ e /∂PDoP becomes very large (see Fig. 8), that is, the curve in Fig. 9 is very flat. This leads to two consequences. First, the SD of EP estimates increases exponentially with measurement error. Second, the mean of EP estimations is larger than the theoretical value. To summary, the best-applied range of our method is θ ∈ [30 • , 60 • ]. In fact, the upper limit of the best-applied range depends on the Brewster angle, so it is material-based, and we only give a general estimate for the Brewster angle of most materials is greater than 60 • . Our method has no lower limit of the suitable range in theory, that is, our method is also applicable to small incident angles in theory. However, because the polarization characteristics are weak at a small incident angle, all the polarization-based methods face the problem of accuracy degradation at a small incident angle. Therefore, the lower limit of our method depends on the accuracy of measurement. If Algorithm 1 is used, it is recommended to use our method only when the incident angle is greater than 30 • . When θ < 30 • , our method is feasible in theory, but it needs higher measurement accuracy or a new measurement method of PDoP/LPR.
Histogram of EP estimations of θ = 30 • , 45 • , and 60 • are shown in Fig. 12. When σ N = 0.5 K, as shown in Fig. 12(a),  estimations are distributed in four clusters (except for a pulse of e = 3), corresponding to four materials. For skin, concrete, and plaster, there is only one peak in each cluster; but for water, there are two peaks. This is because, as analyzed in Section II-B, EP changes slightly with incident angle, so derived EP from different incident angles also has a little difference. When the incident angle is within the suitable range, this difference is very small, as the Mean e of skin, concrete and plaster has shown in Table II. But 60 • runs a little over the suitable range of water, so Mean e of 60 • slightly deviated from that of 30 • and 45 • . Therefore, there are two peaks in the EP estimations of water. However, when the noise increases to 1 K, the SD of estimations is greater than the difference of Mean e , so the two peaks are blurred into one peak, as shown in Fig. 12(b). On balance, high precision of measurement is propitious to material classification, but if the precision is too high, the same material measured at different incident angles may be divided into different categories. Except for a pulse of e = 3, we can see from Fig. 12 that the estimations can be regarded as a Gaussian mixture distribution. The pulse is formed because the PDoP of plaster is 0.978 at θ = 60 • , but due to the measurement error, the measured PDoP sometimes exceeds 1; as shown in Fig. 8, the best estimations of EP are 3 when PDoP> 1, then all these estimations form the pulse in Fig. 12. Therefore, we first eliminate the data of EP= 3, then use the remaining data to fit the Gaussian mixture model where k is the number of categories; α i , μ i , and σ i are the weight, mean, and SD of the ith Gaussian model, respectively. Finally, all the data are classified with the fitted model. Fig. 13 shows the fitted Gaussian mixture model. When σ N = 0.5 K, we give two results with k = 4 and k = 5. We can see that their difference in fitted results is that when k = 5, EP estimations of water are regarded as a combination of two Gaussian distributions. Compared with the fitted result of σ N = 0.5 K, we can see that the μ i in the fitted result of σ N = 1 K is almost unchanged, but the σ i becomes larger.
The material classification results based on the fitted Gaussian mixture model are shown in Tables III-V. We can see that, based on the EP estimations, water is distinguished from the other three materials. This is because, as shown in Table II, the EP difference between water and the other three materials is much greater than the SD of EP estimations. However, as analyzed previously, if the measurement precision is too high (i.e., σ N is small), the same material may be divided into different categories, as  shown in Table IV. The EP difference between skin, concrete, and plaster is about 2, which is 3 ∼ 4 times greater than SD e for σ N = 0.5 K, and 1 ∼ 2 times greater for σ N = 1 K. Therefore, most estimates are classified correctly.

D. Experiment
In the previous section, we analyzed the influence of measurement noise on EP estimation and material classification. It can be seen that EP can be used to characterize material information and serve as the criteria of classification, but the random noise in measurement will cause some proportion of classification errors, especially at a small incident angle. When the type of materials increases, the corresponding EPs are closer to each other, then the proportion of error will further increase. To sum up, the measurement error limits the applicable range of incident angles and determines the sensitivity of classification.
Hu et al. [42] defined a new feature parameter RDoP (Reflecting DoP) and proposed a physically-based method using the characteristic of weak correlation of the emission part and reflection part of MMW radiation as an optimization criterion to obtain the optimal estimation of RDoP, which can theoretically eliminate the influence of most measurement errors, such as the thermal noise of the radiometer, radiative transfer, antenna pattern, and calibration error. Both PDoP and LPR can be calculated from RDoP Since RDoP eliminate the influence of most measurement errors, PDoP and LPR are almost not affected by measurement error.  Based on the above analysis, the measurement method of PDoP and LPR in Algorithm 1 is updated to Algorithm 2.
Since the measurement method of RDoP used in Algorithm 2 has eliminated most measurement errors, the measurement error can be reduced from 13% to 3% compared with traditional methods of (10) and (11), [42]. According to (17), Δ e also has the same degree of improvement. Because Δ e determines the resolution of classification, Algorithm 2 has much more application potential. In addition, Algorithm 2 does not need to measure the physical temperature of the object, does not need calibration, and has lower requirements for the sensitivity and resolution of the radiometer, which is more suitable for practical applications. This section will use Algorithm 2 to estimate the RDoP of water, alcohol, and gasoline, then estimate the EP and classify these materials, so as to show the best sensitivity of our EP-based classification method and its application potential in liquid detection with the existing measurement technology.
The experimental scene is shown in Fig. 14. The radiometer is installed on the turntable, and TB at different incident angles can be measured through the pitch and lift control of the turntable. The measuring objects include water, alcohol, and gasoline in the container, which is surrounded by radar absorbing materials (RAM) to reduce the influence of radiation outside the object. On the premise that the radiometer can receive the detectable radiation of sample, the area of the container should be as small as possible, so as to avoid the smoothing effect of the antenna pattern on the measured feature parameter. We measured the RDoP at five incident angles from 20 • to 60 • , and repeated it ten times at each incident angle to investigate the stability of the measurement results. Table VI shows the mean and SD of RDoP measurements, and Table VII shows the mean and SD of EP estimated by Steps 3∼11 of Algorithm 2. In material classification, we pay more attention to the SD of measurement, which determines the sensitivity of classification. It can be seen that the SD of EP estimations in the experiment is significantly smaller than the simulation results, and we can get very stable EP estimations even at a small incident angle. This is due to the RDoP measurement method used in the experiment, which can eliminate the influence of random error.
The three objects measured in the experiment, respectively, represent three typical characteristics of complex permittivity. The characteristic of water ( = 8.1000 − 13.5616i) is that the imaginary part of its complex permittivity is very large. According to the analysis in Section II-B, EP will also be affected by the incident angle, but within the suitable incident angle range, EP can be considered to be unchanged. The larger the imaginary part, the smaller the suitable range. Since 60 • is at the boundary of the suitable range, the EP at θ = 60 • is slightly larger than other incidence angles. This is consistent with the simulation results in Table II. The characteristic of alcohol is that its imaginary part is small, so it has a wide suitable range of incident angles. However, its Brewster angle is exactly at 60 • . Since the Brewster angle is the place where the error of equivalent reflectivity is the largest, the EP at θ = 60 • will be slightly different from that at other incident angles. The characteristic of gasoline is that its EP is less than 3. According to the analysis in Section III-B, for the objects with EP less than 3, the estimation of EP may be affected by the multisolution problem and get the wrong solution. This is proved when θ = 60 • . Fig. 15 shows the histogram of EP estimations. It can be seen that the EP estimations are clustered into four categories. Because the RDoP measurement method used in the experiment can effectively suppress the random error, there is no overlapping between each cluster, so there will be no classification error. Affected by the multisolution problem, the estimated EP of gasoline at 60 • is wrong, so it is classified into a single category, besides which, water, alcohol, and gasoline are correctly classified into three categories.
To conclude, we have proposed and demonstrated a new material classification method by multipolarization MMW radiometry. The experiment result shows that our method is applicable to objects with different characteristics of complex permittivity. When θ ∈ [20 • , 60 • ], most targets can be accurately classified with high sensitivity. When θ < 20 • , there is almost no polarization difference of object, so our polarization-based method is no longer applicable. When θ > 60 • , there may be two solutions when estimating EP, which may cause classification errors and need to be addressed in the future. In addition, it should be noted that our method cannot distinguish the materials on the same C-curve. As shown in Fig. 5 of [37], different complex permittivities on the same C-curve have very similar  MMW radiation. They can be regarded as one type of material with similar characteristics. Unless their vertical polarization brightness temperature is measured near Brewster angle, they are the same material for the observer. Based on the TB of vertical polarization near the Brewster angle, they can be further distinguished.

A. Radiometry Model of a Stereo Object
The incident plane is defined as the plane containing the SNV (n) and the direction of observation (I). The incident angle in (2) can be expressed by θ = arccos(I · n).
The polarization with electric field perpendicular to the incident plane is called horizontal polarization, and that with electric field parallel to the incident plane is called vertical polarization Above horizontal and vertical polarization are the polarization bases of object radiation in (2). Fig. 16(a) shows the geometry of the radiometry model of a stereo object, which is more complex than that in Fig. 1(a). Natural and artificial objects can be equivalent to the combination of many small surfaces with different SNV, so we only draw one small surface in Fig. 16(a). Since n is unknown to the observer, (22) is not applicable to the radiometry as the polarization bases of measuring TB. Therefore, as shown in Fig. 16(a), the polarization bases of the radiometer are According to (23), the polarization bases of the radiometer are exclusively governed by the observation direction. Different linear polarization measurements can be realized by rotating the antenna around the observation direction. As shown in Fig. 16(b), the polarization rotation angle (PR-angle, β) is defined as the rotation angle of antenna/waveguide relative to V, which determines the electric field direction of received radiation. In general, the two polarization bases, i.e., (h, v) and (H, V), are different, unless z is parallel to the incident plane. We define the relative angle between these two polarization bases as phase angle (P-angle, α 0 ). P-angle depends on both the observation direction and the surface orientation of the object [41]. As shown in Fig. 17, the use of the local coordinate system (V, H, I) can clearly describe the relationship between the incident angle, the P-angle, and the SNV. The incident angle is the angle between the SNV and the observation direction, so it is the zenith angle of n in (V, H, I). Furthermore, both h and v lie on the H − V plane, and h is perpendicular to the plane formed by I and n, so h is perpendicular to the projection of n on the H − V plane. Fig. 17 shows this relationship explicitly, and it is easy to see that the azimuth of n in (V, H, I) is 180 • − α 0 . To sum up, in the local coordinate system, the incident angle determines the zenith angle of SNV, the P-angle determines the azimuth angle of SNV, and the SNV can be expressed as n = (−sinθcosα 0 , sinθsinα 0 , cosθ). Of course, we are more concerned with the SNV in the global coordinate system (x, y, z). As shown in Fig. 16(a), given the observation direction I(θ i , ϕ i ), (V, H, I) can be obtained by first rotating (x, y, z) counterclockwise around z with ϕ i , and then rotating counterclockwise around y with θ i . Therefore, the SNV in (x, y, z) can be expressed as where R z and R y are the rotation matrix of coordinate system from (V, H, I) to (x, y, z) It can be seen from (24) that, since θ i and ϕ i are known to the observer. If we can realize the estimation of P-angle and incident angle, the reconstruction of SNV can be realized.

B. Underdetermined Problem 1) P-Angle Estimation:
The measured TB at different PRangles can be expressed as [41] T B (β) = Qcos[2(β − α 0 )] + I (27) where Q and I are defined as From (27), we can know that T B has a cosine variation with β, and α 0 determines the initial phase of the cosine curve. By measuring T B at three or more PR-angles (i.e., β i , where i ∈ {1, 2, . . . , N}, N ≥ 3), we can get a set of nonlinear equations α 0 , Q, and I can be estimated from (29). To prevent multiple solutions, we set a constraint condition of Q < 0 for T h < T v in most scenarios. Fig. 18 shows a typical cosine curve of T B as a function of PR-angle. For T B | α 0 =α 1 = T B | α 0 =α 1 +180 • , there is a 180 • ambiguity in the estimation of α 0 .
2) Incident Angle Estimation: From (21), we know that the incident angle θ is the zenith angle of n in the coordinate system of (V, H, I). Q and I can be estimated from (29). Then, T h and T v can be calculated from (28). For T obj and T inc can be obtained by measurement, there are there unknowns (i.e., θ, and ) in (2), but (2) contains only two independent equations for p = h and v, respectively. Therefore, the estimation of incident angle is an underdetemined problem.
In Section II, we have proposed the EP to characterize the material information and to replace the complex permittivity in the reflectivity model. Using the simplified reflectivity model, the unknowns in (2) are reduced to two So θ and e can be solved from the above two nonlinear equations. The error caused by this approximation, especially near the Brewster angle, will be analyzed in Section IV-C.

3) SVN Estimation:
After the estimation of P-angle and incident angle, the SNV can be calculated by (24). The estimation method of SNV is summarized as Algorithm 3.

Algorithm 3: Algorithm for SNV Estimation.
Input: Algorithm 3 is only a general framework of our SNV estimation method, and more details will be given in the following section.

C. Error Analysis
The error of SNV estimation mainly comes from two factors: the measurement error of T B , T obj , and T inc ; the approximation error of (30) when estimating θ. Therefore, a three-step analysis is adopted to discuss the error. First, the influence of measurement error on the estimation of α 0 , T h , and T v is analyzed. Second, the approximation error of (30) is analyzed. Finally, the total error of SNV estimation is analyzed.
1) Measurement Error: The sensitivity (ΔT ) indicates the minimum change of TB that can be detected by the radiometer. The difference between T B at different PR-angle should be much greater than the sensitivity, and a large difference is beneficial to the robustness of estimation.
From (29), the difference between two T B measurements (β = β 1 , β 2 respectively) can be expressed as where Δβ = β 1 − β 2 . On the one hand, Δβ should be as larger as possible to increase ΔT B . So, we should evenly select the PR-angle from 0 • to 180 • (e.g., β = 0 • , 60 • , 120 • when N = 3; On the other hand, |Q| should be as large as possible. From (2), Q can be expressed by We can see that a large difference between T obj and T inc is needed to improve the estimation precision. So, the choice of ambiance is important. In outdoor measurements, the sky is a good choice. However, in indoor measurements, an artificial ambiance is needed, just as [43], [44] have done. In addition, as shown in Fig. 4, θ influences the difference between r h and r v , and then influences Q. So, a small θ makes against the robustness of estimation. Based on the above analysis, we designed the Monte Carlo simulations. Table VIII shows the parameters set. θ was set from 10 • to 80 • to analyze its influence. Two groups of β were set to analyze the influence of the number of PR-angles. These parameters were input to (2) and (29) to generated T Bi . To simulate the measurement noise, the white Gaussian noises with SD = 0.5 K were added to T Bi . Finally, the synthetic datasets of (T Bi , β i ) were input to the Levenberg-Marquardt solver to estimate α 0 , T h and T v from (29). Table IX lists the mean and the root-mean-square error (RMSE) of the estimations of 500 independent Monte Carlo Simulations. In this article, to distinguish from the true value, we add · to the original symbol to represent its estimated value (e.g., α 0 represents the estimation of α 0 ). We can see that the accuracy of estimations of α 0 is improved with θ increases. This is consistent with the previous analysis. In contrast, the accuracy of T h and T v estimations are pretty high and almost invariant with incident angle. This is good news for SNV reconstruction, since the accuracy of T h and T v determines the estimation error of the incident angle. Besides, comparing Group 1 and Group 2, it is clear that the accuracy can be improved by increasing the measuring number of P-angles.
2) Approximation Error: Using (30), θ can be solved from the estimated of T h and T v . To make clear the distribution of the solution of (30), we define the error function When we analyzed the approximation error, θ were set differently from 10 • to 80 • , and T obj , T inc , and were set the same as Table I to calculate T h and T v from (1). Then, we went through θ and e to calculate e(θ, e ). The zero points of e(θ, e ) are the solutions of (30). Fig. 19 shows the calculation results and marks the zero point. As shown in Fig. 19, although we cannot obtain the accurate complex permittivity from (30), θ can be effectively estimated. The estimation error is about 0.02 • ∼ 0.18 • for θ i < 60 • and then rises to 0.79 • and 0.45 • for θ i = 70 • and 80 • , respectively. These errors are completely caused by the approximation error of (7) and (8). As shown in Fig. 4, the approximation error reaches its maximum near the Brewster angle (about 70 • − −80 • ), so the estimation error of the incident angle also reaches its maximum near the Brewster angle.
In addition, Fig. 19 exposes the multisolution problem when estimating θ from (30), i.e., there may be two zero points in (33). One of these two solutions appears before the Brewster angle, and the other one appears after the Brewster angle. This problem is especially obvious when θ is large. As an example, the TBs of e = 12.09 and e = 3.71 (the two solutions when θ = 60 • ) are calculated from (2) and then compared with the true TB of e = 10 − 5i, as shown in Fig. 20. We can see that the TB of e = 12.09 almost coincides with the true TB at all incident angles, but only when θ = 75.6 • is the TB of e = 3.71 same as the true TB at θ = 60 • . So e = 12.09 is the right solution of (30). We can find from Fig. 19 that the right solution of e is stable in a small range (about 11.03-12.82) while the wrong one changes sharply with θ (about 2.30-33.12). We can use this characteristic to eliminate the wrong solution.
3) Total Error: After estimating α 0 , T h , and T v from the synthetic datasets of (T Bi , β i ), we input T h , T v , T obj , and T inc to the solver to estimate θ from (30). To simulate the measurement noise, the white Gaussian noises of standard deviation equal to 0.5 K were added to T obj and T inc before they were input to the solver. Table X lists the simulation results of 500 independent Monte Carlo Simulations. The SD θ indicates the error caused by measurement noise, and the |Mean θ − θ| indicates the error caused by the systematic error, i.e., by the approximation error. The RMSE θ indicates the total error. We can see that the total error of θ decreases with the increase of the incident angle. This is because the measurement error dominates the total error when θ < 60 • . However, when θ > 60 • , the approximation error increases with the incident angle until it reaches a maximum value, and then decreases, so the total error increases first and then decreases. We can also see that e is kept within a certain range, and the fluctuation caused by the measurement error is less than 0.4, which is far less than the fluctuation of the wrong solution in Fig. 19. Therefore, the elimination method of the error solution proposed above is feasible.
e SNV is the error of SNV estimation and defined as the angle between the true SNV and its estimation e SNV = arccos(n · n).
Using (24), (34) can be rewritten as Since the rotation of coordinate system does not affect the angle between n and n, R z and R y are ignored in (35). We can see from Table X that the error is significantly larger at θ = 10 • than at other angles. This is due to the large errors of estimations of θ and α 0 at small incident angles. Fortunately, as shown in Table IX, although there is a huge error in the estimation of α 0 when θ i = 10 • , the expectations of e SVN are only 3.34 • and 2.32 • for Groups 1 and 2, respectively. This is because when θ is small, e SVN mainly depends on the error of estimation of θ. In addition, e SVN is also large near the Brewster angle. Compared with Group 1, e SVN of the Group 2 does not decrease significantly. These are because, when θ nears the Brewster angle, the e SVN is mainly determined by the approximation error which is irrelevant to the measuring number of PR-angles. The 180 • ambiguity problem of α 0 estimation and multisolution problem of θ estimation are not considered in Tables IX and X. The solution to these two problems requires an overall consideration of the object characteristics and will be applied in the next section.

D. Simulation
The fast ray-tracing method [45] is used to simulate the MMW radiation imaging of a sphere. We chose the sphere as the object because it is very representative: its surface contains all the incident angles of 0 • ∼ 90 • and all the P-angles of 0 • ∼ 360 • . Fig. 21 shows the imaging scene. A sphere whose radius is 0.5 m is in the center of the field of view (FOV). The FOV and distance of imaging are 20 • × 20 • and 3 m, respectively. T obj , T inc , and are set the same as Table VIII. The antenna pattern smoothing and Gaussian noise are applied to the imaging to simulate the spatial resolution and sensitivity (ΔT ) of a radiometer. The spatial resolution is fixed to 0.4 • . ΔT is set to 0.5 and 2 K, corresponding to low and high noise, respectively. As shown in Fig. 22, 12 images with 100 × 100 pixel are obtained by simulation and divided into two groups: the low noise group and the high noise group. Each group consists of six imagings, and each imaging is obtained at one PR-angle (from 0 • to 150 • with the step of 30 • ). α 0 , T h , and T v of each pixel are obtained by fitting the simulation data to (29). As shown in Fig. 23, there are two α 0   (30), e and θ are solved from T h and T v . There may be two solutions to (30). As analyzed in Section IV-C, e of the true solution is stable in a small range. Using this characteristic, we have eliminated the wrong solution. Fig. 24 shows the estimation results of e and θ. As shown in Figs. 24(b) and(d), most e is near 12; a circle of maximal values appears near the circumference, and then the minimum values appear at the circumference. These are consistent with Fig. 19: the Brewster angle appears near the circumference, where e reaches its maximum.
Because the FOV of imaging is 20 • × 20 • , according to the coordinate system in Fig. 21, the azimuth angle and zenith angle of observation of each pixel can be expressed as where m and n are the row and column of the pixel. Finally, the SNV of each pixel is calculated by (24). For the best visual effect, we evenly take 30 × 30 pixels from the calculations of SNV and show them in Fig. 25. Due to the 180 • ambiguity of α 0 , there are two estimations of SNV for each group. Intuitively comparing (a), (c) and (b), (d) of Fig. 25, we can find that the former is convex while the latter is concave, and both are the possible estimations of SNV. So, we need to use some prior information about objects to eliminate the wrong one. For example, the surface of an object is usually convex, or SNV at the edge of the object is usually outward. so Fig. 25(a) and (c) are the right estimations. Fig. 26 shows the estimation error of SNV calculated by (34). We can see from Fig. 26(a) and (c) that, the error is large near the center of the image, and it gradually decreases as the  pixel moves away from the center, and suddenly increases near the circumference. These are consistent with the analysis in Section IV-C: the incident angle near the center is small, so the error is large; as the incident angle increases when the pixel moves away from the center, the error decreases; At the edge, on the one hand, T B is smaller than the true value due to the antenna pattern smoothing, on the other hand, the approximation error increases due to the incident angle is close to the Brewster angle, so the error suddenly increases. Fig. 26(b) and (d) show the error of wrong estimations of SNV. The wrong estimation of α 0 is 180 • away from the right one, so the spatial relationship of the wrong SNV and real SNV is shown in Fig. 27. Obviously, the error of the wrong SNV estimation is about 2θ, which has been confirmed by Fig. 26(b) and (d).
The simulation has verified the validity of the proposed SNV estimation method. For most SNV estimations, the accuracy is about 1 • with the typical sensitivity of 0.5 K. Our method still keeps a relatively high accuracy of about 2 • with the sensitivity of 2 K, which shows the robustness of our method. However, when the incident angle is small, the error becomes large. This error is mainly caused by the measurement error, which can be reduced by increasing the measurement times or improving the sensitivity of the radiometer.

V. CONCLUSION
A new parameter EP has been created to characterize the material information of the object. By replacing the complex permittivity with EP, we can simplify the reflectivity so as to reduce the dimension of model parameters. Based on the simplified reflectivity model, we have solved the underdetermination problem in the material classification and the SNV estimation. The main works of this article and research prospects are concluded as follows.
Aiming at the instability of feature parameters to represent the material information, this article proposes to use the EP derived from feature parameters to act as the criteria of material classification. The EP hardly changes with the incident angle, which is suitable for object classification and recognition. Simulation and experiment results show that our method is applicable to most objects with different characteristics of complex permittivity, and the applicable range of incident angle can reach [20 • , 60 • ] with the current measurement method and precision.
This article also proposes an SNV reconstruction method that can be realized by monocular radiometry. In our method, the SNV is decomposed into two parameters: P-angle and incident angle. The P-angle can be derived by multipolarization measurement, and the incident angle is derived based on the radiation model with simplified reflectivity. Simulation results have verified the robustness and effectiveness at all incident angles.
The research carried out in this article provides theoretical and methodological support for using MMW radiometry to obtain object information, and is conducive to the application of MMW radiometry in the fields of liquid composition analysis, human security inspection, road condition perception, environmental monitoring, and military target detection. However, there are still some limitations that need further research.
First, although we have demonstrated the classification of three objects with typical complex permittivities, classifying more materials is an important follow-up work, which requires higher classification sensitivity, that is, smaller Δ e , which puts forward higher requirements for measurement accuracy. Based on our method, we are preparing a liquid composition analysis equipment, using noise sources to irradiate and designing the radiation environment, so as to further improve the material resolution. We look forward to identifying different kinds of dangerous liquids and even gasoline of different labels.
Second, the above methods are only suitable for objects with a smooth surface. In the infrared and visible bands, part of the physically-based material classification methods and most of the physically-based 3-D reconstruction methods are also only suitable for smooth objects. Due to the longer wavelength of MMW, the smoothing objects in the MMW band are broader than that in infrared and visible bands. Therefore, our methods still have a wide range of applications. However, if the MMW radiation of objects with different roughness can be accurately modeled, on the one hand, the smooth surface model can be modified for micro rough surfaces to improve the estimation accuracy of interest parameters. On the other hand, the methods in this article can also be extended to objects with a rough surface, so as to expand the scope of application of our method.
Last, our SNV estimation method has not been verified by experiments, this is mainly due to the lack of a high-resolution imaging radiometer. Our method requires the measurement of high-accuracy TB of the objects. However, when the resolution of the radiometer is insufficient, on the one hand, the environmental radiation will enter the radiometer through the sidelobe of the antenna, on the other hand, the radiation of neighboring pixel will also enter the radiometer, which will affect the accuracy of TB measurement and cause the error of SNV reconstruction. In the infrared and visible bands, the model-based 3-D reconstruction of objects is usually based on the assumption of a smooth surface. Due to the shorter wavelength, the requirements for a smooth surface are more stringent than that of the MMW band, but these 3-D reconstruction methods have still been successfully applied in practice because of the high resolution of the camera. Our team is trying to improve the resolution of the radiometer from both the hardware and the image processing, we believe that with the further development of the MMW radiometer in resolution, our SNV reconstruction method will also emerge in various applications.