Composite materials have a high level of uncertainty (intrinsic and non-intrinsic) due to the manufacturing process as well as the placement of different phases of their constituent materials. These uncertainties can be identified in both macro and micro scales. Identifying the behavior of structures made of composite materials without taking into account the uncertainties, whether due to identification or modeling, can lead to unrealistic results, especially in the dynamic behavior of structures. One of these cases is the identification of damage types in composite structures which is usually done by using dynamic responses. Damages in composite materials or structures usually occur during construction or operation. The correct modeling of uncertainty sources is one of the most important factors in identifying the geometry, location, and severity of damages accurately. The uncertainties related to the position and placement of carbon nanotubes (CNTs) can cause noticeable changes in the characteristics of composite materials reinforced with CNTs. For this reason, in the present study, we identified damages in CNT panels by considering all possible sources of uncertainty. A probabilistic multi-stage reliability-based method was proposed in this study to detect damage in these structures. In order to model the intrinsic and non-intrinsic sources of uncertainty, a modified point estimation method (MPEM) was used. In addition, an enhanced differential quadrature (DQ) method was used to model the CNT panels. In each step of the proposed algorithm, the probability of damage in each element of the panels was calculated by analyzing the possible damages. According to the results of the previous step, the elements with a low failure probability were gradually sifted in the next steps. The sieved elements in each step were considered as intact elements in the next step. This systematic filtering of design variables can simultaneously reduce the dimensions and speed up the optimization problem. Finally, the probability of damage was calculated based on the probability density function of various damage severities and positions. The developed approach was applied for damage detection on a laboratory-tested plate to illustrate the efficiency of the proposed method. The effects of using different damage positions and severity levels on the diagnosis results were discussed. The results demonstrated that the number of frequencies and modes of vibration required to identify the position and severity of damages accurately is different according to the damage scenarios and the percentage of uncertainty.