Cerebral blood flow (CBF) offers a quantitative and reliable measurement for brain activity and is increasingly used to study functional networks. However, current methods evaluate inter-regional relations mainly based on CBF temporal dynamics, which suffers from low signal-to-noise ratio and poor temporal resolution. Here we proposed a method to construct functional brain networks by estimating shape similarity (index by Jensen-Shannon divergence) in probability distributions of regional static CBF measured by arterial spin labeling perfusion imaging over a scanning period. Based on CBF data of 30 healthy participants from 10 visits, we found that the CBF networks exhibited non-trivial topological features (e.g., small-world organization, modular architecture, and hubs) and showed low-to-fair test-retest reliability and high between-subject consistency. We further found that interregional CBF similarities were depended on anatomical distance and differed between high- and lower-order subnetworks. Moreover, interregional CBF similarities within high-order subnetworks showed significantly lower reliability than those within low-order subnetworks. Finally, we showed that nodal degree of the CBF networks were related to regional sizes and CBF levels and spatially aligned with maps of the dopamine transporter and metabolic glutamate receptor 5 intensities, expression levels of genes primarily enriched in cholesterol-related pathways and endothelial cells, and meta-analytic activations related to memory, language, and executive function. Altogether, our proposed method provide a novel, relatively reliable, and neurobiologically meaningful means to study functional network organization of the human brain.