As you correctly point out, designing an optimal primer pair for 16S-rRNA sequencing is a tricky affair because even the less-variable regions are not same between different strains and species. Sambo et al (2018) have even developed a bioinformatics software for optimal design of primers for 16S-rRNA sequencing for multiple bacteria.
We propose here a computational method for optimizing the choice of
primer sets, based on multi-objective optimization, which
simultaneously: 1) maximizes efficiency and specificity of target
amplification; 2) maximizes the number of different bacterial 16S
sequences matched by at least one primer; 3) minimizes the differences
in the number of primers matching each bacterial 16S sequence. Our
algorithm can be applied to any desired amplicon length without
affecting computational performance.
There is diversity in the 16S-rRNA genes within the same species i.e. the different copies are not exact duplicates (Větrovský & Baldrian, 2013).
Větrovský & Baldrian (2013)
Interestingly, the different rRNA operons in E.coli have different promoters and are even differentially expressed during stress conditions (Kurylo et al., 2018).
However, Kitahara et al. (2012) found that 16S-rRNA genes isolated from soil samples, in place of the original E.coli gene, could support its growth. In other words E.coli is highly robust to mutations in its 16S-rRNA.
After counterselection, ∼200 clones of KT103-derivatives (carrying
pRB103 whose 16S rRNA gene was substituted with foreign genes) were
obtained, from which 33 nonredundant 16S rRNA genes (A01–H03) were
identified. Through multiple alignment of E. coli 16S rRNA and our
metagenomically retrieved 16S rRNA sequences, it was found that at
least 628 (40.7%) of the 1,542 nucleotides were variable, indicating
marked mutational robustness of the 16S rRNA. Strikingly, the
functional 16S rRNA sequences (except A10 and F02, which were 99.0%
identical to E. coli 16S rRNA) obtained in this study showed only
80.9–89.3% identity to E. coli 16S rRNA, which was well below the value reported thus far (Proteus vulgaris 16S rRNA, 94% identity to E.
coli 16S rRNA)
For your question
But the problem is that it has 7 distinct rRNA operons: rrnA, rrnB,
rrnC, rrnD, rrnE, rrnF, rrnG. Do you have a reference showing a
particular operon used for the 16S position notation?
I don't think any one of them is considered the reference. The complete E.coli 16S-rRNA sequence reported in NCBI seems to consider a "consensus" of different reported sequences:
,  contain updated sequence data for the original work by the
same laboratory . There were too many discrepancies between  and
,  to list each revision in our sites table. The sequence shown
is from . , ,  point to a number of cistron
heterogeneities. There is uncertainty, however, with regard to
assigning these various heterogeneities to specific cistrons. The RNA
method used by , ,  gives the average of all the cistrons
present in the cell . The heterogeneities are classified by their
relative proportions into major, minor and undetermined species. The
sequence shown corresponds to the major species. The heterogeneities
were annotated as variations in the sites table. It is not known which
of the residues 'c' (base 633) or 'a' (base 641) undergoes a deletion,
giving rise to the minor component 'atctg'.  suggests the existence
of one or two mutated cistrons among the known seven cistrons of
ribosomal RNA. With the exception of a single base deletion, this
sequence is identical to the current 16S rDNA sequence for the E.coli
The NCBI page also lists different polymorphisms and base modifications in the 16S-rRNA within and between different strains/"species".
NCBI has several partial sequences as well. When I checked for P.aeruginosa and B.subtilis, I could find only one or two complete sequences (rest were partial). Moreover, each entry indicated the strain that it was obtained from. Therefore I assume there is no single reference sequence.
I guess that people do consider the different variants while doing phylogenetic analysis (or simply just look for conserved "signature" sequences for a given species). I am sure that there are computational algorithms for doing the classification in an optimal manner (see Chatellier et al., 2014). Since, I don't have an expertise in this area I cannot say anything conclusively about the routine practice. In fact there is still ongoing research in improving the analysis (for example, Yang et al., 2016 and Sambo et al., 2018).