What to do with (large) imaginary frequencies for constrained minimum structures?
The first I noticed is you said:


  The two gradient tolerances have been loosened slightly in order to converge the calculations. By using the default gradient tolerances, the calculations simply do not, I have found (especially the maximum gradient).


At first sight, the structure of your molecule should not leads to serious problems, but the above is a bad signal, and turns the results unreliable (at least doubtful). The constrains must be considered in the calculation of vibrational frequencies for they make sense. 

Maybe no one can give you a foolproof steps series to follow. My suggestions are below. 

What to do/check


If comments are right, and Orca does not care about fixing atoms for
the calculation of the vibrational frequencies, then you need to
change the software for this task (because in this case results are not meaningful) (This also answer the question "Is it "safe" to disregard the frequencies that mainly involved fixed atoms?"). The simpler way to check this is by looking at the number of freedom degrees reported (or deduced form the output). I was referring directly to vibrational f.d. because low frequency modes may be treated like rotations. I doubt that Orca lacks of this basic functionality. In such case I bet that this is described their User Manual.
The step you are doing are reasonable. My guess is that maybe you are far from a true minimum, in some unstable geometry, and that is why you have problems with accuracy changes. You said: "Should I converge my geometry tighter to come closer to the true minimum? ". Yes, it is always advisable to get good frequencies. So I would use not only tighter geometry convergence criteria but the TIGHT convergence option for the SCF part.
In the same lines, you can try using modern semi-empirical methods and then make the fine tuning with DFT. Even using a simpler DFT functional.
I agree with comments that you are not using the best basis set out there, but I think that this is not a problem here.
If you want to discard weird behaviors you can try with another functional, but I do not think that this is the problem (because it is known that B3LYP performs well in compounds with the elements involved here).


In brief. I would check if Orca compute correctly the frequencies, if not: change the software. If yes: discard those results, try to converge the structure using cheaper methods (if you are lack of computational resources) from different initial conditions, store those results and start new geometry optimizations and vibrational frequencies calculations starting from those previous converged states. The last step using the tight convergence criteria for SCF and normal convergence criteria for geometries, and using the target basis set (for energy and its derivatives).

Good luck 

Note:


  Once I start computing reaction paths, I plan to try different
  functionals and basis sets to see if this give drastically different
  results.


Why do you expect that something like this would happen? You can save time by using a little large basis set and a better functional once you get the B3LYP/6-31g(d,p) results. It should not take too much time.



EDIT DUE TO COMMENTS

This answers generated some discussion about B3LYP or not B3LYP?

In that respect, my viewpoint is as follows:


B3LYP is not the most accurate functional (there should be no doubt about that). 
It is reasonably good for geometries in systems like yours. Errors in geometries can be larger due to small basis set than due to B3LYP, or at least of the same order.
While there are functional with better results for geometries, you also must take care of how large basis sets those functional needs to be a better option. I mean, maybe using 5$\zeta$ basis sets they are better, but not using a small basis set.  (and maybe not)
The above discussion can be very long, but in any case, you are using approximations that (in my opinion) are much stronger than the errors due to  B3LYP. For example, the permittivity value, or the geometry of the constained parts, etc.. 
Also, even if you could perform a FullCI/CBS you shoud notice that the simulated system is different to the real system, so extremely exacts results in the simulated system may no improve your knowledge about reality.
So, if the functionals are accurate enough, what it is important here is obtaining the corresponding results for the functional/basis set. The functionales can be more or less problematic. B3LYP should not be problematic for your system. My guess is what I said in the answer.
Because of the above, I think that B3LYP is reasonable and enough for your case.