Generic contour smoothing with shading

Update: Collecting the steps in a function:

ClearAll[smoothCP]
smoothCP = Module[{pr = PlotRange @ #, 
  nF =  Nearest[Join @@ Cases[Normal @ #, Line[a_] :> Transpose[#2 /@ Transpose@a], ∞]]},
  # /. GraphicsComplex[a_, b__] :> GraphicsComplex[ 
    If[Or @@ MemberQ @@@ Thread @ {pr, #}, #, First @ nF @ #] & /@ a, b] ] &;

Examples:

p1 = ContourPlot[x^2 + y^2 + RandomReal[{-0.1, 0.1}], {x, -2., 2.}, {y, -2., 2.}, 
     PlotPoints -> 6, Contours -> 4, ContourStyle -> Thick, ImageSize -> 400]; 
Row[{p1, smoothCP[p1, GaussianFilter[#, {5, 5}] &]}]

Row[{p1, smoothCP[p1, GaussianFilter[#, 2.5]&]}] 

Row[{p1, smoothCP[p1, GaussianFilter[#, {25, 25}] &]}] 

Row[{p1, smoothCP[p1, Downsample[#, 4]&]}] 

p2 = ContourPlot[x^2 + y^2 + RandomReal[{-0.1, 0.1}], {x, -2., 2.}, {y, -2., 2.}, 
   PlotPoints -> 6, Contours -> 4, ContourStyle -> Thick, 
   ContourShading -> False, ImageSize -> 400]; 
Row[{p2, smoothCP[p2, GaussianFilter[#, {5, 5}] & ]}] 

lcp = ListContourPlot[Table[Sin[j^2 + i] + 3 RandomReal[{-0.05, 0.05}],
   {i, -Pi, Pi, 0.1}, {j, -Pi, Pi, 0.1}], Contours -> 2, ImageSize -> 400, 
   ColorFunction -> "Rainbow"];
Row[{lcp, smoothCP[lcp, GaussianFilter[#, 10] &]}, Spacer[5]]

Original answer:

pl = ContourPlot[ x^2 + y^2 + RandomReal[{-0.1, 0.1}], {x, -1, 1}, {y, -1, 1}, 
   Contours -> 3, PlotPoints -> 3, ContourStyle -> Thick, ImageSize -> 400];
pnts =  Join @@ Cases[Normal[pl],  Line[a_] :> 
   Transpose[GaussianFilter[#, {5, 5}]& /@ Transpose[a]], ∞]; 
nf = Nearest[pnts];

pl2 = pl /. GraphicsComplex[a_, b__] :>  
  GraphicsComplex[If[Or @@ (Abs[#] == 1. & /@ #), #, nf[#][[1]]] & /@ a, b];
Row[{pl, pl2}]

Mean curvature flow and Laplacian smoothing can also be applied to the graph of a piecewise linear function for smoothing. Here a very brief and imperfect demonstration (using the function GraphDiffusionFlow from the cited post):

RandomSeed[20180506];
mesh2D = DiscretizeRegion[Rectangle[{-1, -1}, {1, 1}], 
   MaxCellMeasure -> {1 -> 0.05}];
vals = With[{pts2D = MeshCoordinates[mesh2D]},
   (pts2D^2).ConstantArray[1., 2] + 
    RandomVariate[NormalDistribution[0, 0.075], Length[pts2D]]
   ];
mesh3D = MeshRegion[Join[pts2D, Partition[vals, 1], 2], 
   MeshCells[mesh2D, 2, "Multicells" -> True]];
smoothedmesh3D = GraphDiffusionFlow[mesh3D, 20, 0.1, 0.8];
g = GraphicsRow[{mesh3D, smoothedmesh3D}, ImageSize -> Large]

This how the resulting height function looks as ContourPlot away from the boundary.

f = Interpolation[pts3D, InterpolationOrder -> 1];
fsmoothed = 
  Interpolation[MeshCoordinates[smoothedmesh3D], 
   InterpolationOrder -> 1];
r = .8;
g = GraphicsRow[{
   ContourPlot[f[u, v], {u, -r, r}, {v, -r, r}, Contours -> 20],
   ContourPlot[fsmoothed[u, v], {u, -r, r}, {v, -r, r}, Contours -> 20]
   }, ImageSize -> Large]

However, the current boundary conditions for the 3D-flow are not well-suited for this task:

r = 1.;
g = GraphicsRow[{
   ContourPlot[f[u, v], {u, -r, r}, {v, -r, r}, Contours -> 20],
   ContourPlot[fsmoothed[u, v], {u, -r, r}, {v, -r, r}, Contours -> 20]
   }, ImageSize -> Large]

There are certainly ways to fix that, though.

Not perfect but a good start, the idea is to use your method to transform boundaries and then use this transformation rules for points used by original output:

pl = ContourPlot[
  x^2 + y^2 + RandomReal[{-0.1, 0.1}], {x, -1, 1}, {y, -1, 1}, 
  Contours -> 3, PlotPoints -> 3, ContourShading -> True, 
  ContourStyle -> Thick];

rules = Flatten@
  Cases[Normal[pl], 
   Line[a__] :> 
    Thread[a -> 
      Transpose@{GaussianFilter[First /@ a, {5, 5}], 
        GaussianFilter[Last /@ a, {5, 5}]}], \[Infinity]];

pl /. rules